*** empty log message ***

mat.cc

@@ -492,6 +492,15 @@ return r;
 }
 
 
+//randomize
+template<>
+void NRMat<double>::randomize(const double &x)
+{
+for(int i=0; i<nn; ++i)
+	for(int j=0; j<mm; ++j)
+		(*this)(i,j) = x*(2.*random()/(1.+RAND_MAX) -1.);
+}
+
 
 
 // Mat *= a

mat.h

@@ -58,6 +58,7 @@ public:
 	inline int getcount() const {return count?*count:0;}
 	NRMat & operator=(const NRMat &rhs);  //assignment
         void clear() {LA_traits<T>::clear((*this)[0],nn*mm);}; //zero out
+	void randomize(const T &x); //fill with random numbers
 	NRMat & operator=(const T &a);    //assign a to diagonal
 	NRMat & operator|=(const NRMat &rhs); //assignment to a new copy
 	NRMat & operator+=(const T &a);   //add diagonal

smat.cc

@@ -132,6 +132,12 @@ const T NRSMat<T>::trace() const
 	return tmp;
 }
 
+template<>
+void NRSMat<double>::randomize(const double &x)
+{
+for(int i=0; i<NN2; ++i) v[i] = x*(2.*random()/(1.+RAND_MAX) -1.);
+}
+
 
 
 // write matrix to the file with specific format

smat.h

@@ -43,6 +43,7 @@ public:
 	NRSMat & operator|=(const NRSMat &rhs);	//assignment to a new copy
 	NRSMat & operator=(const NRSMat &rhs);	//assignment
 	void clear() {LA_traits<T>::clear(v,NN2);}; //zero out
+        void randomize(const T&x);
 	NRSMat & operator=(const T &a);		//assign a to diagonal
         const bool operator!=(const NRSMat &rhs) const {if(nn!=rhs.nn) return 1; return LA_traits<T>::gencmp(v,rhs.v,NN2);} //memcmp for scalars else elementwise
         const bool operator==(const NRSMat &rhs) const {return !(*this != rhs);};

t.cc

@@ -1163,7 +1163,7 @@ cin >>a>>b;
 cout <<a*b;
 }
 
-if(1)
+if(0)
 {
 NRMat<double> a,b;
 cin >>a >>b;
@@ -1171,6 +1171,110 @@ cout <<a.oplus(b);
 cout <<a.otimes(b);
 }
 
