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 2009-09-04 08:09:32 +00:00
parent fec09b23da
commit e82adf91f8
5 changed files with 26 additions and 17 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

2
bitvector.h
View File

@ -53,7 +53,7 @@ public:
void reset(const unsigned int i) {v[i/blockbits] &= ~(1<<(i%blockbits));}; void reset(const unsigned int i) {v[i/blockbits] &= ~(1<<(i%blockbits));};
const bool get(const unsigned int i) {return (v[i/blockbits] >>(i%blockbits))&1;}; const bool get(const unsigned int i) {return (v[i/blockbits] >>(i%blockbits))&1;};
const bool assign(const unsigned int i, const bool r) {if(r) set(i); else reset(i); return r;}; const bool assign(const unsigned int i, const bool r) {if(r) set(i); else reset(i); return r;};
void clear() {memset(v,0,nn*sizeof(bitvector_block));}; void clear() {copyonwrite(); memset(v,0,nn*sizeof(bitvector_block));};
void fill() {memset(v,0xff,nn*sizeof(bitvector_block));}; void fill() {memset(v,0xff,nn*sizeof(bitvector_block));};
bool operator!=(const bitvector &rhs) const; bool operator!=(const bitvector &rhs) const;
bool operator==(const bitvector &rhs) const {return !(*this != rhs);}; bool operator==(const bitvector &rhs) const {return !(*this != rhs);};

15
mat.h
View File

