vector.h

矩阵运算的模板类

		else vec[i]=T();
	}
}


//template <typename T> Vector<T>::Vector(const Vector<T>& A) : vec(A.vec),n(A.n),type(A.type) {
//}


template <typename T> Vector<T>::~Vector() {
}



/////////////////////////////
// operator overloading
/////////////////////////////
//template <typename T> Vector<T>& Vector<T>::operator=(const Vector<T>& A) {
//	if(this == &A) return *this;
//	if(n != A.n) { vec.resize(A.n); n=A.n; }
//	vec=A.vec; type=A.type; return *this;
//}
template <typename T> Vector<T>& Vector<T>::operator=(const T& a) {
	for(int i=1; i<=n; ++i) vec[i] = a; return *this;
}


#ifdef ALLOW_ACCESS_BY_BRACKET
#if defined(_DEBUG)
template <typename T> T& Vector<T>::operator[](int i) {
	if(i<=0 || i>n){
		cerr << "ERROR: index invalid " << i << endl;
		assert(0);
	}
	return vec[i];
}


template <typename T> const T& Vector<T>::operator[](int i) const {
	if(i<=0 || i>n){
		cerr << "ERROR: index invalid " << i << endl;
		assert(0);
	}
	return vec[i];
}


#else // !_DEBUG
template <typename T> inline T& Vector<T>::operator[](int i) {
	return vec[i];
}


template <typename T> inline const T& Vector<T>::operator[](int i) const {
	return vec[i];
}
#endif // _DEBUG
#endif // ALLOW_ACCESS_BY_BRACKET


#if defined(_DEBUG)
template <typename T> const T& Vector<T>::operator()(int i) const {
	if(i<=0 || i>n){
		cerr << "ERROR: index invalid " << i << endl;
		assert(0);
	}
	return vec[i];
}


template <typename T> T& Vector<T>::operator()(int i) {
	if(i<=0 || i>n){
		cerr << "ERROR: index invalid " << i << endl;
		assert(0);
	}
	return vec[i];
}


#else // !_DEBUG
template <typename T> inline const T& Vector<T>::operator()(int i) const {
	return vec[i];
}


template <typename T> inline T& Vector<T>::operator()(int i) {
	return vec[i];
}

#endif // _DEBUG


template <typename T> inline Vector<T> Vector<T>::operator+() {
	return (*this);
}


template <typename T> Vector<T> Vector<T>::operator-() {
	Vector<T> C=(*this); C*=(-1); return C;
}


template <typename T> Vector<T>& Vector<T>::operator+=(const Vector<T>& B) {
	assert(n==B.n);
	assert(type==B.type);
	for(int i=1;i<=n;i++) vec[i]+=B.vec[i]; return (*this);
}


template <typename T> Vector<T>& Vector<T>::operator-=(const Vector<T>& B) {
	assert(n==B.n);
	assert(type==B.type);
	for(int i=1;i<=n;i++) vec[i]-=B.vec[i]; return (*this);
}


template <typename T> Vector<T>& Vector<T>::operator+=(const double b) {
	for(int i=1;i<=n;i++) vec[i]+=b; return (*this);
}


template <typename T> Vector<T>& Vector<T>::operator-=(const double b) {
	for(int i=1;i<=n;i++) vec[i]-=b; return (*this);
}


template <typename T> Vector<T>& Vector<T>::operator*=(const double b) {
	for(int i=1;i<=n;i++) vec[i]*=b; return (*this);
}


template <typename T> Vector<T>& Vector<T>::operator/=(const double b) {
	for(int i=1;i<=n;i++) vec[i]/=b; return (*this);
}


#ifndef DISABLE_COMPLEX
template <typename T> Vector<T>& Vector<T>::operator+=(const complex<T> b) {
	for(int i=1;i<=n;i++) vec[i]+=b; return (*this);
}


template <typename T> Vector<T>& Vector<T>::operator-=(const complex<T> b) {
	for(int i=1;i<=n;i++) vec[i]-=b; return (*this);
}


template <typename T> Vector<T>& Vector<T>::operator*=(const complex<T> b) {
	for(int i=1;i<=n;i++) vec[i]*=b; return (*this);
}