+//test of general square matrix eigenvector derivatives
+if(1)
+{
+const bool biorthonormalize=1;
+
+NRMat<double> a;
+cin >>a;
+if(a.nrows()!=a.ncols()) laerror("must be square matrix");
+int n=a.nrows();
+
+NRMat<double>vl(n,n),vr(n,n);
+NRVec<double> wr(n),wi(n);
+NRMat<double> awork(a);
+gdiagonalize(awork,wr,wi,&vl,&vr,1,0,1,biorthonormalize);
+for(int i=0; i<n; ++i) if(wi[i]) laerror("try another matrix with real eigenvalues");
+
+NRMat<double> eival(0.,n,n);
+eival.diagonalset(wr);
+
+cout <<"test biorthonormality "<< (vl.transpose() * vr - 1.).norm()<<endl;
+
+NRMat<double> hlp;
+hlp = a *vr - vr * eival;
+cout  <<"test right eivectors "<<hlp.norm()<<endl;
+
+hlp= vl.transpose() * a - eival * vl.transpose();
+cout  <<"test left eivectors "<<hlp.norm()<<endl;
+
+hlp =  vl.transpose() * a * vr - eival;
+cout  <<"test eigenvalues "<<hlp.norm()<<endl;
+
+NRMat<double> ader(n,n);
+for(int i=0; i<n; ++i)
+	for(int j=0; j<n; ++j)
+		ader(i,j) = 2.*random()/(1.+RAND_MAX) -1.;
+
+cout <<"eigenvalues\n"<<wr<<endl;
+
+
+
+//compute numerical derivatives of eigenvectors
+double  h=1e-5;
+awork = a+ ader*h;
+NRMat<double> vlx(n,n),vrx(n,n);
+NRVec<double> wrx(n),wix(n);
+gdiagonalize(awork,wrx,wix,&vlx,&vrx,1,0,1,biorthonormalize);
+for(int i=0; i<n; ++i) if(wix[i]) laerror("try another matrix with real eigenvalues");
+
+awork = a - ader*h;
+NRMat<double> vly(n,n),vry(n,n);
+NRVec<double> wry(n),wiy(n);
+gdiagonalize(awork,wry,wiy,&vly,&vry,1,0,1,biorthonormalize);
+for(int i=0; i<n; ++i) if(wiy[i]) laerror("try another matrix with real eigenvalues");
+
+NRMat<double> vld,vrd;
+NRVec<double> wrd;
+vld = (vlx-vly) * (.5/h);
+vrd = (vrx-vry) * (.5/h);
+wrd = (wrx-wry) * (.5/h);
+
+NRMat<double> vlg,vrg;
+NRVec<double> wrg(n);
+
+//compute analytic derivative
+NRMat<double> tmp(n,n);
+tmp.gemm(0.,vl,'t', ader * vr,'n',1.);
+hlp |= tmp;
+cout <<" C~+ VH C = \n"<<tmp<<endl;
+
+tmp.diagonalof(wrg);
+for(int i=0; i<n; ++i)
+        for(int j=0; j<n; ++j)
+		if(i!=j) tmp(i,j) /= (wr[j] - wr[i]); else  tmp(i,j) = 0.;
+cout <<" old X matrix \n"<<tmp<<endl;
+
+NRMat<double> Y = tmp;
+NRMat<double> S = vr.transpose() * vr;
+NRMat<double> tmp2 = S * tmp;
+for(int i=0; i<n; ++i) Y(i,i) -= tmp2(i,i);
+
+cout <<"Y matrix \n"<< Y;
+
+NRMat<double> numX = inverse(vr) * vrd;
+cout <<" numerical X matrix \n"<< numX;
+cout <<" numerical X matrix test = "<< (vr * numX - vrd).norm()<<endl;
+
+vrg = vr * Y;
+
+//and compare
+cout <<"eigenvalue numerical derivative\n"<<wrd<<endl;
+cout <<"eigenvalue analytic derivative\n"<<wrg<<endl;
+cout <<"eigenvalue derivative error = "<<(wrd-wrg).norm()<<endl;
+
+//and for right eigenvectors
+cout <<"right eigenvector  numerical derivative\n"<<vrd<<endl;
+cout <<"right eigenvector  analytic derivative\n"<<vrg<<endl;
+cout <<"right eigenvector derivative error = "<<(vrd-vrg).norm()<<endl;
+
+
+
+
+
+
+}
 
 
 }

vec.cc

@@ -162,6 +162,12 @@ return nn<rhs.nn;
 }
 
 
+template<>
+void NRVec<double>::randomize(const double &x)
+{
+for(int i=0; i<nn; ++i) v[i] = x*(2.*random()/(1.+RAND_MAX) -1.);
+}
+
 
 // axpy call for T = double (not strided)
 template<>

vec.h

@@ -72,6 +72,7 @@ public:
 	NRVec & operator=(const NRVec &rhs);
 	NRVec & operator=(const T &a);  //assign a to every element
 	void clear() {LA_traits<T>::clear(v,nn);}; //zero out
+        void randomize(const T &x);
 	NRVec & operator|=(const NRVec &rhs);
 	const bool operator!=(const NRVec &rhs) const {if(nn!=rhs.nn) return 1; return LA_traits<T>::gencmp(v,rhs.v,nn);} //memcmp for scalars else elementwise
 	const bool operator==(const NRVec &rhs) const {return !(*this != rhs);};