@ -42,6 +42,7 @@ public:
inline NRMat(const T &a, const int n, const int m); inline NRMat(const T &a, const int n, const int m);
NRMat(const T *a, const int n, const int m); NRMat(const T *a, const int n, const int m);
inline NRMat(const NRMat &rhs); inline NRMat(const NRMat &rhs);
NRMat(const typename LA_traits_complex<T>::NRMat_Noncomplex_type &rhs, bool imagpart=false); //construct complex from real
explicit NRMat(const NRSMat<T> &rhs); explicit NRMat(const NRSMat<T> &rhs);
#ifdef MATPTR #ifdef MATPTR
explicit NRMat(const NRVec<T> &rhs, const int n, const int m, const int offset=0) :NRMat(&rhs[0][0] + offset ,n,m) {if (offset < 0 || n*m + offset > rhs.nn) laerror("matrix dimensions and offset incompatible with vector length");}; explicit NRMat(const NRVec<T> &rhs, const int n, const int m, const int offset=0) :NRMat(&rhs[0][0] + offset ,n,m) {if (offset < 0 || n*m + offset > rhs.nn) laerror("matrix dimensions and offset incompatible with vector length");};
@ -57,8 +58,7 @@ public:
const bool operator==(const NRMat &rhs) const {return !(*this != rhs);}; const bool operator==(const NRMat &rhs) const {return !(*this != rhs);};
inline int getcount() const {return count?*count:0;} inline int getcount() const {return count?*count:0;}
NRMat & operator=(const NRMat &rhs); //assignment NRMat & operator=(const NRMat &rhs); //assignment
void clear() {if(nn&&mm) LA_traits<T>::clear((*this)[0],nn*mm);}; //zero out void randomize(const typename LA_traits<T>::normtype &x); //fill with random numbers
void randomize(const T &x); //fill with random numbers
NRMat & operator=(const T &a); //assign a to diagonal NRMat & operator=(const T &a); //assign a to diagonal
NRMat & operator|=(const NRMat &rhs); //assignment to a new copy NRMat & operator|=(const NRMat &rhs); //assignment to a new copy
NRMat & operator+=(const T &a); //add diagonal NRMat & operator+=(const T &a); //add diagonal
@ -88,6 +88,7 @@ public:
const NRMat operator&(const NRMat &rhs) const; // direct sum const NRMat operator&(const NRMat &rhs) const; // direct sum
const NRMat operator|(const NRMat<T> &rhs) const; // direct product const NRMat operator|(const NRMat<T> &rhs) const; // direct product
const NRVec<T> operator*(const NRVec<T> &rhs) const {NRVec<T> result(nn); result.gemv((T)0,*this,'n',(T)1,rhs); return result;}; // Mat * Vec const NRVec<T> operator*(const NRVec<T> &rhs) const {NRVec<T> result(nn); result.gemv((T)0,*this,'n',(T)1,rhs); return result;}; // Mat * Vec
const NRVec<complex<T> > operator*(const NRVec<complex<T> > &rhs) const {NRVec<complex<T> > result(nn); result.gemv((T)0,*this,'n',(T)1,rhs); return result;}; // Mat * Vec
const NRVec<T> rsum() const; //sum of rows const NRVec<T> rsum() const; //sum of rows
const NRVec<T> csum() const; //sum of columns const NRVec<T> csum() const; //sum of columns
const NRVec<T> row(const int i, int l= -1) const; //row of, efficient const NRVec<T> row(const int i, int l= -1) const; //row of, efficient
@ -95,6 +96,7 @@ public:
const T* diagonalof(NRVec<T> &, const bool divide=0, bool cache=false) const; //get diagonal const T* diagonalof(NRVec<T> &, const bool divide=0, bool cache=false) const; //get diagonal
void diagonalset(const NRVec<T> &); //set diagonal elements void diagonalset(const NRVec<T> &); //set diagonal elements
void gemv(const T beta, NRVec<T> &r, const char trans, const T alpha, const NRVec<T> &x) const {r.gemv(beta,*this,trans,alpha,x);}; void gemv(const T beta, NRVec<T> &r, const char trans, const T alpha, const NRVec<T> &x) const {r.gemv(beta,*this,trans,alpha,x);};
void gemv(const T beta, NRVec<complex<T> > &r, const char trans, const T alpha, const NRVec<complex<T> > &x) const {r.gemv(beta,*this,trans,alpha,x);};
inline T* operator[](const int i); //subscripting: pointer to row i inline T* operator[](const int i); //subscripting: pointer to row i
inline const T* operator[](const int i) const; inline const T* operator[](const int i) const;
inline T& operator()(const int i, const int j); // (i,j) subscripts inline T& operator()(const int i, const int j); // (i,j) subscripts
@ -105,6 +107,7 @@ public:
void get(int fd, bool dimensions=1, bool transposed=false); void get(int fd, bool dimensions=1, bool transposed=false);
void put(int fd, bool dimensions=1, bool transposed=false) const; void put(int fd, bool dimensions=1, bool transposed=false) const;
void copyonwrite(); void copyonwrite();
void clear() {if(nn&&mm) {copyonwrite(); LA_traits<T>::clear((*this)[0],nn*mm);}}; //zero out
void resize(int n, int m); void resize(int n, int m);
inline operator T*(); //get a pointer to the data inline operator T*(); //get a pointer to the data
inline operator const T*() const; inline operator const T*() const;
@ -118,7 +121,7 @@ public:
const char transb, const T &alpha);//this = alpha*op( A )*op( B ) + beta*this const char transb, const T &alpha);//this = alpha*op( A )*op( B ) + beta*this
void fprintf(FILE *f, const char *format, const int modulo) const; void fprintf(FILE *f, const char *format, const int modulo) const;
void fscanf(FILE *f, const char *format); void fscanf(FILE *f, const char *format);
const double norm(const T scalar=(T)0) const; const typename LA_traits<T>::normtype norm(const T scalar=(T)0) const;
void axpy(const T alpha, const NRMat &x); // this += a*x void axpy(const T alpha, const NRMat &x); // this += a*x
inline const T amax() const; inline const T amax() const;
const T trace() const; const T trace() const;
@ -258,7 +261,7 @@ template <typename T>
inline T* NRMat<T>::operator[](const int i) inline T* NRMat<T>::operator[](const int i)
{ {
#ifdef DEBUG #ifdef DEBUG
if (*count != 1) laerror("Mat lval use of [] with count > 1"); if (_LA_count_check && *count != 1) laerror("Mat lval use of [] with count > 1");
if (i<0 || i>=nn) laerror("Mat [] out of range"); if (i<0 || i>=nn) laerror("Mat [] out of range");
if (!v) laerror("[] for unallocated Mat"); if (!v) laerror("[] for unallocated Mat");
#endif #endif
@ -287,7 +290,7 @@ template <typename T>
inline T & NRMat<T>::operator()(const int i, const int j) inline T & NRMat<T>::operator()(const int i, const int j)
{ {
#ifdef DEBUG #ifdef DEBUG
if (*count != 1) laerror("Mat lval use of (,) with count > 1"); if (_LA_count_check && *count != 1) laerror("Mat lval use of (,) with count > 1");
if (i<0 || i>=nn &&nn>0 || j<0 || j>=mm && mm>0) laerror("Mat (,) out of range"); if (i<0 || i>=nn &&nn>0 || j<0 || j>=mm && mm>0) laerror("Mat (,) out of range");
if (!v) laerror("(,) for unallocated Mat"); if (!v) laerror("(,) for unallocated Mat");
#endif #endif
@ -626,7 +629,7 @@ public:
inline T& operator() (const int i, const int j) inline T& operator() (const int i, const int j)
{ {
#ifdef DEBUG #ifdef DEBUG
if (*NRMat<T>::count != 1) laerror("Mat lval use of (,) with count > 1"); if (_LA_count_check && *NRMat<T>::count != 1) laerror("Mat lval use of (,) with count > 1");
if (i<1 || i>NRMat<T>::nn || j<1 || j>NRMat<T>::mm) laerror("Mat (,) out of range"); if (i<1 || i>NRMat<T>::nn || j<1 || j>NRMat<T>::mm) laerror("Mat (,) out of range");
if (!NRMat<T>::v) laerror("(,) for unallocated Mat"); if (!NRMat<T>::v) laerror("(,) for unallocated Mat");
#endif #endif

13
smat.h
View File

@ -38,12 +38,12 @@ public:
inline NRSMat(const T &a, const int n); //Initialize to constant inline NRSMat(const T &a, const int n); //Initialize to constant
inline NRSMat(const T *a, const int n); // Initialize to array inline NRSMat(const T *a, const int n); // Initialize to array
inline NRSMat(const NRSMat &rhs); // Copy constructor inline NRSMat(const NRSMat &rhs); // Copy constructor
NRSMat(const typename LA_traits_complex<T>::NRSMat_Noncomplex_type &rhs, bool imagpart=false); //construct complex from real
explicit NRSMat(const NRMat<T> &rhs); // symmetric part of general matrix explicit NRSMat(const NRMat<T> &rhs); // symmetric part of general matrix
explicit NRSMat(const NRVec<T> &rhs, const int n); //construct matrix from vector explicit NRSMat(const NRVec<T> &rhs, const int n); //construct matrix from vector
NRSMat & operator|=(const NRSMat &rhs); //assignment to a new copy NRSMat & operator|=(const NRSMat &rhs); //assignment to a new copy
NRSMat & operator=(const NRSMat &rhs); //assignment NRSMat & operator=(const NRSMat &rhs); //assignment
void clear() {LA_traits<T>::clear(v,NN2);}; //zero out void randomize(const typename LA_traits<T>::normtype &x);
void randomize(const T&x);
NRSMat & operator=(const T &a); //assign a to diagonal NRSMat & operator=(const T &a); //assign a to diagonal
const bool operator!=(const NRSMat &rhs) const {if(nn!=rhs.nn) return 1; return LA_traits<T>::gencmp(v,rhs.v,NN2);} //memcmp for scalars else elementwise const bool operator!=(const NRSMat &rhs) const {if(nn!=rhs.nn) return 1; return LA_traits<T>::gencmp(v,rhs.v,NN2);} //memcmp for scalars else elementwise
const bool operator==(const NRSMat &rhs) const {return !(*this != rhs);}; const bool operator==(const NRSMat &rhs) const {return !(*this != rhs);};
@ -66,8 +66,10 @@ public:
const T dot(const NRSMat &rhs) const; // Smat.Smat//@@@for complex do conjugate const T dot(const NRSMat &rhs) const; // Smat.Smat//@@@for complex do conjugate
const T dot(const NRVec<T> &rhs) const; //Smat(as vec).vec //@@@for complex do conjugate const T dot(const NRVec<T> &rhs) const; //Smat(as vec).vec //@@@for complex do conjugate
const NRVec<T> operator*(const NRVec<T> &rhs) const {NRVec<T> result(nn); result.gemv((T)0,*this,'n',(T)1,rhs); return result;}; // Mat * Vec const NRVec<T> operator*(const NRVec<T> &rhs) const {NRVec<T> result(nn); result.gemv((T)0,*this,'n',(T)1,rhs); return result;}; // Mat * Vec
const NRVec<complex<T> > operator*(const NRVec<complex<T> > &rhs) const {NRVec<complex<T> > result(nn); result.gemv((T)0,*this,'n',(T)1,rhs); return result;}; // Mat * Vec
const T* diagonalof(NRVec<T> &, const bool divide=0, bool cache=false) const; //get diagonal const T* diagonalof(NRVec<T> &, const bool divide=0, bool cache=false) const; //get diagonal
void gemv(const T beta, NRVec<T> &r, const char trans, const T alpha, const NRVec<T> &x) const {r.gemv(beta,*this,trans,alpha,x);}; void gemv(const T beta, NRVec<T> &r, const char trans, const T alpha, const NRVec<T> &x) const {r.gemv(beta,*this,trans,alpha,x);};
void gemv(const T beta, NRVec<complex<T> > &r, const char trans, const T alpha, const NRVec<complex<T> > &x) const {r.gemv(beta,*this,trans,alpha,x);};
inline const T& operator[](const int ij) const; inline const T& operator[](const int ij) const;
inline T& operator[](const int ij); inline T& operator[](const int ij);
inline const T& operator()(const int i, const int j) const; inline const T& operator()(const int i, const int j) const;
@ -76,13 +78,14 @@ public:
inline int ncols() const; inline int ncols() const;
inline int size() const; inline int size() const;
inline bool transp(const int i, const int j) const {return i>j;} //this can be used for compact storage of matrices, which are actually not symmetric, but one triangle of them is redundant inline bool transp(const int i, const int j) const {return i>j;} //this can be used for compact storage of matrices, which are actually not symmetric, but one triangle of them is redundant
const double norm(const T scalar=(T)0) const; const typename LA_traits<T>::normtype norm(const T scalar=(T)0) const;
void axpy(const T alpha, const NRSMat &x); // this+= a*x void axpy(const T alpha, const NRSMat &x); // this+= a*x
inline const T amax() const; inline const T amax() const;
const T trace() const; const T trace() const;
void get(int fd, bool dimensions=1, bool transp=0); void get(int fd, bool dimensions=1, bool transp=0);
void put(int fd, bool dimensions=1, bool transp=0) const; void put(int fd, bool dimensions=1, bool transp=0) const;
void copyonwrite(); void copyonwrite();
void clear() {copyonwrite(); LA_traits<T>::clear(v,NN2);}; //zero out
void resize(const int n); void resize(const int n);
inline operator T*(); //get a pointer to the data inline operator T*(); //get a pointer to the data
inline operator const T*() const; //get a pointer to the data inline operator const T*() const; //get a pointer to the data
@ -283,7 +286,7 @@ template <typename T>
inline T & NRSMat<T>::operator[](const int ij) inline T & NRSMat<T>::operator[](const int ij)
{ {
#ifdef DEBUG #ifdef DEBUG
if (*count != 1) laerror("lval [] with count > 1 in Smat"); if (_LA_count_check && *count != 1) laerror("lval [] with count > 1 in Smat");
if (ij<0 || ij>=NN2) laerror("SMat [] out of range"); if (ij<0 || ij>=NN2) laerror("SMat [] out of range");
if (!v) laerror("[] for unallocated Smat"); if (!v) laerror("[] for unallocated Smat");
#endif #endif
@ -359,7 +362,7 @@ template <typename T>
inline T & NRSMat<T>::operator()(const int i, const int j) inline T & NRSMat<T>::operator()(const int i, const int j)
{ {
#ifdef DEBUG #ifdef DEBUG
if (*count != 1) laerror("lval (i,j) with count > 1 in Smat"); if (_LA_count_check && *count != 1) laerror("lval (i,j) with count > 1 in Smat");
if (i<0 || i>=nn || j<0 || j>=nn) laerror("SMat (i,j) out of range"); if (i<0 || i>=nn || j<0 || j>=nn) laerror("SMat (i,j) out of range");
if (!v) laerror("(i,j) for unallocated Smat"); if (!v) laerror("(i,j) for unallocated Smat");
#endif #endif

2
sparsemat.h
View File

@ -88,7 +88,6 @@ public:
SparseMat & operator+=(const SparseMat &rhs); SparseMat & operator+=(const SparseMat &rhs);
SparseMat & addtriangle(const SparseMat &rhs, const bool lower, const char sign); SparseMat & addtriangle(const SparseMat &rhs, const bool lower, const char sign);
SparseMat & join(SparseMat &rhs); //more efficient +=, rhs will be emptied SparseMat & join(SparseMat &rhs); //more efficient +=, rhs will be emptied
void clear() {unsort();deletelist();}
SparseMat & operator-=(const SparseMat &rhs); SparseMat & operator-=(const SparseMat &rhs);
inline const SparseMat operator+(const T &rhs) const {return SparseMat(*this) += rhs;} inline const SparseMat operator+(const T &rhs) const {return SparseMat(*this) += rhs;}
inline const SparseMat operator-(const T &rhs) const {return SparseMat(*this) -= rhs;} inline const SparseMat operator-(const T &rhs) const {return SparseMat(*this) -= rhs;}
@ -133,6 +132,7 @@ public:
inline bool issymmetric() const {return symmetric;} inline bool issymmetric() const {return symmetric;}
unsigned int length() const; unsigned int length() const;
void copyonwrite(); void copyonwrite();
void clear() {copyonwrite();unsort();deletelist();}
void simplify(); void simplify();
const T trace() const; const T trace() const;
const T norm(const T scalar=(T)0) const; //is const only mathematically, not in internal implementation - we have to simplify first const T norm(const T scalar=(T)0) const; //is const only mathematically, not in internal implementation - we have to simplify first

11
vec.h
View File

@ -63,6 +63,7 @@ public:
inline NRVec(const T *a, const int n); inline NRVec(const T *a, const int n);
inline NRVec(T *a, const int n, bool skeleton); inline NRVec(T *a, const int n, bool skeleton);
inline NRVec(const NRVec &rhs); inline NRVec(const NRVec &rhs);
NRVec(const typename LA_traits_complex<T>::NRVec_Noncomplex_type &rhs, bool imagpart=false); //construct complex from real
inline explicit NRVec(const NRSMat<T> & S); inline explicit NRVec(const NRSMat<T> & S);
#ifdef MATPTR #ifdef MATPTR
explicit NRVec(const NRMat<T> &rhs) : NRVec(&rhs[0][0],rhs.nrows()*rhs.ncols()) {}; explicit NRVec(const NRMat<T> &rhs) : NRVec(&rhs[0][0],rhs.nrows()*rhs.ncols()) {};
@ -71,8 +72,7 @@ public:
#endif #endif
NRVec & operator=(const NRVec &rhs); NRVec & operator=(const NRVec &rhs);
NRVec & operator=(const T &a); //assign a to every element NRVec & operator=(const T &a); //assign a to every element
void clear() {LA_traits<T>::clear(v,nn);}; //zero out void randomize(const typename LA_traits<T>::normtype &x);
void randomize(const T &x);
NRVec & operator|=(const NRVec &rhs); NRVec & operator|=(const NRVec &rhs);
const bool operator!=(const NRVec &rhs) const {if(nn!=rhs.nn) return 1; return LA_traits<T>::gencmp(v,rhs.v,nn);} //memcmp for scalars else elementwise const bool operator!=(const NRVec &rhs) const {if(nn!=rhs.nn) return 1; return LA_traits<T>::gencmp(v,rhs.v,nn);} //memcmp for scalars else elementwise
const bool operator==(const NRVec &rhs) const {return !(*this != rhs);}; const bool operator==(const NRVec &rhs) const {return !(*this != rhs);};
@ -99,6 +99,8 @@ public:
void gemv(const T beta, const NRMat<T> &a, const char trans, const T alpha, const NRVec &x); void gemv(const T beta, const NRMat<T> &a, const char trans, const T alpha, const NRVec &x);
void gemv(const T beta, const NRSMat<T> &a, const char trans /*just for compatibility*/, const T alpha, const NRVec &x); void gemv(const T beta, const NRSMat<T> &a, const char trans /*just for compatibility*/, const T alpha, const NRVec &x);
void gemv(const T beta, const SparseMat<T> &a, const char trans, const T alpha, const NRVec &x,const bool treat_as_symmetric=false); void gemv(const T beta, const SparseMat<T> &a, const char trans, const T alpha, const NRVec &x,const bool treat_as_symmetric=false);
void gemv(const typename LA_traits_complex<T>::Component_type beta, const typename LA_traits_complex<T>::NRMat_Noncomplex_type &a, const char trans, const typename LA_traits_complex<T>::Component_type alpha, const NRVec &x);
void gemv(const typename LA_traits_complex<T>::Component_type beta, const typename LA_traits_complex<T>::NRSMat_Noncomplex_type &a, const char trans, const typename LA_traits_complex<T>::Component_type alpha, const NRVec &x);
const NRVec operator*(const NRMat<T> &mat) const {NRVec<T> result(mat.ncols()); result.gemv((T)0,mat,'t',(T)1,*this); return result;}; const NRVec operator*(const NRMat<T> &mat) const {NRVec<T> result(mat.ncols()); result.gemv((T)0,mat,'t',(T)1,*this); return result;};
const NRVec operator*(const NRSMat<T> &mat) const {NRVec<T> result(mat.ncols()); result.gemv((T)0,mat,'t',(T)1,*this); return result;}; const NRVec operator*(const NRSMat<T> &mat) const {NRVec<T> result(mat.ncols()); result.gemv((T)0,mat,'t',(T)1,*this); return result;};
const NRVec operator*(const SparseMat<T> &mat) const {NRVec<T> result(mat.ncols()); result.gemv((T)0,mat,'t',(T)1,*this); return result;}; const NRVec operator*(const SparseMat<T> &mat) const {NRVec<T> result(mat.ncols()); result.gemv((T)0,mat,'t',(T)1,*this); return result;};
@ -115,11 +117,12 @@ public:
void axpy(const T alpha, const NRVec &x); // this+= a*x void axpy(const T alpha, const NRVec &x); // this+= a*x
void axpy(const T alpha, const T *x, const int stride=1); // this+= a*x void axpy(const T alpha, const T *x, const int stride=1); // this+= a*x
void copyonwrite(); void copyonwrite();
void clear() {copyonwrite(); LA_traits<T>::clear(v,nn);}; //zero out
void resize(const int n); void resize(const int n);
void get(int fd, bool dimensions=1, bool transp=0); void get(int fd, bool dimensions=1, bool transp=0);
void put(int fd, bool dimensions=1, bool transp=0) const; void put(int fd, bool dimensions=1, bool transp=0) const;
NRVec & normalize(); NRVec & normalize();
inline const double norm() const; inline const typename LA_traits<T>::normtype norm() const;
inline const T amax() const; inline const T amax() const;
inline const NRVec unitvector() const; inline const NRVec unitvector() const;
void fprintf(FILE *f, const char *format, const int modulo) const; void fprintf(FILE *f, const char *format, const int modulo) const;
@ -474,7 +477,7 @@ template <typename T>
inline T & NRVec<T>::operator[](const int i) inline T & NRVec<T>::operator[](const int i)
{ {
#ifdef DEBUG #ifdef DEBUG
if(*count != 1) laerror("possible lval [] with count > 1"); if(_LA_count_check && *count != 1) laerror("possible lval [] with count > 1");
if(i < 0 || i >= nn) laerror("NRVec out of range"); if(i < 0 || i >= nn) laerror("NRVec out of range");
if(!v) laerror("[] on unallocated NRVec"); if(!v) laerror("[] on unallocated NRVec");
#endif #endif