template <typename T> Vector<T>& Vector<T>::operator/=(const complex<T> b) {
	for(int i=1;i<=n;i++) vec[i]/=b; return (*this);
}
#endif


/////////////////////////////
// member functions
/////////////////////////////
template <typename T> inline int Vector<T>::Size() const { 
	return n;
}


template <typename T> inline VectorType Vector<T>::Type() const { 
	return type;
}


template <typename T> inline int Vector<T>::SetType(VectorType type_) {
	type=type_;
	return 1;
}


template <typename T> Vector<T> Vector<T>::t() const {
	Vector<T> C(Size(),(Type()==ROW_VECTOR)?COL_VECTOR:ROW_VECTOR,T());
	for(int i=1;i<=Size();i++) C[i]=((*this)[i]);
	return C;
}


template <typename T> Vector<T> Vector<T>::ct() const {
	Vector<T> C(Size(),(Type()==ROW_VECTOR)?COL_VECTOR:ROW_VECTOR,T());
	for(int i=1;i<=Size();i++) C[i]=conj((*this)[i]);
	return C;
}


template <typename T> inline int Vector<T>::Resize(int new_n) { 
	vec.resize(new_n+1);
	n=new_n;
	return 1;
}


template <typename T> int Vector<T>::Add(const T& a){ 
	int n0=n; 
	Resize(n0+1); 
	vec[n0+1]=a; 
	n=n0+1; 
	return n; 
}


template <typename T> int Vector<T>::Erase(int c){ 
	vec.erase(vec.begin()+c); 
	n--;
	return 1;
}


template <typename T> inline void Vector<T>::Clear(){ 
	vec.clear(); 
}


template <typename T> inline bool Vector<T>::Empty() const{ 
	return (Size()==0);
}


template <typename T> T Vector<T>::Sum_() const {
	T sum=0; 
	for(int i=1;i<=n;i++) sum+=vec[i]; 
	return sum; 
}


template <typename T> T Min(const Vector<T>& A,int* idx=0){
	return A.Min_(idx); 
}


template <typename T> T Vector<T>::Max_(int* idx) const {
	T maxA=-mtl_numeric_limits<T>::max();
	if(idx) *idx=0;
	for(int i=1;i<=n;i++){
		if(maxA<vec[i]){
			maxA=vec[i];
			if(idx) *idx=i;
		}
	}
	return T(maxA);
}


template <typename T> T Vector<T>::Min_(int* idx) const {
	T minA=mtl_numeric_limits<T>::max();
	if(idx) *idx=0;
	for(int i=1;i<=n;i++){
		if(minA>vec[i]){
			minA=vec[i];
			if(idx) *idx=i;
		}
	}
	return T(minA);
}


template <typename T> int Vector<T>::Print(ostream& s, const char *name) {
	if(name != 0 && name[0] != 0) s << name << " ";
	s << 1 << " " << n << endl;
	//s << (*this) << endl;
	mtl_ostream(s,(*this));
	return 1;
}


template <typename T> int Vector<T>::Print(ostream& s, string& name) {
	return Print(s,name.c_str());
}


template <typename T> int Vector<T>::Read(istream& s, const char *name) {
	string dummy;
	int itmp, nn;
	s >> dummy;
	if(name != 0 && name[0] != 0){
		assert(name==dummy);
	}
	s >> itmp >> nn;
	assert(itmp==1);
	Resize(nn);
	s >> (*this);
	return 1;
}


template <typename T> int Vector<T>::Read(istream& s, string& name) {
	return Read(s,name.c_str());
}


template <typename T> Vector<T> Vector<T>::SubVector(int n1, int n2){
	assert(n2>=n1);
	Vector<T> C(n2-n1+1); for(int i=1;i<=n2-n1+1;i++) C[i]=vec[i-1+n1];
	return C;
}


template <typename T> int Vector<T>::Save(const string& fname){
	FILE* fp=fopen(fname.c_str(),"wb");
	if(!fp) return 0;
	fprintf(fp,"FMAT 1 %d\n",n);
	for(int j=1;j<=n;j++){
		if((j-1)%10!=0) fprintf(fp," ");
		fprintf(fp,"%f",vec[j]);
		if(j%10==0) fprintf(fp,"\n");
	}
	if(n%10!=0) fprintf(fp,"\n");
	fclose(fp);
	return n;
}


