00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018 #ifndef _davidson_h
00019 #define _davidson_h
00020 #include "vec.h"
00021 #include "smat.h"
00022 #include "mat.h"
00023 #include "sparsemat.h"
00024 #include "nonclass.h"
00025 #include "auxstorage.h"
00026
00027 namespace LA {
00028
00029
00030
00031
00032
00033
00034
00035
00036 template <typename T, typename Matrix>
00037 extern void davidson(const Matrix &bigmat, NRVec<T> &eivals, NRVec<T> *eivecs, const char *eivecsfile,
00038 int nroots=1, const bool verbose=0, const double eps=1e-6,
00039 const bool incore=1, int maxit=100, const int maxkrylov = 500,
00040 void (*initguess)(NRVec<T> &)=NULL);
00041
00042
00043
00044
00045
00046
00047
00048
00049
00050
00051 template <typename T, typename Matrix>
00052 void davidson(const Matrix &bigmat, NRVec<T> &eivals, NRVec<T> *eivecs, const char *eivecsfile,
00053 int nroots, const bool verbose, const double eps,
00054 const bool incore, int maxit, const int maxkrylov,
00055 void (*initguess)(NRVec<T> &))
00056 {
00057 bool flag=0;
00058 int n=bigmat.nrows();
00059 if ( n!= (int)bigmat.ncols()) laerror("non-square matrix in davidson");
00060 if(eivals.size()<nroots) laerror("too small eivals dimension in davidson");
00061
00062 NRVec<T> vec1(n),vec2(n);
00063 NRMat<T> smallH(maxkrylov,maxkrylov),smallS(maxkrylov,maxkrylov),smallV;
00064 NRVec<T> r(maxkrylov);
00065 NRVec<T> *v0,*v1;
00066 AuxStorage<T> *s0,*s1;
00067
00068 if(incore)
00069 {
00070 v0 = new NRVec<T>[maxkrylov];
00071 v1 = new NRVec<T>[maxkrylov];
00072 }
00073 else
00074 {
00075 s0 = new AuxStorage<T>;
00076 s1 = new AuxStorage<T>;
00077 }
00078
00079 int i,j;
00080
00081 if(maxkrylov<maxit) maxit=maxkrylov;
00082 if(nroots>=maxkrylov) nroots =maxkrylov-1;
00083 int nroot=0;
00084 int oldnroot;
00085 smallS=0;
00086 smallH=0;
00087
00088
00089
00090 if(initguess) (*initguess)(vec1);
00091 else
00092 {
00093 const T *diagonal = bigmat.diagonalof(vec2,false,true);
00094 T t=1e100; int i,j;
00095 vec1=0;
00096 for(i=0, j= -1; i<n; ++i) if(diagonal[i]<t) {t=diagonal[i]; j=i;}
00097 vec1[j]=1;
00098 }
00099
00100
00101 bigmat.gemv(0,vec2,'n',1,vec1);
00102 smallH(0,0) = vec1*vec2;
00103 smallS(0,0) = vec1*vec1;
00104 int krylovsize = 0;
00105 if(incore) v0[0]=vec1; else s0->put(vec1,0);
00106 if(incore) v1[0]=vec2; else s1->put(vec2,0);
00107
00108
00109
00110 int it;
00111 for(it=0; it<maxit; ++it)
00112 {
00113 if(it>0)
00114 {
00115
00116 if(incore) v0[krylovsize]=vec1; else s0->put(vec1,krylovsize);
00117 for(j=0; j<krylovsize; ++j)
00118 {
00119 if(!incore) s0->get(vec2,j);
00120 smallS(krylovsize,j) = smallS(j,krylovsize) = vec1*(incore?v0[j]:vec2);
00121 }
00122 smallS(krylovsize,krylovsize) = vec1*vec1;
00123 bigmat.gemv(0,vec2,'n',1,vec1);
00124 if(incore) v1[krylovsize]=vec2; else s1->put(vec2,krylovsize);
00125
00126
00127 smallH(krylovsize,krylovsize) = vec1*vec2;
00128 for(j=0; j<krylovsize; ++j)
00129 {
00130 if(!incore) s0->get(vec1,j);
00131 smallH(j,krylovsize) = (incore?v0[j]:vec1)*vec2;
00132 if(bigmat.issymmetric()) smallH(krylovsize,j) = smallH(j,krylovsize);
00133 }
00134 if(!bigmat.issymmetric())
00135 {
00136 if(!incore) s0->get(vec1,krylovsize);
00137 for(j=0; j<krylovsize; ++j)
00138 {
00139 if(!incore) s1->get(vec2,j);
00140 smallH(krylovsize,j) = incore? v1[j]*v0[krylovsize] :vec1*vec2;
00141 }
00142 }
00143 }
00144 smallV=smallH;
00145 NRMat<T> smallSwork=smallS;
00146 if(bigmat.issymmetric())
00147 diagonalize(smallV,r,1,1,krylovsize+1,&smallSwork,1);
00148 else
00149 {
00150 NRVec<T> ri(krylovsize+1),beta(krylovsize+1);
00151 NRMat<T> scratch;
00152 scratch=smallV;
00153 gdiagonalize(scratch, r, ri,NULL, &smallV, 1, krylovsize+1, 2, 0, &smallSwork, &beta);
00154 for(int i=0; i<=krylovsize; ++i) {r[i]/=beta[i]; ri[i]/=beta[i];}
00155 }
00156
00157 T eival_n=r[nroot];
00158 oldnroot=nroot;
00159 typename LA_traits<T>::normtype test=std::abs(smallV(krylovsize,nroot));
00160 if(test<eps) nroot++;
00161 if(it==0) nroot = 0;
00162 for(int iroot=0; iroot<=std::min(krylovsize,nroots-1); ++iroot)
00163 {
00164 test = std::abs(smallV(krylovsize,iroot));
00165 if(test>eps) nroot=std::min(nroot,iroot);
00166 if(verbose && iroot<=std::max(oldnroot,nroot))
00167 {
00168 std::cout <<"Davidson: iter="<<it <<" dim="<<krylovsize<<" root="<<iroot<<" energy="<<r[iroot]<<"\n";
00169 std::cout.flush();
00170 }
00171 }
00172
00173 if(verbose && oldnroot!=nroot) {std::cout <<"root no. "<<oldnroot<<" converged\n"; std::cout.flush();}
00174 if (nroot>=nroots) goto converged;
00175 if (it==maxit-1) break;
00176
00177 if (krylovsize==maxkrylov)
00178 {
00179 if(nroot!=0) {flag=1; goto finished;}
00180 smallH=0;
00181 smallS=0;
00182 vec1=0;
00183 for(i=0; i<=krylovsize; ++i)
00184 {
00185 if(!incore) s0->get(vec2,i);
00186 vec1.axpy(smallV(i,0),incore?v0[i]:vec2);
00187 }
00188 s0->put(vec1,0);
00189 vec1.normalize();
00190 krylovsize = 0;
00191 continue;
00192 }
00193
00194
00195 vec1=0;
00196 for(j=0; j<=krylovsize; ++j)
00197 {
00198 if(!incore) s0->get(vec2,j);
00199 vec1.axpy(-r[nroot]*smallV(j,nroot),incore?v0[j]:vec2);
00200 if(!incore) s1->get(vec2,j);
00201 vec1.axpy(smallV(j,nroot),incore?v1[j]:vec2);
00202 }
00203
00204 {
00205 const T *diagonal = bigmat.diagonalof(vec2,false,true);
00206 eival_n = r[nroot];
00207 for(i=0; i<n; ++i)
00208 {
00209 T denom = diagonal[i] - eival_n;
00210 denom = denom<0?-std::max(0.1,std::abs(denom)):std::max(0.1,std::abs(denom));
00211 vec1[i] /= denom;
00212 }
00213 }
00214
00215
00216 vec1.normalize();
00217 for(j=0; j<=krylovsize; ++j)
00218 {
00219 typename LA_traits<T>::normtype vnorm;
00220 if(!incore) s0->get(vec2,j);
00221 do {
00222 T ab = vec1*(incore?v0[j]:vec2) /smallS(j,j);
00223 vec1.axpy(-ab,incore?v0[j]:vec2);
00224 vnorm = vec1.norm();
00225 vec1 *= (1./vnorm);
00226 } while (vnorm<0.99);
00227 }
00228
00229
00230
00231
00232 ++krylovsize;
00233 }
00234 flag=1;
00235 goto finished;
00236
00237 converged:
00238 AuxStorage<typename LA_traits<T>::elementtype> *ev;
00239 if(eivecsfile) ev = new AuxStorage<typename LA_traits<T>::elementtype>(eivecsfile);
00240 if(verbose) {std::cout << "Davidson converged in "<<it<<" iterations.\n"; std::cout.flush();}
00241 for(nroot=0; nroot<nroots; ++nroot)
00242 {
00243 eivals[nroot]=r[nroot];
00244 if(eivecs)
00245 {
00246 vec1=0;
00247 for(j=0; j<=krylovsize; ++j )
00248 {
00249 if(!incore) s0->get(vec2,j);
00250 vec1.axpy(smallV(j,nroot),incore?v0[j]:vec2);
00251 }
00252 vec1.normalize();
00253 if(eivecs) eivecs[nroot]|=vec1;
00254 if(eivecsfile)
00255 {
00256 ev->put(vec1,nroot);
00257 }
00258 }
00259 }
00260
00261 if(eivecsfile) delete ev;
00262
00263 finished:
00264 if(incore) {delete[] v0; delete[] v1;}
00265 else {delete s0; delete s1;}
00266
00267 if(flag) laerror("no convergence in davidson");
00268 }
00269
00270 }
00271 #endif