From b0a83b6d6ebd34ed9bf9d79ef569b39d39f8ad18 Mon Sep 17 00:00:00 2001 From: jiri Date: Sun, 10 Sep 2006 20:06:44 +0000 Subject: [PATCH] *** empty log message *** --- fourindex.h | 4 +-- mat.cc | 22 ++++++++-------- mat.h | 1 + matexp.h | 14 +++++----- nonclass.cc | 2 ++ nonclass.h | 2 ++ smat.cc | 22 ++++++++-------- sparsemat.cc | 12 ++++----- t.cc | 72 ++++++++++++++++++++++++++++++---------------------- vec.cc | 53 ++++++++++++++++++++++++++++++++------ 10 files changed, 130 insertions(+), 74 deletions(-) diff --git a/fourindex.h b/fourindex.h index ec6dbad..8fda646 100644 --- a/fourindex.h +++ b/fourindex.h @@ -622,7 +622,7 @@ public: }; template -fourindex_dense::fourindex_dense(const fourindex &rhs) : NRSMat((T)0,rhs.size()*(rhs.size()+1)/2) +fourindex_dense::fourindex_dense(const fourindex &rhs) : NRSMat((T)0,rhs.size()*(rhs.size()+1)/2) { if(rhs.getsymmetry() != twoelectronrealmullikan ) laerror("fourindex_dense symmetry mismatch"); typename fourindex::iterator p; @@ -635,7 +635,7 @@ for(p=rhs.begin(); p!= rhs.end(); ++p) (*this)(p->index.indiv.i,p->index.indiv.j } template -fourindex_dense::fourindex_dense(const fourindex_ext &rhs) : NRSMat((T)0,rhs.size()*(rhs.size()+1)/2) +fourindex_dense::fourindex_dense(const fourindex_ext &rhs) : NRSMat((T)0,rhs.size()*(rhs.size()+1)/2) { if(rhs.getsymmetry() != twoelectronrealmullikan ) laerror("fourindex_dense symmetry mismatch"); typename fourindex_ext::iterator p; diff --git a/mat.cc b/mat.cc index 924e442..eacf3dc 100644 --- a/mat.cc +++ b/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; -template class NRMat >; -template class NRMat; -template class NRMat; -template class NRMat; -template class NRMat; -template class NRMat; -template class NRMat; - /* * Templates first, specializations for BLAS next @@ -1113,3 +1102,14 @@ istream& operator>>(istream &s, NRMat &x) +////////////////////////////////////////////////////////////////////////////// +//// forced instantization in the corresponding object file +template class NRMat; +template class NRMat >; +template class NRMat; +template class NRMat; +template class NRMat; +template class NRMat; +template class NRMat; +template class NRMat; + diff --git a/mat.h b/mat.h index 47136d6..8136d73 100644 --- a/mat.h +++ b/mat.h @@ -560,6 +560,7 @@ public: NRMat_from1(): NRMat() {}; explicit NRMat_from1(const int n): NRMat(n) {}; NRMat_from1(const NRMat &rhs): NRMat(rhs) {}; //be able to convert the parent class transparently to this + NRMat_from1(const int n, const int m): NRMat(n,m) {}; NRMat_from1(const T &a, const int n, const int m): NRMat(a,n,m) {}; NRMat_from1(const T *a, const int n, const int m): NRMat(a,n,m) {}; diff --git a/matexp.h b/matexp.h index 269b0a2..66a7f97 100644 --- a/matexp.h +++ b/matexp.h @@ -109,7 +109,7 @@ return z; //general BCH expansion (can be written more efficiently in a specialization for matrices) template -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 "<::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::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 -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::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 diff --git a/nonclass.cc b/nonclass.cc index 44b4189..47606bc 100644 --- a/nonclass.cc +++ b/nonclass.cc @@ -617,7 +617,9 @@ NRMat matrixfunction(NRMat a, complex int n = a.nrows(); NRMat< complex > u(n, n), v(n, n); NRVec< complex > w(n); +/* NRMat > a0=complexify(a); +*/ gdiagonalize(a, w, &u, &v);//a gets destroyed, eigenvectors are rows NRVec< complex > z = diagofproduct(u, v, 1, 1); /* diff --git a/nonclass.h b/nonclass.h index e216609..00c5cda 100644 --- a/nonclass.h +++ b/nonclass.h @@ -70,6 +70,7 @@ extern const NRVec diagofproduct(const NRMat &a, const NRMat &b,\ bool trb=0, bool conjb=0); \ extern T trace2(const NRMat &a, const NRMat &b, bool trb=0); \ extern T trace2(const NRSMat &a, const NRSMat &b, const bool diagscaled=0);\ +extern T trace2(const NRSMat &a, const NRMat &b, const bool diagscaled=0);\ extern void linear_solve(NRMat &a, NRMat *b, double *det=0,int n=0); \ extern void linear_solve(NRSMat &a, NRMat *b, double *det=0, int n=0); \ extern void linear_solve(NRMat &a, NRVec &b, double *det=0, int n=0); \ @@ -108,6 +109,7 @@ inline NRMat realsqrt(const NRMat &a) { return realmatrixfuncti inline NRMat realsqrtinv(const NRMat &a) { return realmatrixfunction(a,&sqrtinv); } inline NRMat log(const NRSMat &a) { return matrixfunction(a,&log); } extern NRMat log(const NRMat &a); +extern NRMat exp0(const NRMat &a); extern const NRMat realpart(const NRMat< complex >&); diff --git a/smat.cc b/smat.cc index 3269926..1d27661 100644 --- a/smat.cc +++ b/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; -template class NRSMat< complex >; -template class NRSMat; -template class NRSMat; -template class NRSMat; -template class NRSMat; -template class NRSMat; - /* @@ -411,3 +401,15 @@ INSTANTIZE(char) INSTANTIZE(unsigned int) INSTANTIZE(unsigned long) + + +////////////////////////////////////////////////////////////////////////////// +////// forced instantization in the corresponding object file +template class NRSMat; +template class NRSMat< complex >; +template class NRSMat; +template class NRSMat; +template class NRSMat; +template class NRSMat; +template class NRSMat; + diff --git a/sparsemat.cc b/sparsemat.cc index 8a12b79..e9ab6b6 100644 --- a/sparsemat.cc +++ b/sparsemat.cc @@ -17,11 +17,6 @@ static inline void SWAP(T &a, T &b) -////////////////////////////////////////////////////////////////////////////// -//// forced instantization in the corresponding object file -template class SparseMat; -template class SparseMat >; - export template @@ -1327,13 +1322,13 @@ template unsigned int SparseMat::length() const; \ template void SparseMat::deletelist(); \ template void SparseMat::simplify(); \ template void SparseMat::copylist(const matel *l); \ -template void SparseMat::add(const SPMatindex n, const SPMatindex m, const T elem); \ template SparseMat & SparseMat::operator=(const SparseMat &rhs); \ template SparseMat & SparseMat::operator+=(const SparseMat &rhs); \ template SparseMat & SparseMat::operator-=(const SparseMat &rhs); \ template SparseMat::SparseMat(const NRMat &rhs); \ template SparseMat::SparseMat(const NRSMat &rhs); \ template void SparseMat::transposeme(); \ +template const T* SparseMat::diagonalof(NRVec &r, const bool divide, bool cache) const; \ template SparseMat & SparseMat::operator*=(const T a); \ template void SparseMat::setunsymmetric(); \ template SparseMat & SparseMat::operator=(const T a); \ @@ -1362,3 +1357,8 @@ INSTANTIZE(double) INSTANTIZE(complex) //some functions are not OK for hermitean matrices, needs a revision!!! +////////////////////////////////////////////////////////////////////////////// +//// forced instantization in the corresponding object file +template class SparseMat; +template class SparseMat >; + diff --git a/t.cc b/t.cc index 9cdce7e..1a18ac0 100644 --- a/t.cc +++ b/t.cc @@ -52,7 +52,8 @@ NRVec x(1.,10); NRVec y(2.,10); NRVec 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 a(n,n); @@ -655,39 +655,49 @@ for(int i=0;i b=exp(a); -cout < a,b; -cin >>b; -int n=b.nrows(); -cout <<"difference from identity = "< b=exp0(a); +cout <<"b=exp(a) matrix\n"< x(0.,n,n),x0; - double r; -int i=0; -do - { - x0=x; - NRMat y=exp(x*-.5); - x+= y*b*y; - x-= 1.; - x=(x-x.transpose())*.5; - cout <<"matrix x\n"<1e-10); -cout <<"result\n"< c=log(b); //matrixfunction(a,&mycident,1); -cout < d=exp(c); -cout <<"exp(log(x))\n"< e=exp(c); +cout <<"e=exp(c)\n"< b=exp0(a); +cout <<"exp0 took "< c=exp(a,true); +cout <<" horner exp took "< NRVec & NRVec::normalize() {laerror("normalize() impossible for integer types"); return *this;} +template<> +NRVec & NRVec::normalize() {laerror("normalize() impossible for integer types"); return *this;} template<> @@ -314,6 +309,15 @@ void NRVec::gemv(const unsigned long beta, laerror("not yet implemented"); } +template<> +void NRVec::gemv(const unsigned int beta, + const NRSMat &A, const char trans, + const unsigned int alpha, const NRVec &x) +{ +laerror("not yet implemented"); +} + + template<> void NRVec::gemv(const unsigned char beta, const NRSMat &A, const char trans, @@ -355,6 +359,15 @@ void NRVec::gemv(const unsigned long beta, laerror("not yet implemented"); } +template<> +void NRVec::gemv(const unsigned int beta, + const NRMat &A, const char trans, + const unsigned int alpha, const NRVec &x) +{ +laerror("not yet implemented"); +} + + template<> void NRVec::gemv(const unsigned char beta, const NRMat &A, const char trans, @@ -407,6 +420,16 @@ void NRVec::gemv(const unsigned long beta, laerror("not yet implemented"); } +template<> +void NRVec::gemv(const unsigned int beta, + const SparseMat &A, const char trans, + const unsigned int alpha, const NRVec &x, bool s) +{ +laerror("not yet implemented"); +} + + + template<> void NRVec::gemv(const unsigned char beta, const SparseMat &A, const char trans, @@ -501,3 +524,17 @@ if(to == -1) to=nn-1; if(direction) return memqsort<1,NRVec,int,int>(*this,perm,from,to); else return memqsort<0,NRVec,int,int>(*this,perm,from,to); } + + +////////////////////////////////////////////////////////////////////////////// +//// forced instantization in the corespoding object file + + +template class NRVec; +template class NRVec >; +template class NRVec; +template class NRVec; +template class NRVec; +template class NRVec; +template class NRVec; +