00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028 #if !defined(SLR_HEADER)
00029 #define SLR_HEADER
00030
00031
00032 #include <unistd.h>
00033
00034 #include "realtype.h"
00035
00036 namespace LR {
00037 class VarVector;
00040 template<bool MultByS,bool SwapXY>
00041 class VarVectorProxyOp {
00042 public:
00043 const VarVector& vec;
00044 ergo_real scalar;
00045 explicit VarVectorProxyOp(const VarVector& a, ergo_real s=1.0) : vec(a), scalar(s) {}
00046 };
00047
00048
00052 class VarVector {
00053 ergo_real *v;
00054 public:
00055 char *fName;
00057 int nvar;
00058 unsigned onDisk:1;
00059 unsigned inMemory:1;
00061 VarVector(const VarVector& a) : fName(NULL), nvar(a.nvar),
00062 onDisk(0), inMemory(1) {
00063 if(nvar) {
00064 v = new ergo_real[nvar*2];
00065 memcpy(v, a.v, 2*nvar*sizeof(ergo_real));
00066 } else v = NULL;
00067 }
00068
00069 VarVector() : v(NULL), fName(NULL), nvar(0), onDisk(0), inMemory(1) {}
00070
00072 VarVector(int nbast, int nocc, const ergo_real *full_mat)
00073 : v(NULL), fName(NULL), nvar(0), onDisk(0), inMemory(1) {
00074 setFromFull(nbast, nocc, full_mat);
00075 }
00076
00077 ~VarVector() {
00078 if(v) delete v;
00079 if(fName) {
00080 unlink(fName);
00081 delete fName;
00082 }
00083 }
00084 ergo_real* x() const { return v; }
00085 ergo_real* y() const { return v + nvar; }
00086 void symorth(void);
00087 void setSize(int n) {
00088 if(nvar != n) {
00089 if(v) delete v; v = new ergo_real[2*n];
00090 nvar = n;
00091 onDisk = 0;
00092 }
00093 }
00094 int size() const { return nvar; }
00095 void print(const char *comment = NULL) {
00096 if(comment) puts(comment);
00097 for(int i=0; i<nvar*2; i++) printf("%15.10f\n", (double)v[i]);
00098 }
00099 void setFromFull(int nbast, int nocc, const ergo_real *full_mat);
00100 void setFull(int nbast, int nocc, ergo_real *full_mat) const;
00101 const ergo_real& operator[](int i) const { return v[i]; }
00102 ergo_real& operator[](int i) { return v[i]; }
00103 VarVector& operator=(ergo_real scalar) {
00104 for(int i=0; i<2*nvar; i++) v[i] = scalar;
00105 onDisk = 0;
00106 return *this;
00107 };
00108 VarVector& operator=(const VarVector& b) {
00109 if(this == &b) return *this;
00110 if(nvar != b.nvar) {
00111 nvar = b.nvar;
00112 if(v) delete v;
00113 v = new ergo_real[2*nvar];
00114 }
00115 memcpy(v, b.v, 2*nvar*sizeof(v[0]));
00116 onDisk = 0;
00117 return *this;
00118 }
00119
00120 VarVector&
00121 operator=(const VarVectorProxyOp<false, false >& proxy) {
00122 if (&proxy.vec == this)
00123 throw "VarVector self-assignment not-implemented";
00124 if(nvar != proxy.vec.nvar) {
00125 if(v) delete v;
00126 nvar = proxy.vec.nvar;
00127 v = new ergo_real[2*nvar];
00128 }
00129
00130 for(int i=0; i<2*nvar; i++)
00131 v[i] = proxy.scalar*proxy.vec[i];
00132 onDisk = 0;
00133 return *this;
00134 }
00135 VarVector&
00136 operator=(const VarVectorProxyOp<false, true >& proxy) {
00137 if (&proxy.vec == this)
00138 throw "VarVector self-assignment not-implemented";
00139 if(nvar != proxy.vec.nvar) {
00140 if(v) delete v;
00141 nvar = proxy.vec.nvar;
00142 v = new ergo_real[2*nvar];
00143 }
00144 for(int i=0; i<nvar; i++) {
00145 v[i] = proxy.scalar*proxy.vec[i+nvar];
00146 v[i+nvar] = proxy.scalar*proxy.vec[i];
00147 }
00148 onDisk = 0;
00149 return *this;
00150 }
00151
00153 void load(const char* tmpdir);
00155 void save(const char* tmpdir);
00157 void release(const char* tmpdir);
00158
00159 friend ergo_real operator*(const VarVector& a, const VarVector& b);
00160 friend ergo_real operator*(const VarVector &a,
00161 const VarVectorProxyOp<false,false>& b);
00162 friend ergo_real operator*(const VarVector &a,
00163 const VarVectorProxyOp<true,false>& b);
00164 };
00165
00169 class E2Evaluator {
00170 public:
00171 virtual bool transform(const ergo_real *dmat, ergo_real *fmat) = 0;
00172 virtual ~E2Evaluator() {}
00173 };
00174
00175
00177 class VarVectorCollection {
00178 VarVector *vecs;
00179 unsigned *ages;
00180 unsigned currentAge;
00181 int nVecs;
00182 int nAllocated;
00183 bool diskMode;
00184 public:
00185 explicit VarVectorCollection(int nSize=0) : vecs(NULL), ages(NULL), currentAge(0),
00186 nVecs(0), nAllocated(0), diskMode(false) {
00187 if(nSize) setSize(nSize);
00188 }
00189 ~VarVectorCollection();
00190 void setSize(int sz);
00191 VarVector& operator[](int i);
00192 int size() const { return nVecs; }
00193 bool getDiskMode() const { return diskMode; }
00194 void setDiskMode(bool x) { diskMode = x; }
00196 void release();
00198 void releaseAll();
00199 static const char *tmpdir;
00200 };
00201
00202
00204 class OneElOperator {
00205 public:
00206 virtual void getOper(ergo_real *result) = 0;
00207 virtual ~OneElOperator() {}
00208 };
00209
00211 class SmallMatrix {
00212 ergo_real *mat;
00213 int nsize;
00214 protected:
00215 struct RowProxy {
00216 ergo_real *toprow;
00217 explicit RowProxy(ergo_real *r) : toprow(r) {}
00218 ergo_real& operator[](int j) const {
00219
00220 return toprow[j]; }
00221 };
00222 public:
00223 explicit SmallMatrix(int sz) : mat(new ergo_real[sz*sz]), nsize(sz) {}
00224 ~SmallMatrix() { delete mat; }
00225 const RowProxy operator[](int i) {
00226
00227 return RowProxy(mat + i*nsize); }
00228 void expand(int newSize);
00229 };
00230
00231
00232
00235 class LRSolver {
00236 public:
00237
00238 LRSolver(int nbast, int nocc,
00239 const ergo_real *fock_matrix,
00240 const ergo_real *s);
00241 virtual ~LRSolver() {
00242 if(cmo) delete cmo;
00243 if(fdiag) delete fdiag;
00244 delete xSub;
00245 }
00246
00253 virtual bool getResidual(VarVectorCollection& residualv) = 0;
00254
00258 virtual int getInitialGuess(VarVectorCollection& vecs) = 0;
00259
00262 virtual ergo_real getPreconditionerShift(int i) const = 0;
00263
00265 virtual void increaseSubspaceLimit(int newSize);
00266
00268 bool solve(E2Evaluator& e, bool diskMode = false);
00269 void computeExactE2Diag(E2Evaluator& e2);
00270 ergo_real convThreshold;
00271 int maxSubspaceSize;
00272 protected:
00273 static const int MVEC = 200;
00274 VarVector e2diag;
00275 int subspaceSize;
00276 SmallMatrix eSub;
00277 SmallMatrix sSub;
00278 ergo_real *xSub;
00282 void getAvMinusFreqSv(ergo_real f, ergo_real *weights,
00283 VarVector& r);
00284
00289 void projectOnSubspace(const VarVector& full, ergo_real *w);
00291 void buildVector(const ergo_real *w, VarVector& full) ;
00292
00294 void operToVec(OneElOperator& oper, VarVector& res) const;
00295
00299 ergo_real setE2diag(int nbast, int nocc,
00300 const ergo_real *fock_matrix,
00301 const ergo_real *s);
00302
00303 int nbast;
00304 int nocc;
00305 VarVectorCollection vects;
00306 virtual void addToSpace(VarVectorCollection& vecs, E2Evaluator& e2);
00307 void mo2ao(int nbast, const ergo_real *mo, ergo_real *ao) const;
00308 void ao2mo(int nbast, const ergo_real *ao, ergo_real *mo) const;
00309 private:
00313 VarVectorCollection Avects;
00314 ergo_real *fdiag;
00316 ergo_real *cmo;
00318 void load_F_MO(ergo_real *fmat) const;
00319 bool lintrans(E2Evaluator& e2, const VarVector& v, VarVector& Av) const;
00320 };
00321
00322
00325 class SetOfEqSolver : public LRSolver {
00326 ergo_real frequency;
00327 VarVector rhs;
00328 public:
00331 SetOfEqSolver(int nbast, int nocc,
00332 const ergo_real *fock_matrix,
00333 const ergo_real *s,
00334 ergo_real freq)
00335 : LRSolver(nbast, nocc, fock_matrix, s), frequency(freq),
00336 rhsSub(new ergo_real[MVEC]) {};
00337 void setRHS(OneElOperator& op);
00338 virtual ~SetOfEqSolver() { delete rhsSub; }
00339 virtual ergo_real getPreconditionerShift(int) const { return frequency; }
00340 virtual int getInitialGuess(VarVectorCollection& vecs);
00341 virtual bool getResidual(VarVectorCollection& residualv);
00343 virtual void increaseSubspaceLimit(int newSize);
00344 ergo_real getPolarisability(OneElOperator& oper) ;
00345 protected:
00346 ergo_real *rhsSub;
00347 virtual void addToSpace(VarVectorCollection& vecs, E2Evaluator& e2);
00348 ergo_real multiplyXtimesVec(const VarVector& rhs) ;
00349 ergo_real xTimesRHS;
00350 };
00351
00352
00353
00355 class EigenSolver : public LRSolver {
00356 ergo_real *ritzVals;
00357 ergo_real *transMoms2;
00358 int nStates;
00359 int nConverged;
00360 ergo_real *last_ev;
00361 public:
00362 EigenSolver(int nbast, int nocc,
00363 const ergo_real *fock_matrix,
00364 const ergo_real *overlap, int n)
00365 : LRSolver(nbast, nocc, NULL, NULL),
00366 ritzVals(new ergo_real[MVEC]), transMoms2(new ergo_real[MVEC]),
00367 nStates(n), nConverged(0), last_ev(NULL) {
00368 ritzVals[0] = setE2diag(nbast, nocc, fock_matrix, overlap);
00369 for(int i=1; i<n; i++) ritzVals[i] = ritzVals[0];
00370 }
00371 virtual ~EigenSolver() {
00372 if(last_ev) delete last_ev;
00373 delete ritzVals;
00374 delete transMoms2;
00375 }
00376 virtual ergo_real getPreconditionerShift(int i) const {
00377 return ritzVals[nConverged+i];
00378 }
00379 virtual int getInitialGuess(VarVectorCollection& vecs);
00380 virtual bool getResidual(VarVectorCollection& residualv);
00382 virtual void increaseSubspaceLimit(int newSize);
00383
00384 ergo_real getFreq(int i) const { return ritzVals[i]; }
00385 void computeMoments(OneElOperator& dipx,
00386 OneElOperator& dipy,
00387 OneElOperator& dipz);
00388 ergo_real getTransitionMoment2(int i) const { return transMoms2[i]; }
00389 };
00390
00391 }
00392 #endif