jiri/LA_library
1
0
Fork 0
Code Issues Pull Requests Projects Releases Wiki Activity

*** empty log message ***

Browse Source
This commit is contained in:
jiri 2006-09-18 21:46:45 +00:00
parent 7618bbd325
commit aaa74b8d40
3 changed files with 52 additions and 23 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

25
diis.h
View File

@ -9,7 +9,9 @@
#include "la_traits.h" #include "la_traits.h"
#include "auxstorage.h" #include "auxstorage.h"
// T is some solution vector in form of NRVec, NRMat, or NRSMat over double or complex<double> fields // Typically, T is some solution vector in form of NRVec, NRMat, or NRSMat over double or complex<double> fields
// actually it can be anything what has operator=(), dot() , axpy(), norm() and copyonwrite(), and LA_traits<T>::normtype and elementtype
// and get() and put() if external storage is requested
template<typename T> template<typename T>
class DIIS class DIIS
@ -23,7 +25,9 @@ class DIIS
AuxStorage<Te> *st; AuxStorage<Te> *st;
T *stor; T *stor;
public: public:
DIIS() {dim=0;}; //for array of diis
DIIS(const int n, const bool core=1); DIIS(const int n, const bool core=1);
void setup(const int n, const bool core=1);
~DIIS(); ~DIIS();
typename LA_traits<T>::normtype extrapolate(T &vec); //vec is input/output; returns square residual norm typename LA_traits<T>::normtype extrapolate(T &vec); //vec is input/output; returns square residual norm
}; };
@ -35,9 +39,21 @@ st=incore?NULL: new AuxStorage<Te>;
stor= incore? new T[dim] : NULL; stor= incore? new T[dim] : NULL;
bmat= (Te)0; for(int i=1; i<n; ++i) bmat(0,i) = (Te)-1; bmat= (Te)0; for(int i=1; i<n; ++i) bmat(0,i) = (Te)-1;
aktdim=cyclicshift=0; aktdim=cyclicshift=0;
} }
template<typename T>
void DIIS<T>::setup(const int n, const bool core)
{
dim=n;
incore=core;
bmat.resize(n);
st=incore?NULL: new AuxStorage<Te>;
stor= incore? new T[dim] : NULL;
bmat= (Te)0; for(int i=1; i<n; ++i) bmat(0,i) = (Te)-1;
aktdim=cyclicshift=0;
}
template<typename T> template<typename T>
DIIS<T>::~DIIS() DIIS<T>::~DIIS()
{ {
@ -49,6 +65,7 @@ if(stor) delete[] stor;
template<typename T> template<typename T>
typename LA_traits<T>::normtype DIIS<T>::extrapolate(T &vec) typename LA_traits<T>::normtype DIIS<T>::extrapolate(T &vec)
{ {
if(!dim) laerror("attempt to extrapolate from uninitialized DIIS");
//if dim exceeded, shift //if dim exceeded, shift
if(aktdim==dim) if(aktdim==dim)
{ {
@ -68,12 +85,12 @@ if(aktdim==1) return (typename LA_traits<T>::normtype)1000000000;
//calculate difference; //calculate difference;
vec.copyonwrite(); vec.copyonwrite();
if(incore) vec -= stor[(aktdim-2+cyclicshift)%dim]; if(incore) vec.axpy((Te)-1,stor[(aktdim-2+cyclicshift)%dim]);
else else
{ {
T tmp=vec; T tmp=vec;
st->get(tmp,(aktdim-2+cyclicshift)%dim); st->get(tmp,(aktdim-2+cyclicshift)%dim);
vec -= tmp; vec.axpy((Te)-1,tmp);
} }
//calculate overlaps of differences (if storage is cheap, they could rather be stored than recomputed) //calculate overlaps of differences (if storage is cheap, they could rather be stored than recomputed)

12
fourindex.h
View File

@ -600,8 +600,17 @@ istream& operator>>(istream &s, fourindex<I,T> &x)
//note - loops for the twoelectronrealmullikan integral to be unique and in canonical order //note - loops for the twoelectronrealmullikan integral to be unique and in canonical order
// i=1..n, j=1..i, k=1..i, l=1..(i==k?j:k) // i=1..n, j=1..i, k=1..i, l=1..(i==k?j:k)
//general template declaration
template<fourindexsymtype S, class T, class DUMMY> class fourindex_dense; template<fourindexsymtype S, class T, class DUMMY> class fourindex_dense;
//traits class
template<fourindexsymtype S, class T, class DUMMY>
struct LA_traits<fourindex_dense<S,T,DUMMY> > {
typedef T elementtype;
typedef typename LA_traits<T>::normtype normtype;
};
//make it as a derived class in order to be able to use it in a base class context - "supermatrix" operations //make it as a derived class in order to be able to use it in a base class context - "supermatrix" operations
template<class T, class I> template<class T, class I>
class fourindex_dense<twoelectronrealmullikan,T,I> : public NRSMat<T> { class fourindex_dense<twoelectronrealmullikan,T,I> : public NRSMat<T> {
@ -768,7 +777,4 @@ if(a<b) {minus++; unsigned int t=a; a=b; b=t;}
#endif /*_fourindex_included*/ #endif /*_fourindex_included*/

38
la_traits.h
View File

@ -25,7 +25,9 @@ extern "C" {
//forward declarations //forward declarations
template<typename C> class NRVec; template<typename C> class NRVec;
template<typename C> class NRMat; template<typename C> class NRMat;
template<typename C> class NRMat_from1;
template<typename C> class NRSMat; template<typename C> class NRSMat;
template<typename C> class NRSMat_from1;
template<typename C> class SparseMat; template<typename C> class SparseMat;
//for general sortable classes //for general sortable classes
@ -174,30 +176,34 @@ static void clear(C *dest, unsigned int n) {for(unsigned int i=0; i<n; ++i) dest
//non-scalar types defined in this library //non-scalar types defined in this library
generate_traits(NRMat) generate_traits(NRMat)
generate_traits(NRMat_from1)
generate_traits(NRVec) generate_traits(NRVec)
generate_traits(SparseMat) generate_traits(SparseMat)
#undef generate_traits #undef generate_traits
//smat //smat
template<typename C> #define generate_traits_smat(X) \
struct LA_traits_aux<NRSMat<C>, scalar_false> { template<typename C> \
typedef C elementtype; struct LA_traits_aux<X<C>, scalar_false> { \
typedef NRMat<C> producttype; typedef C elementtype; \
typedef typename LA_traits<C>::normtype normtype; typedef NRMat<C> producttype; \
static bool gencmp(const C *x, const C *y, int n) {for(int i=0; i<n; ++i) if(x[i]!=y[i]) return true; return false;} typedef typename LA_traits<C>::normtype normtype; \
static inline bool bigger(const C &x, const C &y) {return x>y;} static bool gencmp(const C *x, const C *y, int n) {for(int i=0; i<n; ++i) if(x[i]!=y[i]) return true; return false;} \
static inline bool smaller(const C &x, const C &y) {return x<y;} static inline bool bigger(const C &x, const C &y) {return x>y;} \
static inline normtype norm (const NRSMat<C> &x) {return x.norm();} static inline bool smaller(const C &x, const C &y) {return x<y;} \
static inline void axpy (NRSMat<C>&s, const NRSMat<C> &x, const C c) {s.axpy(c,x);} static inline normtype norm (const X<C> &x) {return x.norm();} \
static void put(int fd, const C &x, bool dimensions=1, bool transp=0) {x.put(fd,dimensions);} static inline void axpy (X<C>&s, const X<C> &x, const C c) {s.axpy(c,x);} \
static void get(int fd, C &x, bool dimensions=1, bool transp=0) {x.get(fd,dimensions);} static void put(int fd, const C &x, bool dimensions=1, bool transp=0) {x.put(fd,dimensions);} \
static void multiput(unsigned int n,int fd, const C *x, bool dimensions=1) {for(unsigned int i=0; i<n; ++i) x[i].put(fd,dimensions);} static void get(int fd, C &x, bool dimensions=1, bool transp=0) {x.get(fd,dimensions);} \
static void multiget(unsigned int n,int fd, C *x, bool dimensions=1) {for(unsigned int i=0; i<n; ++i) x[i].get(fd,dimensions);} static void multiput(unsigned int n,int fd, const C *x, bool dimensions=1) {for(unsigned int i=0; i<n; ++i) x[i].put(fd,dimensions);} \
static void copy(C *dest, C *src, unsigned int n) {for(unsigned int i=0; i<n; ++i) dest[i]=src[i];} static void multiget(unsigned int n,int fd, C *x, bool dimensions=1) {for(unsigned int i=0; i<n; ++i) x[i].get(fd,dimensions);} \
static void clear(C *dest, unsigned int n) {for(unsigned int i=0; i<n; ++i) dest[i].clear();} static void copy(C *dest, C *src, unsigned int n) {for(unsigned int i=0; i<n; ++i) dest[i]=src[i];} \
static void clear(C *dest, unsigned int n) {for(unsigned int i=0; i<n; ++i) dest[i].clear();} \
}; };
generate_traits_smat(NRSMat)
generate_traits_smat(NRSMat_from1)
//the final traits class //the final traits class
template<typename C> template<typename C>