template <typename T> int Vector<T>::SaveMatlab(const string& fname){
	FILE* fp=fopen(fname.c_str(),"wb");
	if(!fp) return 0;
	for(int j=1;j<=n;j++){
		fprintf(fp,"%f ",vec[j]);
	}
	fclose(fp);
	return n;
}


template <typename T> int Vector<T>::IsOne(double tolerance) const {
	for (int i=1; i<=n; i++) if(Abs(vec[i]-1)>tolerance) return 0;
	return 1;
}


template <typename T> int Vector<T>::IsZero(double tolerance) const {
	for (int i=1; i<=n; i++) if(Abs(vec[i])>tolerance) return 0;
	return 1;
}


template <typename T> int Vector<T>::IsUnit(double tolerance) const {
	double norm=Norm(*this);
	if(Abs(norm-1)>tolerance) return 0;
	return 1;
}



/////////////////////////////
// basic friend functions
/////////////////////////////
template <typename T> Vector<T> mtl_plus(const Vector<T>& A, const Vector<T>& B) {
	assert(A.n==B.n);
	assert(A.type==B.type);
	Vector<T> C=A; C+=B; return C;
}


template <typename T> Vector<T> mtl_minus(const Vector<T>& A, const Vector<T>& B) {
	assert(A.n==B.n);
	assert(A.type==B.type);
	Vector<T> C=A; C-=B; return C;
}


template <typename T> Vector<T> mtl_plus(const Vector<T>& A, const double b) {
	Vector<T> C=A; C+=b; return C;
}


template <typename T> Vector<T> mtl_minus(const Vector<T>& A, const double b) {
	Vector<T> C=A; C-=b; return C;
}


template <typename T> Vector<T> mtl_times(const Vector<T>& A, const double b) {
	Vector<T> C=A; C*=b; return C;
}


template <typename T> Vector<T> mtl_divide(const Vector<T>& A, const double b) {
	Vector<T> C=A; C/=b; return C;
}


template <typename T> Vector<T> mtl_plus(const double b, const Vector<T>& A) {
	Vector<T> C=A; C+=b; return C;
}


template <typename T> Vector<T> mtl_minus(const double b, const Vector<T>& A) {
	Vector<T> C=-A; C+=b; return C;
}


template <typename T> Vector<T> mtl_times(const double b, const Vector<T>& A) {
	Vector<T> C=A; C*=b; return C;
}


template <typename T> Vector<T> mtl_divide(const double b, const Vector<T>& A) {
	Vector<T> C(A.n); for(int i=1;i<=A.n;i++) C.vec[i]=b/A.vec[i]; return C;
}


#ifndef DISABLE_COMPLEX
template <typename T> Vector<T> mtl_plus(const Vector<T>& A, const complex<T> b) {
	Vector<T> C=A; C+=b; return C;
}


template <typename T> Vector<T> mtl_minus(const Vector<T>& A, const complex<T> b) {
	Vector<T> C=A; C-=b; return C;
}


template <typename T> Vector<T> mtl_times(const Vector<T>& A, const complex<T> b) {
	Vector<T> C=A; C*=b; return C;
}


template <typename T> Vector<T> mtl_divide(const Vector<T>& A, const complex<T> b) {
	Vector<T> C=A; C/=b; return C;
}


template <typename T> Vector<T> mtl_plus(const complex<T> b, const Vector<T>& A) {
	Vector<T> C=A; C+=b; return C;
}


template <typename T> Vector<T> mtl_minus(const complex<T> b, const Vector<T>& A) {
	Vector<T> C=-A; C+=b; return C;
}


template <typename T> Vector<T> mtl_times(const complex<T> b, const Vector<T>& A) {
	Vector<T> C=A; C*=b; return C;
}


template <typename T> Vector<T> mtl_divide(const complex<T> b, const Vector<T>& A) {
	Vector<T> C(A.n); for(int i=1;i<=A.n;i++) C.vec[i]=b/A.vec[i]; return C;
}
#endif


template <typename T> Vector<T> mtl_times(const Vector<T>& A, const Vector<T>& B){
	assert(A.Size()==B.Size());
	assert(A.type==B.type);