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
00029
00038 #if !defined(SLR_HEADER)
00039 #define SLR_HEADER
00040
00041
00042 #include <unistd.h>
00043
00044 #include "realtype.h"
00045
00046 namespace LR {
00047 class VarVector;
00050 template<bool MultByS,bool SwapXY>
00051 class VarVectorProxyOp {
00052 public:
00053 const VarVector& vec;
00054 ergo_real scalar;
00055 explicit VarVectorProxyOp(const VarVector& a, ergo_real s=1.0) : vec(a), scalar(s) {}
00056 };
00057
00058
00062 class VarVector {
00063 ergo_real *v;
00064 public:
00065 char *fName;
00067 int nvar;
00068 unsigned onDisk:1;
00069 unsigned inMemory:1;
00071 VarVector(const VarVector& a) : fName(NULL), nvar(a.nvar),
00072 onDisk(0), inMemory(1) {
00073 if(nvar) {
00074 v = new ergo_real[nvar*2];
00075 memcpy(v, a.v, 2*nvar*sizeof(ergo_real));
00076 } else v = NULL;
00077 }
00078
00079 VarVector() : v(NULL), fName(NULL), nvar(0), onDisk(0), inMemory(1) {}
00080
00082 VarVector(int nbast, int nocc, const ergo_real *full_mat)
00083 : v(NULL), fName(NULL), nvar(0), onDisk(0), inMemory(1) {
00084 setFromFull(nbast, nocc, full_mat);
00085 }
00086
00087 ~VarVector() {
00088 if(v) delete v;
00089 if(fName) {
00090 unlink(fName);
00091 delete fName;
00092 }
00093 }
00094 ergo_real* x() const { return v; }
00095 ergo_real* y() const { return v + nvar; }
00096 void symorth(void);
00097 void setSize(int n) {
00098 if(nvar != n) {
00099 if(v)
00100 delete v;
00101 v = new ergo_real[2*n];
00102 nvar = n;
00103 onDisk = 0;
00104 }
00105 }
00106 int size() const { return nvar; }
00107 void print(const char *comment = NULL) {
00108 if(comment) puts(comment);
00109 for(int i=0; i<nvar*2; i++) printf("%15.10f\n", (double)v[i]);
00110 }
00111 void setFromFull(int nbast, int nocc, const ergo_real *full_mat);
00112 void setFull(int nbast, int nocc, ergo_real *full_mat) const;
00113 const ergo_real& operator[](int i) const { return v[i]; }
00114 ergo_real& operator[](int i) { return v[i]; }
00115 VarVector& operator=(ergo_real scalar) {
00116 for(int i=0; i<2*nvar; i++) v[i] = scalar;
00117 onDisk = 0;
00118 return *this;
00119 };
00120 VarVector& operator=(const VarVector& b) {
00121 if(this == &b) return *this;
00122 if(nvar != b.nvar) {
00123 nvar = b.nvar;
00124 if(v) delete v;
00125 v = new ergo_real[2*nvar];
00126 }
00127 memcpy(v, b.v, 2*nvar*sizeof(v[0]));
00128 onDisk = 0;
00129 return *this;
00130 }
00131
00132 VarVector&
00133 operator=(const VarVectorProxyOp<false, false >& proxy) {
00134 if (&proxy.vec == this)
00135 throw "VarVector self-assignment not-implemented";
00136 if(nvar != proxy.vec.nvar) {
00137 if(v) delete v;
00138 nvar = proxy.vec.nvar;
00139 v = new ergo_real[2*nvar];
00140 }
00141
00142 for(int i=0; i<2*nvar; i++)
00143 v[i] = proxy.scalar*proxy.vec[i];
00144 onDisk = 0;
00145 return *this;
00146 }
00147 VarVector&
00148 operator=(const VarVectorProxyOp<false, true >& proxy) {
00149 if (&proxy.vec == this)
00150 throw "VarVector self-assignment not-implemented";
00151 if(nvar != proxy.vec.nvar) {
00152 if(v) delete v;
00153 nvar = proxy.vec.nvar;
00154 v = new ergo_real[2*nvar];
00155 }
00156 for(int i=0; i<nvar; i++) {
00157 v[i] = proxy.scalar*proxy.vec[i+nvar];
00158 v[i+nvar] = proxy.scalar*proxy.vec[i];
00159 }
00160 onDisk = 0;
00161 return *this;
00162 }
00163
00165 void load(const char* tmpdir);
00167 void save(const char* tmpdir);
00169 void release(const char* tmpdir);
00170
00171 friend ergo_real operator*(const VarVector& a, const VarVector& b);
00172 friend ergo_real operator*(const VarVector &a,
00173 const VarVectorProxyOp<false,false>& b);
00174 friend ergo_real operator*(const VarVector &a,
00175 const VarVectorProxyOp<true,false>& b);
00176 };
00177
00181 class E2Evaluator {
00182 public:
00183 virtual bool transform(const ergo_real *dmat, ergo_real *fmat) = 0;
00184 virtual ~E2Evaluator() {}
00185 };
00186
00187
00189 class VarVectorCollection {
00190 VarVector *vecs;
00191 unsigned *ages;
00192 unsigned currentAge;
00193 int nVecs;
00194 int nAllocated;
00195 bool diskMode;
00196 public:
00197 explicit VarVectorCollection(int nSize=0) : vecs(NULL), ages(NULL), currentAge(0),
00198 nVecs(0), nAllocated(0), diskMode(false) {
00199 if(nSize) setSize(nSize);
00200 }
00201 ~VarVectorCollection();
00202 void setSize(int sz);
00203 VarVector& operator[](int i);
00204 int size() const { return nVecs; }
00205 bool getDiskMode() const { return diskMode; }
00206 void setDiskMode(bool x) { diskMode = x; }
00208 void release();
00210 void releaseAll();
00211 static const char *tmpdir;
00212 };
00213
00214
00216 class OneElOperator {
00217 public:
00218 virtual void getOper(ergo_real *result) = 0;
00219 virtual ~OneElOperator() {}
00220 };
00221
00223 class SmallMatrix {
00224 ergo_real *mat;
00225 int nsize;
00226 protected:
00227 struct RowProxy {
00228 ergo_real *toprow;
00229 explicit RowProxy(ergo_real *r) : toprow(r) {}
00230 ergo_real& operator[](int j) const {
00231
00232 return toprow[j]; }
00233 };
00234 public:
00235 explicit SmallMatrix(int sz) : mat(new ergo_real[sz*sz]), nsize(sz) {}
00236 ~SmallMatrix() { delete [] mat; }
00237 const RowProxy operator[](int i) {
00238
00239 return RowProxy(mat + i*nsize); }
00240 void expand(int newSize);
00241 };
00242
00243
00244
00247 class LRSolver {
00248 public:
00249
00250 LRSolver(int nbast, int nocc,
00251 const ergo_real *fock_matrix,
00252 const ergo_real *s);
00253 virtual ~LRSolver() {
00254 if(cmo) delete cmo;
00255 if(fdiag) delete fdiag;
00256 delete [] xSub;
00257 }
00258
00265 virtual bool getResidual(VarVectorCollection& residualv) = 0;
00266
00270 virtual int getInitialGuess(VarVectorCollection& vecs) = 0;
00271
00274 virtual ergo_real getPreconditionerShift(int i) const = 0;
00275
00277 virtual void increaseSubspaceLimit(int newSize);
00278
00280 bool solve(E2Evaluator& e, bool diskMode = false);
00281 void computeExactE2Diag(E2Evaluator& e2);
00282 ergo_real convThreshold;
00283 int maxSubspaceSize;
00284 protected:
00285 static const int MVEC = 200;
00286 VarVector e2diag;
00287 int subspaceSize;
00288 SmallMatrix eSub;
00289 SmallMatrix sSub;
00290 ergo_real *xSub;
00294 void getAvMinusFreqSv(ergo_real f, ergo_real *weights,
00295 VarVector& r);
00296
00301 void projectOnSubspace(const VarVector& full, ergo_real *w);
00303 void buildVector(const ergo_real *w, VarVector& full) ;
00304
00306 void operToVec(OneElOperator& oper, VarVector& res) const;
00307
00311 ergo_real setE2diag(int nbast, int nocc,
00312 const ergo_real *fock_matrix,
00313 const ergo_real *s);
00314
00315 int nbast;
00316 int nocc;
00317 VarVectorCollection vects;
00318 virtual void addToSpace(VarVectorCollection& vecs, E2Evaluator& e2);
00319 void mo2ao(int nbast, const ergo_real *mo, ergo_real *ao) const;
00320 void ao2mo(int nbast, const ergo_real *ao, ergo_real *mo) const;
00321 private:
00325 VarVectorCollection Avects;
00326 ergo_real *fdiag;
00328 ergo_real *cmo;
00330 void load_F_MO(ergo_real *fmat) const;
00331 bool lintrans(E2Evaluator& e2, const VarVector& v, VarVector& Av) const;
00332 };
00333
00334
00337 class SetOfEqSolver : public LRSolver {
00338 ergo_real frequency;
00339 VarVector rhs;
00340 public:
00343 SetOfEqSolver(int nbast, int nocc,
00344 const ergo_real *fock_matrix,
00345 const ergo_real *s,
00346 ergo_real freq)
00347 : LRSolver(nbast, nocc, fock_matrix, s), frequency(freq),
00348 rhsSub(new ergo_real[MVEC]) {};
00349 void setRHS(OneElOperator& op);
00350 virtual ~SetOfEqSolver() { delete [] rhsSub; }
00351 virtual ergo_real getPreconditionerShift(int) const { return frequency; }
00352 virtual int getInitialGuess(VarVectorCollection& vecs);
00353 virtual bool getResidual(VarVectorCollection& residualv);
00355 virtual void increaseSubspaceLimit(int newSize);
00356 ergo_real getPolarisability(OneElOperator& oper) ;
00357 protected:
00358 ergo_real *rhsSub;
00359 virtual void addToSpace(VarVectorCollection& vecs, E2Evaluator& e2);
00360 ergo_real multiplyXtimesVec(const VarVector& rhs) ;
00361 ergo_real xTimesRHS;
00362 };
00363
00364
00365
00367 class EigenSolver : public LRSolver {
00368 ergo_real *ritzVals;
00369 ergo_real *transMoms2;
00370 int nStates;
00371 int nConverged;
00372 ergo_real *last_ev;
00373 public:
00374 EigenSolver(int nbast, int nocc,
00375 const ergo_real *fock_matrix,
00376 const ergo_real *overlap, int n)
00377 : LRSolver(nbast, nocc, NULL, NULL),
00378 ritzVals(new ergo_real[MVEC]), transMoms2(new ergo_real[MVEC]),
00379 nStates(n), nConverged(0), last_ev(NULL) {
00380 ritzVals[0] = setE2diag(nbast, nocc, fock_matrix, overlap);
00381 for(int i=1; i<n; i++) ritzVals[i] = ritzVals[0];
00382 }
00383 virtual ~EigenSolver() {
00384 if(last_ev) delete last_ev;
00385 delete [] ritzVals;
00386 delete [] transMoms2;
00387 }
00388 virtual ergo_real getPreconditionerShift(int i) const {
00389 return ritzVals[nConverged+i];
00390 }
00391 virtual int getInitialGuess(VarVectorCollection& vecs);
00392 virtual bool getResidual(VarVectorCollection& residualv);
00394 virtual void increaseSubspaceLimit(int newSize);
00395
00396 ergo_real getFreq(int i) const { return ritzVals[i]; }
00397 void computeMoments(OneElOperator& dipx,
00398 OneElOperator& dipy,
00399 OneElOperator& dipz);
00400 ergo_real getTransitionMoment2(int i) const { return transMoms2[i]; }
00401 };
00402
00403 }
00404 #endif