*** empty log message ***

fourindex.h

@@ -622,7 +622,7 @@ public:
 };
 
 template<class T, class I> 
-fourindex_dense<twoelectronrealmullikan,T,I>::fourindex_dense<twoelectronrealmullikan,T,I>(const fourindex<I,T> &rhs) : NRSMat<T>((T)0,rhs.size()*(rhs.size()+1)/2)
+fourindex_dense<twoelectronrealmullikan,T,I>::fourindex_dense(const fourindex<I,T> &rhs) : NRSMat<T>((T)0,rhs.size()*(rhs.size()+1)/2)
 {
 if(rhs.getsymmetry() != twoelectronrealmullikan ) laerror("fourindex_dense symmetry mismatch");
 typename fourindex<I,T>::iterator p;
@@ -635,7 +635,7 @@ for(p=rhs.begin(); p!= rhs.end(); ++p) (*this)(p->index.indiv.i,p->index.indiv.j
 }
 
 template<class T, class I>
-fourindex_dense<twoelectronrealmullikan,T,I>::fourindex_dense<twoelectronrealmullikan,T,I>(const fourindex_ext<I,T> &rhs) : NRSMat<T>((T)0,rhs.size()*(rhs.size()+1)/2)
+fourindex_dense<twoelectronrealmullikan,T,I>::fourindex_dense(const fourindex_ext<I,T> &rhs) : NRSMat<T>((T)0,rhs.size()*(rhs.size()+1)/2)
 {
 if(rhs.getsymmetry() != twoelectronrealmullikan ) laerror("fourindex_dense symmetry mismatch");
 typename fourindex_ext<I,T>::iterator p;

mat.cc

@@ -11,17 +11,6 @@ extern ssize_t write(int, const void *, size_t);
 // TODO :
 //
 
-//////////////////////////////////////////////////////////////////////////////
-//// forced instantization in the corresponding object file
-template class NRMat<double>;
-template class NRMat<complex<double> >;
-template class NRMat<int>;
-template class NRMat<short>;
-template class NRMat<char>;
-template class NRMat<unsigned char>;
-template class NRMat<unsigned int>;
-template class NRMat<unsigned long>;
-
 
 /*
  *  Templates first, specializations for BLAS next
@@ -1113,3 +1102,14 @@ istream& operator>>(istream  &s, NRMat<T> &x)
 
 
 
+//////////////////////////////////////////////////////////////////////////////
+//// forced instantization in the corresponding object file
+template class NRMat<double>;
+template class NRMat<complex<double> >;
+template class NRMat<int>;
+template class NRMat<short>;
+template class NRMat<char>;
+template class NRMat<unsigned char>;
+template class NRMat<unsigned int>;
+template class NRMat<unsigned long>;
+

mat.h

@@ -560,6 +560,7 @@ public:
         NRMat_from1(): NRMat<T>() {};
         explicit NRMat_from1(const int n): NRMat<T>(n) {};
         NRMat_from1(const NRMat<T> &rhs): NRMat<T>(rhs) {}; //be able to convert the parent class  transparently to this
+        NRMat_from1(const int n, const int m): NRMat<T>(n,m) {};
         NRMat_from1(const T &a, const int n, const int m): NRMat<T>(a,n,m) {};
         NRMat_from1(const T *a, const int n, const int m): NRMat<T>(a,n,m) {};
 

matexp.h

@@ -109,7 +109,7 @@ return z;
 
 //general BCH expansion (can be written more efficiently in a specialization for matrices)
 template<class T>
-const T BCHexpansion (const T &h, const T &t, const int n, const bool verbose=1)\
+const T BCHexpansion (const T &h, const T &t, const int n, const bool verbose=0)\
 {
 T result=h;
 double factor=1.;
@@ -118,7 +118,7 @@ for(int i=1; i<=n; ++i)
 	{
 	factor/=i;
 	z= z*t-t*z;
-	if(verbose) cerr << "BCH contribution at order "<<i<<" : "<<z.norm()<<endl;
+	if(verbose) cerr << "BCH contribution at order "<<i<<" : "<<z.norm()*factor<<endl;
 	result+= z*factor; 
 	}
 return result;
@@ -184,12 +184,12 @@ static double exptaylor[]={
 4.1103176233121648441e-19,
 0.};
 double mnorm= LA_traits<T>::norm(x);
-power=nextpow2(mnorm);
+power=nextpow2(mnorm); 
 double scale=exp(-log(2.)*power);
 
 
 //find how long taylor expansion will be necessary
-const double precision=1e-16;
+const double precision=1e-14; //decreasing brings nothing
 double s,t;
 s=mnorm*scale;
 int n=0;
@@ -202,6 +202,7 @@ while(t*exptaylor[n]>precision);//taylor 0 will terminate in any case
 
 
 
+
 int i; //adjust the coefficients in order to avoid scaling the argument
 NRVec<typename LA_traits<T>::elementtype> taylor2(n+1);
 for(i=0,t=1.;i<=n;i++)
@@ -214,8 +215,9 @@ return taylor2;
 
 
 
+//it seems that we do not gain anything by polynom vs polynom0, check the m-optimization!
 template<class T>
-const T exp(const T &x, const bool simple=false)
+const T exp(const T &x, const bool horner=true)
 {
 int power;
 
@@ -223,7 +225,7 @@ int power;
 NRVec<typename LA_traits<T>::elementtype> taylor2=exp_aux(x,power);
 
 
-T r= simple?polynom0(x,taylor2):polynom(x,taylor2); 
+T r= horner?polynom0(x,taylor2):polynom(x,taylor2); 
 //for accuracy summing from the smallest terms up would be better, but this is more efficient for matrices
 
 //power the result back

nonclass.cc

@@ -617,7 +617,9 @@ NRMat<double> matrixfunction(NRMat<double> a, complex<double>
 	int n = a.nrows();
 	NRMat< complex<double> > u(n, n), v(n, n);
 	NRVec< complex<double> > w(n);
+/*
 NRMat<complex<double> > a0=complexify(a);
+*/
 	gdiagonalize(a, w, &u, &v);//a gets destroyed, eigenvectors are rows
 	NRVec< complex<double> > z = diagofproduct(u, v, 1, 1);
 /*

nonclass.h

@@ -70,6 +70,7 @@ extern const  NRVec<T> diagofproduct(const NRMat<T> &a, const NRMat<T> &b,\
 		bool trb=0, bool conjb=0); \
 extern T trace2(const NRMat<T> &a, const NRMat<T> &b, bool trb=0); \
 extern T trace2(const NRSMat<T> &a, const NRSMat<T> &b, const bool diagscaled=0);\
+extern T trace2(const NRSMat<T> &a, const NRMat<T> &b, const bool diagscaled=0);\
 extern void linear_solve(NRMat<T> &a, NRMat<T> *b, double *det=0,int n=0); \
 extern void linear_solve(NRSMat<T> &a, NRMat<T> *b, double *det=0, int n=0); \
 extern void linear_solve(NRMat<T> &a, NRVec<T> &b, double *det=0, int n=0); \
@@ -108,6 +109,7 @@ inline NRMat<double>  realsqrt(const NRMat<double> &a) { return realmatrixfuncti
 inline NRMat<double>  realsqrtinv(const NRMat<double> &a) { return realmatrixfunction(a,&sqrtinv); }
 inline NRMat<double>  log(const NRSMat<double> &a) { return matrixfunction(a,&log); }
 extern NRMat<double> log(const NRMat<double> &a);
+extern NRMat<double> exp0(const NRMat<double> &a);
 
 
 extern const NRMat<double> realpart(const NRMat< complex<double> >&);

smat.cc

@@ -12,16 +12,6 @@ extern ssize_t write(int, const void *, size_t);
 // specialize unary minus
 
 
-//////////////////////////////////////////////////////////////////////////////
-////// forced instantization in the corresponding object file
-template class NRSMat<double>;
-template class NRSMat< complex<double> >;
-template class NRSMat<int>;
-template class NRSMat<short>;
-template class NRSMat<char>;
-template class NRSMat<unsigned int>;
-template class NRSMat<unsigned long>;
-
 
 
 /*
@@ -411,3 +401,15 @@ INSTANTIZE(char)
 INSTANTIZE(unsigned int)
 INSTANTIZE(unsigned long)
 
+
+
+//////////////////////////////////////////////////////////////////////////////
+////// forced instantization in the corresponding object file
+template class NRSMat<double>;
+template class NRSMat< complex<double> >;
+template class NRSMat<int>;
+template class NRSMat<short>;
+template class NRSMat<char>;
+template class NRSMat<unsigned int>;
+template class NRSMat<unsigned long>;
+

sparsemat.cc

@@ -17,11 +17,6 @@ static inline void SWAP(T &a, T &b)
 
 
 
-//////////////////////////////////////////////////////////////////////////////
-//// forced instantization in the corresponding object file
-template class SparseMat<double>;
-template class SparseMat<complex<double> >;
-
 
 
 export template <class T>
@@ -1327,13 +1322,13 @@ template unsigned int SparseMat<T>::length() const; \
 template void SparseMat<T>::deletelist(); \
 template void SparseMat<T>::simplify(); \
 template void SparseMat<T>::copylist(const matel<T> *l); \
-template void SparseMat<T>::add(const SPMatindex n, const SPMatindex m, const T elem); \
 template SparseMat<T> & SparseMat<T>::operator=(const SparseMat<T> &rhs); \
 template SparseMat<T> & SparseMat<T>::operator+=(const SparseMat<T> &rhs); \
 template SparseMat<T> & SparseMat<T>::operator-=(const SparseMat<T> &rhs); \
 template SparseMat<T>::SparseMat(const NRMat<T> &rhs); \
 template SparseMat<T>::SparseMat(const NRSMat<T> &rhs); \
 template void SparseMat<T>::transposeme(); \
+template const T* SparseMat<T>::diagonalof(NRVec<T> &r, const bool divide, bool cache) const; \
 template SparseMat<T> & SparseMat<T>::operator*=(const T a); \
 template void SparseMat<T>::setunsymmetric(); \
 template SparseMat<T> & SparseMat<T>::operator=(const T a); \
@@ -1362,3 +1357,8 @@ INSTANTIZE(double)
 INSTANTIZE(complex<double>) //some functions are not OK for hermitean matrices, needs a revision!!!
 
 
+//////////////////////////////////////////////////////////////////////////////
+//// forced instantization in the corresponding object file
+template class SparseMat<double>;
+template class SparseMat<complex<double> >;
+

t.cc

@@ -52,7 +52,8 @@ NRVec<double> x(1.,10);
 NRVec<double> y(2.,10);
 NRVec<double> z(-2.,10);
 
-cout.setf(ios::fixed);
+cout.setf(ios::scientific);
+//cc:cout.setf(ios::fixed);
 cout.precision(12);
 
 
@@ -646,7 +647,6 @@ cout <<"test orthogonality\n" << vl.transpose() * vr;
 
 if(0)
 {
-/*
 int n;
 cin>>n;
 NRMat<double> a(n,n);
@@ -655,39 +655,49 @@ for(int i=0;i<n;++i) for(int j=0;j<i;++j)
         a(i,j)= random()/(1.+RAND_MAX);
         a(j,i)= -a(i,j);
         }
-NRMat<double> b=exp(a);
+cout <<"a matrix \n"<<a;
-cout <<a;
+cout<<"EXP\n";
-*/
+NRMat<double> b=exp0(a);
-NRMat<double> a,b;
+cout <<"b=exp(a) matrix\n"<<b;
-cin >>b;
-int n=b.nrows();
-cout <<"difference from identity = "<<b.norm(1.)<<endl;
 
-NRMat<double> x(0.,n,n),x0;
+cout <<"LOG\n";
-	double r;
-int i=0;
-do
-	{
-	x0=x;
-	NRMat<double> y=exp(x*-.5);
-	x+= y*b*y; 
-	x-= 1.;
-	x=(x-x.transpose())*.5;
-	cout <<"matrix x\n"<<x;
-	cout <<"iter "<<i <<" residue "<< (r=(exp(x)-b).norm())<<endl;
-	cout <<"iter "<<i <<" conv "<<(r=(x-x0).norm())<<endl;
-	++i;
-	} while(abs(r)>1e-10);
-cout <<"result\n"<<x<<endl;
-cout <<"exp(result)"<<exp(x)<<endl;
 NRMat<double> c=log(b); //matrixfunction(a,&mycident,1);
-cout <<c;
+cout <<"c=log(exp(a))\n"<<c <<"error: "<<(c-a).norm()<<endl;
-NRMat<double> d=exp(c);
+cout <<"EXP-MY\n";
-cout <<"exp(log(x))\n"<<d;
+NRMat<double> e=exp(c);
-cout<<(d-b).norm()<<endl;
+cout <<"e=exp(c)\n"<<e;
+cout<<"error2 = "<<(e-b).norm()<<endl;
+}
+
+if(1)
+{
+int n;
+double f;
+cin>>n >>f ;
+NRMat<double> a(n,n);
+for(int i=0;i<n;++i) for(int j=0;j<n;++j)
+        {
+        a(i,j)= f*random()/(1.+RAND_MAX);
+        }
+//cout <<"a matrix \n"<<a;
+//cout<<"EXP\n";
+double t=clock()/((double) (CLOCKS_PER_SEC));
+NRMat<double> b=exp0(a);
+cout <<"exp0 took "<<clock()/((double) (CLOCKS_PER_SEC))-t<<endl;
+//cout <<"b=exp0(a) matrix\n"<<b;
+t=clock()/((double) (CLOCKS_PER_SEC));
+NRMat<double> c=exp(a,true);
+cout <<" horner exp  took "<<clock()/((double) (CLOCKS_PER_SEC))-t<<endl;
+//cout <<"exp(a)\n"<<c;
+cout<<"error1 = "<<(c-b).norm()/b.norm()<<endl;
+t=clock()/((double) (CLOCKS_PER_SEC));
+c=exp(a,false);
+cout <<" tricky exp  took "<<clock()/((double) (CLOCKS_PER_SEC))-t<<endl;
+cout<<"error2 = "<<(c-b).norm()/b.norm()<<endl;
 
 }
 
+
 if(0)
 {
 int n;
@@ -1115,7 +1125,7 @@ cout <<b*b;
 cout <<(b*b-a).norm();
 }
 
-if(1)
+if(0)
 {
 NRSMat<double> a;
 NRMat<double> b;

vec.cc

@@ -25,20 +25,13 @@ template void NRVec<T>::get(int fd, bool dim, bool transp); \
 INSTANTIZE(double)
 INSTANTIZE(complex<double>)
 INSTANTIZE(int)
-INSTANTIZE(unsigned int)
 INSTANTIZE(short)
 INSTANTIZE(unsigned short)
 INSTANTIZE(char)
 INSTANTIZE(unsigned char)
+INSTANTIZE(unsigned int)
 INSTANTIZE(unsigned long)
 
-template class NRVec<double>;
-template class NRVec<complex<double> >;
-template class NRVec<char>;
-template class NRVec<short>;
-template class NRVec<int>;
-template class NRVec<unsigned long>;
-
 
 
 
@@ -288,6 +281,8 @@ template<>
 NRVec<char> &  NRVec<char>::normalize() {laerror("normalize() impossible for integer types"); return *this;}
 template<>
 NRVec<unsigned long> &  NRVec<unsigned long>::normalize() {laerror("normalize() impossible for integer types"); return *this;}
+template<>
+NRVec<unsigned int> &  NRVec<unsigned int>::normalize() {laerror("normalize() impossible for integer types"); return *this;}
 
 
 template<>
@@ -314,6 +309,15 @@ void NRVec<unsigned long>::gemv(const unsigned long beta,
 laerror("not yet implemented");
 }
 
+template<>
+void NRVec<unsigned int>::gemv(const unsigned int beta,
+                const NRSMat<unsigned int> &A, const char trans,
+                const unsigned int alpha, const NRVec &x)
+{
+laerror("not yet implemented");
+}
+
+
 template<>
 void NRVec<unsigned char>::gemv(const unsigned char beta,
                 const NRSMat<unsigned char> &A, const char trans,
@@ -355,6 +359,15 @@ void NRVec<unsigned long>::gemv(const unsigned long beta,
 laerror("not yet implemented");
 }
 
+template<>
+void NRVec<unsigned int>::gemv(const unsigned int beta,
+                const NRMat<unsigned int> &A, const char trans,
+                const unsigned int alpha, const NRVec &x)
+{
+laerror("not yet implemented");
+}
+
+
 template<>
 void NRVec<unsigned char>::gemv(const unsigned char beta,
                 const NRMat<unsigned char> &A, const char trans,
@@ -407,6 +420,16 @@ void NRVec<unsigned long>::gemv(const unsigned long beta,
 laerror("not yet implemented");
 }
 
+template<>
+void NRVec<unsigned int>::gemv(const unsigned int beta,
+                const SparseMat<unsigned int> &A, const char trans,
+                const unsigned int alpha, const NRVec &x, bool s)
+{
+laerror("not yet implemented");
+}
+
+
+
 template<>
 void NRVec<unsigned char>::gemv(const unsigned char beta,
                 const SparseMat<unsigned char> &A, const char trans,
@@ -501,3 +524,17 @@ if(to == -1) to=nn-1;
 if(direction) return memqsort<1,NRVec<T>,int,int>(*this,perm,from,to);
 else return memqsort<0,NRVec<T>,int,int>(*this,perm,from,to);
 }
+
+
+//////////////////////////////////////////////////////////////////////////////
+//// forced instantization in the corespoding object file
+
+
+template class NRVec<double>;
+template class NRVec<complex<double> >;
+template class NRVec<char>;
+template class NRVec<short>;
+template class NRVec<int>;
+template class NRVec<unsigned int>;
+template class NRVec<unsigned long>;
+