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

*** empty log message ***

This commit is contained in:
jiri 2009-11-12 21:01:19 +00:00
parent f44662bdab
commit 7f7c4aa553
33 changed files with 457 additions and 309 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

5
auxstorage.h
View File

@ -17,6 +17,7 @@
*/ */
#ifndef _AUXSTORAGE_H_ #ifndef _AUXSTORAGE_H_
#define _AUXSTORAGE_H_ #define _AUXSTORAGE_H_
#include "vec.h" #include "vec.h"
#include "mat.h" #include "mat.h"
#include "smat.h" #include "smat.h"
@ -27,6 +28,8 @@
#include <fcntl.h> #include <fcntl.h>
#include <unistd.h> #include <unistd.h>
namespace LA {
//CAUTION: //CAUTION:
//it will not work if T is itself a class with dynamically allocated components //it will not work if T is itself a class with dynamically allocated components
//it cannot be implemented for SparseMat, which lacks fixed record length //it cannot be implemented for SparseMat, which lacks fixed record length
@ -136,6 +139,6 @@ if(0>write(fd,&x(0,0),recl)) {perror(""); laerror("write failed in AuxStorage");
} }
}//namespace
#endif #endif

5
bisection.h
View File

@ -17,6 +17,9 @@
*/ */
#ifndef _BISECTION_H #ifndef _BISECTION_H
#define _BISECTION_H #define _BISECTION_H
namespace LA {
//general bisection search between dm and hm //general bisection search between dm and hm
//returns dm-1 on failure, otherwise number between dm and hm //returns dm-1 on failure, otherwise number between dm and hm
//cmp returns 0 on equal, >0 if first > second argument //cmp returns 0 on equal, >0 if first > second argument
@ -75,5 +78,5 @@ int interpolation_find(INDEX dm0, INDEX high, const SUBJECT *x, const SUBJECT *b
} }
}//namespace
#endif #endif

7
bitvector.cc
View File

@ -18,14 +18,16 @@
#include "bitvector.h" #include "bitvector.h"
namespace LA {
//inefficient I/O operators //inefficient I/O operators
ostream & operator<<(ostream &s, const bitvector &x) std::ostream & operator<<(std::ostream &s, const bitvector &x)
{ {
for(unsigned int i=0; i<x.size(); ++i) s<< x[i]; for(unsigned int i=0; i<x.size(); ++i) s<< x[i];
return s; return s;
} }
istream & operator>>(istream &s, bitvector &x) std::istream & operator>>(std::istream &s, bitvector &x)
{ {
bool a; bool a;
x.copyonwrite(); x.copyonwrite();
@ -168,3 +170,4 @@ if(modulo)
return s+word_popul(a); return s+word_popul(a);
} }
}//namespace

7
bitvector.h
View File

@ -21,6 +21,8 @@
#include "vec.h" #include "vec.h"
namespace LA {
//compressed storage of large bit vectors //compressed storage of large bit vectors
//any reasonable compiler changes the dividions and modulos to shift/and instructions //any reasonable compiler changes the dividions and modulos to shift/and instructions
@ -75,7 +77,8 @@ public:
}; };
extern ostream & operator<<(ostream &s, const bitvector &x); extern std::ostream & operator<<(std::ostream &s, const bitvector &x);
extern istream & operator>>(istream &s, bitvector &x); extern std::istream & operator>>(std::istream &s, bitvector &x);
}//namespace
#endif #endif

14
conjgrad.h
View File

@ -15,6 +15,8 @@
You should have received a copy of the GNU General Public License You should have received a copy of the GNU General Public License
along with this program. If not, see <http://www.gnu.org/licenses/>. along with this program. If not, see <http://www.gnu.org/licenses/>.
*/ */
#ifndef _CONJGRAD_H_
#define _CONJGRAD_H_
#include "vec.h" #include "vec.h"
#include "smat.h" #include "smat.h"
#include "mat.h" #include "mat.h"
@ -22,6 +24,8 @@
#include "nonclass.h" #include "nonclass.h"
#include <iomanip> #include <iomanip>
namespace LA {
//conjugate gradient solution of a linear system //conjugate gradient solution of a linear system
//matrix can be any class which has nrows(), ncols(), diagonalof() and gemv() available //matrix can be any class which has nrows(), ncols(), diagonalof() and gemv() available
@ -63,10 +67,10 @@ for(int iter=0; iter<= itmax; iter++)
double err=p.norm(); double err=p.norm();
if(verbose) if(verbose)
{ {
cout << "conjgrad: iter= "<<iter<<" err= "<< std::cout << "conjgrad: iter= "<<iter<<" err= "<<
setiosflags(ios::scientific)<<setprecision(8) <<err<< std::setiosflags(std::ios::scientific)<<std::setprecision(8) <<err<<
resetiosflags(ios::scientific)<<setprecision(12)<<"\n"; std::resetiosflags(std::ios::scientific)<<std::setprecision(12)<<"\n";
cout.flush(); std::cout.flush();
} }
if(err <= tol) if(err <= tol)
{ {
@ -94,3 +98,5 @@ if(!issquare) delete r;
return false; return false;
} }
}//namespace
#endif

25
davidson.h
View File

@ -24,6 +24,8 @@
#include "nonclass.h" #include "nonclass.h"
#include "auxstorage.h" #include "auxstorage.h"
namespace LA {
//Davidson diagonalization of real symmetric matrix (modified Lanczos) //Davidson diagonalization of real symmetric matrix (modified Lanczos)
//matrix can be any class which has nrows(), ncols(), diagonalof(), issymmetric(), and gemv() available //matrix can be any class which has nrows(), ncols(), diagonalof(), issymmetric(), and gemv() available
@ -146,21 +148,21 @@ else
T eival_n=r[nroot]; T eival_n=r[nroot];
oldnroot=nroot; oldnroot=nroot;
typename LA_traits<T>::normtype test=abs(smallV(krylovsize,nroot)); typename LA_traits<T>::normtype test=std::abs(smallV(krylovsize,nroot));
if(test<eps) nroot++; if(test<eps) nroot++;
if(it==0) nroot = 0; if(it==0) nroot = 0;
for(int iroot=0; iroot<=min(krylovsize,nroots-1); ++iroot) for(int iroot=0; iroot<=std::min(krylovsize,nroots-1); ++iroot)
{ {
test = abs(smallV(krylovsize,iroot)); test = std::abs(smallV(krylovsize,iroot));
if(test>eps) nroot=min(nroot,iroot); if(test>eps) nroot=std::min(nroot,iroot);
if(verbose && iroot<=max(oldnroot,nroot)) if(verbose && iroot<=std::max(oldnroot,nroot))
{ {
cout <<"Davidson: iter="<<it <<" dim="<<krylovsize<<" root="<<iroot<<" energy="<<r[iroot]<<"\n"; std::cout <<"Davidson: iter="<<it <<" dim="<<krylovsize<<" root="<<iroot<<" energy="<<r[iroot]<<"\n";
cout.flush(); std::cout.flush();
} }
} }
if(verbose && oldnroot!=nroot) {cout <<"root no. "<<oldnroot<<" converged\n"; cout.flush();} if(verbose && oldnroot!=nroot) {std::cout <<"root no. "<<oldnroot<<" converged\n"; std::cout.flush();}
if (nroot>=nroots) goto converged; if (nroot>=nroots) goto converged;
if (it==maxit-1) break; //not converged if (it==maxit-1) break; //not converged
@ -195,7 +197,7 @@ eival_n = r[nroot];
for(i=0; i<n; ++i) for(i=0; i<n; ++i)
{ {
T denom = diagonal[i] - eival_n; T denom = diagonal[i] - eival_n;
denom = denom<0?-max(0.1,abs(denom)):max(0.1,abs(denom)); denom = denom<0?-std::max(0.1,std::abs(denom)):std::max(0.1,std::abs(denom));
vec1[i] /= denom; vec1[i] /= denom;
} }
@ -224,7 +226,7 @@ goto finished;
converged: converged:
AuxStorage<typename LA_traits<T>::elementtype> *ev; AuxStorage<typename LA_traits<T>::elementtype> *ev;
if(eivecsfile) ev = new AuxStorage<typename LA_traits<T>::elementtype>(eivecsfile); if(eivecsfile) ev = new AuxStorage<typename LA_traits<T>::elementtype>(eivecsfile);
if(verbose) {cout << "Davidson converged in "<<it<<" iterations.\n"; cout.flush();} if(verbose) {std::cout << "Davidson converged in "<<it<<" iterations.\n"; std::cout.flush();}
for(nroot=0; nroot<nroots; ++nroot) for(nroot=0; nroot<nroots; ++nroot)
{ {
eivals[nroot]=r[nroot]; eivals[nroot]=r[nroot];
@ -254,6 +256,5 @@ else {delete s0; delete s1;}
if(flag) laerror("no convergence in davidson"); if(flag) laerror("no convergence in davidson");
} }
}//namespace
#endif #endif

7
diis.h
View File

@ -26,8 +26,10 @@
#include "la_traits.h" #include "la_traits.h"
#include "auxstorage.h" #include "auxstorage.h"
namespace LA {
//Pulay memorial book remarks - for numerical stabilization small addition to diagonal //Pulay memorial book remarks - for numerical stabilization small addition to diagonal
#define DIISEPS 1e-9 const double DIISEPS=1e-9;
// Typically, 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=(const T&), clear(), dot() , axpy(), norm() and copyonwrite(), and LA_traits<T>::normtype and elementtype // actually it can be anything what has operator=(const T&), clear(), dot() , axpy(), norm() and copyonwrite(), and LA_traits<T>::normtype and elementtype
@ -150,7 +152,7 @@ rhs= (Ue)0; rhs[0]= (Ue)-1;
NRSMat<Te> amat=bmat; NRSMat<Te> amat=bmat;
linear_solve(amat,rhs,NULL,aktdim+1); linear_solve(amat,rhs,NULL,aktdim+1);
} }
if(verbose) cout <<"DIIS coefficients: "<<rhs<<endl; if(verbose) std::cout <<"DIIS coefficients: "<<rhs<<std::endl;
//build the new linear combination //build the new linear combination
vec.clear(); vec.clear();
@ -170,5 +172,6 @@ else
return norm; return norm;
} }
}//namespace
#endif #endif

28
fourindex.h
View File

@ -26,9 +26,13 @@
#include <sys/stat.h> #include <sys/stat.h>
#include <unistd.h> #include <unistd.h>
#include <sys/stat.h> #include <sys/stat.h>
#include "la.h"
#include "laerror.h" #include "laerror.h"
#include "vec.h"
#include "smat.h"
#include "mat.h"
#include "nonclass.h"
namespace LA {
static unsigned int hcd0(unsigned int big,unsigned int small) static unsigned int hcd0(unsigned int big,unsigned int small)
{ {
@ -591,7 +595,7 @@ return n;
} }
template <class I, class T> template <class I, class T>
ostream& operator<<(ostream &s, const fourindex_ext<I,T> &x) std::ostream& operator<<(std::ostream &s, const fourindex_ext<I,T> &x)
{ {
int n; int n;
n=x.size(); n=x.size();
@ -609,7 +613,7 @@ ostream& operator<<(ostream &s, const fourindex_ext<I,T> &x)
template <class I, class T> template <class I, class T>
ostream& operator<<(ostream &s, const fourindex<I,T> &x) std::ostream& operator<<(std::ostream &s, const fourindex<I,T> &x)
{ {
int n; int n;
n=x.size(); n=x.size();
@ -625,7 +629,7 @@ ostream& operator<<(ostream &s, const fourindex<I,T> &x)
} }
template <class I, class T> template <class I, class T>
istream& operator>>(istream &s, fourindex<I,T> &x) std::istream& operator>>(std::istream &s, fourindex<I,T> &x)
{ {
typename LA_traits_io<I>::IOtype i,j,k,l; typename LA_traits_io<I>::IOtype i,j,k,l;
typename LA_traits_io<T>::IOtype elem; typename LA_traits_io<T>::IOtype elem;
@ -680,7 +684,7 @@ public:
const T& operator() (unsigned int i, unsigned int j, unsigned int k, unsigned int l) const; const T& operator() (unsigned int i, unsigned int j, unsigned int k, unsigned int l) const;
void resize(const int n) {(*this).NRSMat<T>::resize(n*(n+1)/2);}; void resize(const int n) {(*this).NRSMat<T>::resize(n*(n+1)/2);};
void putext(int f, T thr=1e-15); void putext(int f, T thr=1e-15);
int nbas() const {return (int)sqrt(2*(*this).nrows());}; int nbas() const {return (int)std::sqrt(2*(*this).nrows());};
}; };
@ -795,10 +799,10 @@ if (!NRMat<T>::v) laerror("access to unallocated fourindex_dense");
return (*this).NRMat<T>::operator() ((j-1)*noca+i-1,(b-1)*nvra+a-1); return (*this).NRMat<T>::operator() ((j-1)*noca+i-1,(b-1)*nvra+a-1);
} }
void print(ostream &out) const void print(std::ostream &out) const
{ {
unsigned int i,j,a,b; unsigned int i,j,a,b;
for(i=1; i<=noca; ++i) for(j=1; j<=nocb; ++j) for(a=1; a<=nvra; ++a) for(b=1; b<=nvrb; ++b) out << i<<" "<<j<<" "<<a<<" "<<b<<" "<<(*this)(i,j,a,b)<<endl; for(i=1; i<=noca; ++i) for(j=1; j<=nocb; ++j) for(a=1; a<=nvra; ++a) for(b=1; b<=nvrb; ++b) out << i<<" "<<j<<" "<<a<<" "<<b<<" "<<(*this)(i,j,a,b)<<std::endl;
} }
}; };
@ -854,10 +858,10 @@ if(a<b) {minus++; unsigned int t=a; a=b; b=t;}
} }
void print(ostream &out) const void print(std::ostream &out) const
{ {
unsigned int i,j,a,b; unsigned int i,j,a,b;
for(i=1; i<=nocc; ++i) for(j=1; j<i; ++j) for(a=1; a<=nvrt; ++a) for(b=1; b<a; ++b) out << i<<" "<<j<<" "<<a<<" "<<b<<" "<<(*this)(i,j,a,b)<<endl; for(i=1; i<=nocc; ++i) for(j=1; j<i; ++j) for(a=1; a<=nvrt; ++a) for(b=1; b<a; ++b) out << i<<" "<<j<<" "<<a<<" "<<b<<" "<<(*this)(i,j,a,b)<<std::endl;
} }
@ -886,14 +890,14 @@ public:
void add_unique(unsigned int i, unsigned int j, unsigned int k, unsigned int l, T elem); void add_unique(unsigned int i, unsigned int j, unsigned int k, unsigned int l, T elem);
const T& operator() (unsigned int i, unsigned int j, unsigned int k, unsigned int l) const; const T& operator() (unsigned int i, unsigned int j, unsigned int k, unsigned int l) const;
void resize(const int n) {nbas=n; (*this).NRSMat<T>::resize(n*(n-1)/2);}; void resize(const int n) {nbas=n; (*this).NRSMat<T>::resize(n*(n-1)/2);};
void print(ostream &out) const void print(std::ostream &out) const
{ {
unsigned int i,j,k,l; unsigned int i,j,k,l;
for(i=1; i<=nbas; ++i) for(i=1; i<=nbas; ++i)
for(k=1;k<i; ++k) for(k=1;k<i; ++k)
for(j=1; j<=i; ++j) for(j=1; j<=i; ++j)
for(l=1; l<j && (j==i ? l<=k : 1); ++l) for(l=1; l<j && (j==i ? l<=k : 1); ++l)
cout << i<<" "<<k<<" "<<j<<" "<<l<<" "<<(*this)(i,k,j,l)<<endl; std::cout << i<<" "<<k<<" "<<j<<" "<<l<<" "<<(*this)(i,k,j,l)<<std::endl;
} }
@ -992,6 +996,6 @@ for(p= const_cast<fourindex_ext<I,T> *>(&rhs)->pbegin(); p.notend(); ++p)
} }
}//namespace
#endif /*_fourindex_included*/ #endif /*_fourindex_included*/

21
gmres.h
View File

@ -16,6 +16,8 @@
You should have received a copy of the GNU General Public License You should have received a copy of the GNU General Public License
along with this program. If not, see <http://www.gnu.org/licenses/>. along with this program. If not, see <http://www.gnu.org/licenses/>.
*/ */
#ifndef _GMRES_H
#define _GMRES_H
#include "vec.h" #include "vec.h"
#include "smat.h" #include "smat.h"
#include "mat.h" #include "mat.h"
@ -24,6 +26,8 @@
#include <iomanip> #include <iomanip>
#include "auxstorage.h" #include "auxstorage.h"
namespace LA {
//GMRES solution of a linear system //GMRES solution of a linear system
//matrix can be any class which has nrows(), ncols(), diagonalof() and gemv() available //matrix can be any class which has nrows(), ncols(), diagonalof() and gemv() available
@ -135,7 +139,7 @@ for (int l=0;l<neustart;l++) // main loop for restarts
H(k+1,k) = tmp= z.norm(); H(k+1,k) = tmp= z.norm();
if(tmp < 1.e-2*eps ) if(tmp < 1.e-2*eps )
{ {
if(verbose) cerr <<("gmres restart performed\n"); if(verbose) std::cerr <<("gmres restart performed\n");
// Abbruchbedingung, konstruiere x_k // Abbruchbedingung, konstruiere x_k
for (int i=0;i<k;i++) for (int i=0;i<k;i++)
{ {
@ -190,12 +194,12 @@ for (int l=0;l<neustart;l++) // main loop for restarts
//Schritt 6: Konvergenz? //Schritt 6: Konvergenz?
if(verbose) if(verbose)
{ {
cout << "gmres iter "<<l<<" "<<k<<" resid " std::cout << "gmres iter "<<l<<" "<<k<<" resid "
<<setw(0)<<setiosflags(ios::scientific)<<setprecision(8) <<std::setw(0)<<std::setiosflags(std::ios::scientific)<<std::setprecision(8)
<<abs(d[k+1])<< " thr "<<eps*beta_0<< " reduction " <<std::abs(d[k+1])<< " thr "<<eps*beta_0<< " reduction "
<<setw(5)<<setprecision(2)<<resetiosflags(ios::scientific) <<std::setw(5)<<std::setprecision(2)<<std::resetiosflags(std::ios::scientific)
<<(d_alt - abs(d[k+1]))/d_alt*100<< "\n" <<setprecision(12); <<(d_alt - std::abs(d[k+1]))/d_alt*100<< "\n" <<std::setprecision(12);
cout.flush(); std::cout.flush();
} }
d_alt = abs(d[k+1]); d_alt = abs(d[k+1]);
@ -258,3 +262,6 @@ if(!incore) delete st;
return !flag; return !flag;
} }
}//namespace
#endif

25
la.h
View File

@ -20,10 +20,29 @@
//this should be the single include file for the end user //this should be the single include file for the end user
// //
#ifdef USE_TRACEBACK
#include "vec.h" #include "traceback.h"
#include "smat.h" #endif
#include "la_traits.h"
#include "laerror.h"
#include "auxstorage.h"
#include "bisection.h"
#include "bitvector.h"
#include "conjgrad.h"
#include "davidson.h"
#include "diis.h"
#include "fourindex.h"
#include "gmres.h"
#include "mat.h" #include "mat.h"
#include "matexp.h"
#include "noncblas.h"
#include "nonclass.h" #include "nonclass.h"
#include "permutation.h"
#include "qsort.h"
#include "smat.h"
#include "sparsemat.h"
#include "sparsesmat.h"
#include "vec.h"
#endif /* _LA_H_ */ #endif /* _LA_H_ */

28
la_traits.h
View File

@ -27,14 +27,17 @@
#define _LA_TRAITS_INCL #define _LA_TRAITS_INCL
extern bool _LA_count_check;
using namespace std;
#include <stdio.h> #include <stdio.h>
#include <string.h> #include <string.h>
#include <iostream> #include <iostream>
#include <complex> #include <complex>
//using namespace std;
#define complex std::complex
#include "laerror.h" #include "laerror.h"
#ifdef NONCBLAS #ifdef NONCBLAS
@ -53,6 +56,9 @@ extern "C" {
} }
#endif #endif
namespace LA {
extern bool _LA_count_check;
//forward declarations //forward declarations
template<typename C> class NRVec; template<typename C> class NRVec;
@ -61,6 +67,7 @@ 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 NRSMat_from1;
template<typename C> class SparseMat; template<typename C> class SparseMat;
template<typename C> class SparseSMat;
typedef class {} Dummy_type; typedef class {} Dummy_type;
@ -95,9 +102,13 @@ SPECIALIZE_COMPLEX(complex<float>)
SPECIALIZE_COMPLEX(char) SPECIALIZE_COMPLEX(char)
SPECIALIZE_COMPLEX(unsigned char) SPECIALIZE_COMPLEX(unsigned char)
SPECIALIZE_COMPLEX(short) SPECIALIZE_COMPLEX(short)
SPECIALIZE_COMPLEX(unsigned short)
SPECIALIZE_COMPLEX(int) SPECIALIZE_COMPLEX(int)
SPECIALIZE_COMPLEX(unsigned int) SPECIALIZE_COMPLEX(unsigned int)
SPECIALIZE_COMPLEX(long)
SPECIALIZE_COMPLEX(unsigned long) SPECIALIZE_COMPLEX(unsigned long)
SPECIALIZE_COMPLEX(long long)
SPECIALIZE_COMPLEX(unsigned long long)
//for general sortable classes //for general sortable classes
@ -194,10 +205,11 @@ struct LA_traits_aux<complex<C>, scalar_true> {
typedef complex<C> elementtype; typedef complex<C> elementtype;
typedef complex<C> producttype; typedef complex<C> producttype;
typedef C normtype; typedef C normtype;
static inline C sqrabs(const complex<C> x) { return x.real()*x.real()+x.imag()*x.imag();}
static inline bool gencmp(const complex<C> *x, const complex<C> *y, int n) {return memcmp(x,y,n*sizeof(complex<C>));} static inline bool gencmp(const complex<C> *x, const complex<C> *y, int n) {return memcmp(x,y,n*sizeof(complex<C>));}
static bool bigger(const complex<C> &x, const complex<C> &y) {laerror("complex comparison undefined"); return false;} static bool bigger(const complex<C> &x, const complex<C> &y) {laerror("complex comparison undefined"); return false;}
static bool smaller(const complex<C> &x, const complex<C> &y) {laerror("complex comparison undefined"); return false;} static bool smaller(const complex<C> &x, const complex<C> &y) {laerror("complex comparison undefined"); return false;}
static inline normtype norm (const complex<C> &x) {return abs(x);} static inline normtype norm (const complex<C> &x) {return std::abs(x);}
static inline void axpy (complex<C> &s, const complex<C> &x, const complex<C> &c) {s+=x*c;} static inline void axpy (complex<C> &s, const complex<C> &x, const complex<C> &c) {s+=x*c;}
static inline void get(int fd, complex<C> &x, bool dimensions=0, bool transp=0) {if(sizeof(complex<C>)!=read(fd,&x,sizeof(complex<C>))) laerror("read error");} static inline void get(int fd, complex<C> &x, bool dimensions=0, bool transp=0) {if(sizeof(complex<C>)!=read(fd,&x,sizeof(complex<C>))) laerror("read error");}
static inline void put(int fd, const complex<C> &x, bool dimensions=0, bool transp=0) {if(sizeof(complex<C>)!=write(fd,&x,sizeof(complex<C>))) laerror("write error");} static inline void put(int fd, const complex<C> &x, bool dimensions=0, bool transp=0) {if(sizeof(complex<C>)!=write(fd,&x,sizeof(complex<C>))) laerror("write error");}
@ -206,6 +218,7 @@ static void multiput(unsigned int n, int fd, const complex<C> *x, bool dimension
static void copy(complex<C> *dest, complex<C> *src, unsigned int n) {memcpy(dest,src,n*sizeof(complex<C>));} static void copy(complex<C> *dest, complex<C> *src, unsigned int n) {memcpy(dest,src,n*sizeof(complex<C>));}
static void clear(complex<C> *dest, unsigned int n) {memset(dest,0,n*sizeof(complex<C>));} static void clear(complex<C> *dest, unsigned int n) {memset(dest,0,n*sizeof(complex<C>));}
static void copyonwrite(complex<C> &x) {}; static void copyonwrite(complex<C> &x) {};
static void clearme(complex<C> &x) {x=0;};
}; };
//non-complex scalars //non-complex scalars
@ -214,10 +227,11 @@ struct LA_traits_aux<C, scalar_true> {
typedef C elementtype; typedef C elementtype;
typedef C producttype; typedef C producttype;
typedef C normtype; typedef C normtype;
static inline C sqrabs(const C x) { return x*x;}
static inline bool gencmp(const C *x, const C *y, int n) {return memcmp(x,y,n*sizeof(C));} static inline bool gencmp(const C *x, const C *y, int n) {return memcmp(x,y,n*sizeof(C));}
static inline bool bigger(const C &x, const C &y) {return x>y;} static inline bool bigger(const C &x, const C &y) {return x>y;}
static inline bool smaller(const C &x, const C &y) {return x<y;} static inline bool smaller(const C &x, const C &y) {return x<y;}
static inline normtype norm (const C &x) {return abs(x);} static inline normtype norm (const C &x) {return std::abs(x);}
static inline void axpy (C &s, const C &x, const C &c) {s+=x*c;} static inline void axpy (C &s, const C &x, const C &c) {s+=x*c;}
static inline void put(int fd, const C &x, bool dimensions=0, bool transp=0) {if(sizeof(C)!=write(fd,&x,sizeof(C))) laerror("write error");} static inline void put(int fd, const C &x, bool dimensions=0, bool transp=0) {if(sizeof(C)!=write(fd,&x,sizeof(C))) laerror("write error");}
static inline void get(int fd, C &x, bool dimensions=0, bool transp=0) {if(sizeof(C)!=read(fd,&x,sizeof(C))) laerror("read error");} static inline void get(int fd, C &x, bool dimensions=0, bool transp=0) {if(sizeof(C)!=read(fd,&x,sizeof(C))) laerror("read error");}
@ -226,6 +240,7 @@ static void multiget(unsigned int n, int fd, C *x, bool dimensions=0) {if((ssize
static void copy(C *dest, C *src, unsigned int n) {memcpy(dest,src,n*sizeof(C));} static void copy(C *dest, C *src, unsigned int n) {memcpy(dest,src,n*sizeof(C));}
static void clear(C *dest, unsigned int n) {memset(dest,0,n*sizeof(C));} static void clear(C *dest, unsigned int n) {memset(dest,0,n*sizeof(C));}
static void copyonwrite(C &x) {}; static void copyonwrite(C &x) {};
static void clearme(complex<C> &x) {x=0;};
}; };
@ -252,6 +267,7 @@ static void multiget(unsigned int n,int fd, C *x, bool dimensions=1) {for(unsign
static void copy(C *dest, C *src, unsigned int n) {for(unsigned int i=0; i<n; ++i) dest[i]=src[i];} \ 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();}\ static void clear(C *dest, unsigned int n) {for(unsigned int i=0; i<n; ++i) dest[i].clear();}\
static void copyonwrite(X<C> &x) {x.copyonwrite();}\ static void copyonwrite(X<C> &x) {x.copyonwrite();}\
static void clearme(X<C> &x) {x.clear();}\
}; };
@ -260,6 +276,7 @@ generate_traits(NRMat)
generate_traits(NRMat_from1) generate_traits(NRMat_from1)
generate_traits(NRVec) generate_traits(NRVec)
generate_traits(SparseMat) generate_traits(SparseMat)
generate_traits(SparseSMat) //product leading to non-symmetric result not implemented
#undef generate_traits #undef generate_traits
@ -282,6 +299,7 @@ static void multiget(unsigned int n,int fd, C *x, bool dimensions=1) {for(unsign
static void copy(C *dest, C *src, unsigned int n) {for(unsigned int i=0; i<n; ++i) dest[i]=src[i];} \ 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();} \ static void clear(C *dest, unsigned int n) {for(unsigned int i=0; i<n; ++i) dest[i].clear();} \
static void copyonwrite(X<C> &x) {x.copyonwrite();} \ static void copyonwrite(X<C> &x) {x.copyonwrite();} \
static void clearme(X<C> &x) {x.clear();} \
}; };
generate_traits_smat(NRSMat) generate_traits_smat(NRSMat)
@ -292,4 +310,6 @@ generate_traits_smat(NRSMat_from1)
template<typename C> template<typename C>
struct LA_traits : LA_traits_aux<C, typename isscalar<C>::scalar_type> {}; struct LA_traits : LA_traits_aux<C, typename isscalar<C>::scalar_type> {};
}//namespace
#endif #endif

4
laerror.cc
View File

@ -30,6 +30,8 @@
#include "traceback.h" #include "traceback.h"
#endif #endif
namespace LA {
bool _LA_count_check=true; bool _LA_count_check=true;
extern "C" void _findme(void) {}; //for autoconf test we need a function with C linkage extern "C" void _findme(void) {}; //for autoconf test we need a function with C linkage
@ -86,3 +88,5 @@ int cblas_errprn(int ierr, int info, char *form, ...)
laerror(msg0); laerror(msg0);
return 0; return 0;
} }
}//namespace

4
laerror.h
View File

@ -17,9 +17,10 @@
*/ */
#ifndef _LAERROR_H_ #ifndef _LAERROR_H_
#define _LAERROR_H_ #define _LAERROR_H_
#include <iostream> #include <iostream>
namespace LA {
//exception class for laerror //exception class for laerror
class LAerror class LAerror
{ {
@ -36,5 +37,6 @@ s << x.msg;
return s; return s;
} }
}//namespace
#endif #endif

52
mat.cc
View File

@ -31,6 +31,7 @@ extern ssize_t write(int, const void *, size_t);
// TODO : // TODO :
// //
namespace LA {
/* /*
* Templates first, specializations for BLAS next * Templates first, specializations for BLAS next
@ -439,9 +440,9 @@ NRSMat<complex<double> > r(mm,mm);
int i,j; int i,j;
for(i=0; i<mm; ++i) for(j=0; j<=i; ++j) for(i=0; i<mm; ++i) for(j=0; j<=i; ++j)
#ifdef MATPTR #ifdef MATPTR
cblas_zdotc_sub(nn, v[0]+i , mm,v[0]+j, mm, (void *)(&r(i,j))); cblas_zdotc_sub(nn, v[0]+i , mm,v[0]+j, mm, &r(i,j));
#else #else
cblas_zdotc_sub(nn, v+i , mm,v+j, mm, (void *)(&r(i,j))); cblas_zdotc_sub(nn, v+i , mm,v+j, mm, &r(i,j));
#endif #endif
return r; return r;
} }
@ -486,9 +487,9 @@ NRSMat<complex<double> > r(nn,nn);
int i,j; int i,j;
for(i=0; i<nn; ++i) for(j=0; j<=i; ++j) for(i=0; i<nn; ++i) for(j=0; j<=i; ++j)
#ifdef MATPTR #ifdef MATPTR
cblas_zdotc_sub(mm, v[i] , 1,v[j], 1, (void *)(&r(i,j))); cblas_zdotc_sub(mm, v[i] , 1,v[j], 1, &r(i,j));
#else #else
cblas_zdotc_sub(mm, v+i*mm , 1,v+j*mm, 1, (void *)(&r(i,j))); cblas_zdotc_sub(mm, v+i*mm , 1,v+j*mm, 1, &r(i,j));
#endif #endif
return r; return r;
} }
@ -548,7 +549,7 @@ NRMat< complex<double> > &
NRMat< complex<double> >::operator*=(const complex<double> &a) NRMat< complex<double> >::operator*=(const complex<double> &a)
{ {
copyonwrite(); copyonwrite();
cblas_zscal(nn*mm, &a, (void *)(*this)[0], 1); cblas_zscal(nn*mm, &a, (*this)[0], 1);
return *this; return *this;
} }
@ -591,7 +592,7 @@ NRMat< complex<double> >::operator+=(const NRMat< complex<double> > &rhs)
laerror("Mat += Mat of incompatible matrices"); laerror("Mat += Mat of incompatible matrices");
#endif #endif
copyonwrite(); copyonwrite();
cblas_zaxpy(nn*mm, &CONE, (void *)rhs[0], 1, (void *)(*this)[0], 1); cblas_zaxpy(nn*mm, &CONE, rhs[0], 1, (*this)[0], 1);
return *this; return *this;
} }
@ -639,7 +640,7 @@ NRMat< complex<double> >::operator-=(const NRMat< complex<double> > &rhs)
laerror("Mat -= Mat of incompatible matrices"); laerror("Mat -= Mat of incompatible matrices");
#endif #endif
copyonwrite(); copyonwrite();
cblas_zaxpy(nn*mm, &CMONE, (void *)rhs[0], 1, (void *)(*this)[0], 1); cblas_zaxpy(nn*mm, &CMONE, rhs[0], 1, (*this)[0], 1);
return *this; return *this;
} }
@ -696,12 +697,12 @@ NRMat< complex<double> >::operator+=(const NRSMat< complex<double> > &rhs)
const complex<double> *p = rhs; const complex<double> *p = rhs;
copyonwrite(); copyonwrite();
for (int i=0; i<nn; i++) { for (int i=0; i<nn; i++) {
cblas_zaxpy(i+1, (void *)&CONE, (void *)p, 1, (void *)(*this)[i], 1); cblas_zaxpy(i+1, &CONE, p, 1, (*this)[i], 1);
p += i+1; p += i+1;
} }
p = rhs; p++; p = rhs; p++;
for (int i=1; i<nn; i++) { for (int i=1; i<nn; i++) {
cblas_zaxpy(i, (void *)&CONE, (void *)p, 1, (void *)((*this)[i]+i), nn); cblas_zaxpy(i, &CONE, p, 1, (*this)[0]+i, nn);
p += i+1; p += i+1;
} }
return *this; return *this;
@ -725,7 +726,7 @@ NRMat<T> & NRMat<T>::operator+=(const NRSMat<T> &rhs)
} }
p = rhs; p++; p = rhs; p++;
for (int i=1; i<nn; i++) { for (int i=1; i<nn; i++) {
for(int j=0; j<i; ++j) *((*this)[i]+i+nn*j) += p[j]; for(int j=0; j<i; ++j) *((*this)[j]+i) += p[j];
p += i+1; p += i+1;
} }
return *this; return *this;
@ -765,12 +766,12 @@ NRMat< complex<double> >::operator-=(const NRSMat< complex<double> > &rhs)
const complex<double> *p = rhs; const complex<double> *p = rhs;
copyonwrite(); copyonwrite();
for (int i=0; i<nn; i++) { for (int i=0; i<nn; i++) {
cblas_zaxpy(i+1, (void *)&CMONE, (void *)p, 1, (void *)(*this)[i], 1); cblas_zaxpy(i+1, &CMONE, p, 1, (*this)[i], 1);
p += i+1; p += i+1;
} }
p = rhs; p++; p = rhs; p++;
for (int i=1; i<nn; i++) { for (int i=1; i<nn; i++) {
cblas_zaxpy(i, (void *)&CMONE, (void *)p, 1, (void *)((*this)[i]+i), nn); cblas_zaxpy(i, &CMONE, p, 1, (*this)[0]+i, nn);
p += i+1; p += i+1;
} }
return *this; return *this;
@ -792,7 +793,7 @@ NRMat<T> & NRMat<T>::operator-=(const NRSMat<T> &rhs)
} }
p = rhs; p++; p = rhs; p++;
for (int i=1; i<nn; i++) { for (int i=1; i<nn; i++) {
for(int j=0; j<i; ++j) *((*this)[i]+i+nn*j) -= p[j]; for(int j=0; j<i; ++j) *((*this)[j]+i) -= p[j];
p += i+1; p += i+1;
} }
return *this; return *this;
@ -821,8 +822,8 @@ NRMat< complex<double> >::dot(const NRMat< complex<double> > &rhs) const
if(nn!=rhs.nn || mm!= rhs.mm) laerror("Mat.Mat incompatible matrices"); if(nn!=rhs.nn || mm!= rhs.mm) laerror("Mat.Mat incompatible matrices");
#endif #endif
complex<double> dot; complex<double> dot;
cblas_zdotc_sub(nn*mm, (void *)(*this)[0], 1, (void *)rhs[0], 1, cblas_zdotc_sub(nn*mm, (*this)[0], 1, rhs[0], 1,
(void *)(&dot)); &dot);
return dot; return dot;
} }
@ -851,8 +852,8 @@ NRMat< complex<double> >::operator*(const NRMat< complex<double> > &rhs) const
#endif #endif
NRMat< complex<double> > result(nn, rhs.mm); NRMat< complex<double> > result(nn, rhs.mm);
cblas_zgemm(CblasRowMajor, CblasNoTrans, CblasNoTrans, nn, rhs.mm, mm, cblas_zgemm(CblasRowMajor, CblasNoTrans, CblasNoTrans, nn, rhs.mm, mm,
(const void *)(&CONE),(const void *)(*this)[0], mm, (const void *)rhs[0], &CONE,(*this)[0], mm, rhs[0],
rhs.mm, (const void *)(&CZERO), (void *)result[0], rhs.mm); rhs.mm, &CZERO, result[0], rhs.mm);
return result; return result;
} }
@ -936,8 +937,8 @@ NRMat< complex<double> >::operator*(const NRSMat< complex<double> > &rhs) const
#endif #endif
NRMat< complex<double> > result(nn, rhs.ncols()); NRMat< complex<double> > result(nn, rhs.ncols());
for (int i=0; i<nn; i++) for (int i=0; i<nn; i++)
cblas_zhpmv(CblasRowMajor, CblasLower, mm, (void *)&CONE, (void *)&rhs[0], cblas_zhpmv(CblasRowMajor, CblasLower, mm, &CONE, &rhs[0],
(void *)(*this)[i], 1, (void *)&CZERO, (void *)result[i], 1); (*this)[i], 1, &CZERO, result[i], 1);
return result; return result;
} }
@ -969,7 +970,7 @@ NRMat< complex<double> >::transpose(bool conj) const
{ {
NRMat< complex<double> > result(mm,nn); NRMat< complex<double> > result(mm,nn);
for (int i=0; i<nn; i++) for (int i=0; i<nn; i++)
cblas_zcopy(mm, (void *)(*this)[i], 1, (void *)(result[0]+i), nn); cblas_zcopy(mm, (*this)[i], 1, (result[0]+i), nn);
if (conj) cblas_dscal(mm*nn, -1.0, (double *)(result[0])+1, 2); if (conj) cblas_dscal(mm*nn, -1.0, (double *)(result[0])+1, 2);
return result; return result;
} }
@ -1039,7 +1040,7 @@ const double NRMat<double>::norm(const double scalar) const
if (i==j) tmp -= scalar; if (i==j) tmp -= scalar;
sum += tmp*tmp; sum += tmp*tmp;
} }
return sqrt(sum); return std::sqrt(sum);
} }
@ -1060,7 +1061,7 @@ const double NRMat< complex<double> >::norm(const complex<double> scalar) const
if (i==j) tmp -= scalar; if (i==j) tmp -= scalar;
sum += tmp.real()*tmp.real()+tmp.imag()*tmp.imag(); sum += tmp.real()*tmp.real()+tmp.imag()*tmp.imag();
} }
return sqrt(sum); return std::sqrt(sum);
} }
@ -1087,7 +1088,7 @@ void NRMat< complex<double> >::axpy(const complex<double> alpha,
if (nn!=mat.nn || mm!=mat.mm) laerror("zaxpy of incompatible matrices"); if (nn!=mat.nn || mm!=mat.mm) laerror("zaxpy of incompatible matrices");
#endif #endif
copyonwrite(); copyonwrite();
cblas_zaxpy(nn*mm, (void *)&alpha, mat, 1, (void *)(*this)[0], 1); cblas_zaxpy(nn*mm, &alpha, mat, 1, (*this)[0], 1);
} }
@ -1186,7 +1187,7 @@ if(rowcol) //vectors are rows
NRVec<double> tmp = *metric * (*this).row(j); NRVec<double> tmp = *metric * (*this).row(j);
double norm = cblas_ddot(mm,(*this)[j],1,tmp,1); double norm = cblas_ddot(mm,(*this)[j],1,tmp,1);
if(norm<=0.) laerror("zero vector in orthonormalize or nonpositive metric"); if(norm<=0.) laerror("zero vector in orthonormalize or nonpositive metric");
cblas_dscal(mm,1./sqrt(norm),(*this)[j],1); cblas_dscal(mm,1./std::sqrt(norm),(*this)[j],1);
} }
} }
else //vectors are columns else //vectors are columns
@ -1203,7 +1204,7 @@ else //vectors are columns
NRVec<double> tmp = *metric * (*this).column(j); NRVec<double> tmp = *metric * (*this).column(j);
double norm = cblas_ddot(nn,&(*this)[0][j],mm,tmp,1); double norm = cblas_ddot(nn,&(*this)[0][j],mm,tmp,1);
if(norm<=0.) laerror("zero vector in orthonormalize or nonpositive metric"); if(norm<=0.) laerror("zero vector in orthonormalize or nonpositive metric");
cblas_dscal(nn,1./sqrt(norm),&(*this)[0][j],mm); cblas_dscal(nn,1./std::sqrt(norm),&(*this)[0][j],mm);
} }
} }
} }
@ -1261,3 +1262,4 @@ template class NRMat<unsigned int>;
template class NRMat<unsigned long>; template class NRMat<unsigned long>;
template class NRMat<unsigned long long>; template class NRMat<unsigned long long>;
}//namespace

14
mat.h
View File

@ -19,9 +19,9 @@
*/ */
#ifndef _LA_MAT_H_ #ifndef _LA_MAT_H_
#define _LA_MAT_H_ #define _LA_MAT_H_
#include "la_traits.h" #include "la_traits.h"
namespace LA {
template <typename T> template <typename T>
class NRMat { class NRMat {
protected: protected:
@ -91,6 +91,7 @@ public:
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<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
void orthonormalize(const bool rowcol, const NRSMat<T> *metric=NULL);//orthonormalize (true - vectors are rows)
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
const NRVec<T> column(const int j, int l= -1) const {if(l<0) l=nn; NRVec<T> r(l); for(int i=0; i<l; ++i) r[i]= (*this)(i,j); return r;}; //column of, general but not very efficient const NRVec<T> column(const int j, int l= -1) const {if(l<0) l=nn; NRVec<T> r(l); for(int i=0; i<l; ++i) r[i]= (*this)(i,j); return r;}; //column of, general but not very efficient
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
@ -140,11 +141,13 @@ public:
#endif #endif
}; };
}//namespace
//due to mutual includes this has to be after full class declaration //due to mutual includes this has to be after full class declaration
#include "vec.h" #include "vec.h"
#include "smat.h" #include "smat.h"
#include "sparsemat.h" #include "sparsemat.h"
namespace LA {
// ctors // ctors
template <typename T> template <typename T>
NRMat<T>::NRMat(const int n, const int m) : nn(n), mm(m), count(new int) NRMat<T>::NRMat(const int n, const int m) : nn(n), mm(m), count(new int)
@ -374,9 +377,9 @@ template<>
inline const complex<double> NRMat< complex<double> >::amax() const inline const complex<double> NRMat< complex<double> >::amax() const
{ {
#ifdef MATPTR #ifdef MATPTR
return v[0][cblas_izamax(nn*mm, (void *)v[0], 1)]; return v[0][cblas_izamax(nn*mm, v[0], 1)];
#else #else
return v[cblas_izamax(nn*mm, (void *)v, 1)]; return v[cblas_izamax(nn*mm, v, 1)];
#endif #endif
} }
@ -576,7 +579,7 @@ return r;
// I/O // I/O
template <typename T> template <typename T>
ostream& operator<<(ostream &s, const NRMat<T> &x) std::ostream& operator<<(std::ostream &s, const NRMat<T> &x)
{ {
int i,j,n,m; int i,j,n,m;
n=x.nrows(); n=x.nrows();
@ -590,7 +593,7 @@ ostream& operator<<(ostream &s, const NRMat<T> &x)
} }
template <typename T> template <typename T>
istream& operator>>(istream &s, NRMat<T> &x) std::istream& operator>>(std::istream &s, NRMat<T> &x)
{ {
int i,j,n,m; int i,j,n,m;
s >> n >> m; s >> n >> m;
@ -655,4 +658,5 @@ NRVECMAT_OPER(Mat,*)
NRVECMAT_OPER2(Mat,+) NRVECMAT_OPER2(Mat,+)
NRVECMAT_OPER2(Mat,-) NRVECMAT_OPER2(Mat,-)
}//namespace
#endif /* _LA_MAT_H_ */ #endif /* _LA_MAT_H_ */

34
matexp.h
View File

@ -21,10 +21,11 @@
//of matrix-matrix multiplications on cost of additions and memory //of matrix-matrix multiplications on cost of additions and memory
// the polynom and exp routines will work on any type, for which traits class // the polynom and exp routines will work on any type, for which traits class
// is defined containing definition of an element type, norm and axpy operation // is defined containing definition of an element type, norm and axpy operation
#include "la_traits.h" #include "la_traits.h"
#include "laerror.h" #include "laerror.h"
#include <math.h>
namespace LA {
template<class T,class R> template<class T,class R>
const T polynom0(const T &x, const NRVec<R> &c) const T polynom0(const T &x, const NRVec<R> &c)
@ -40,6 +41,7 @@ else
z=x*c[order]; z=x*c[order];
for(i=order-1; i>=0; i--) for(i=order-1; i>=0; i--)
{ {
//std::cerr<<"TEST polynom0 "<<i<<'\n';
if(i<order-1) z=y*x; if(i<order-1) z=y*x;
y=z+c[i]; y=z+c[i];
} }
@ -138,7 +140,7 @@ for(int i=1; i<=n; ++i)
{ {
factor/=i; factor/=i;
z= z*t-t*z; z= z*t-t*z;
if(verbose) cerr << "BCH contribution at order "<<i<<" : "<<z.norm()*factor<<endl; if(verbose) std::cerr << "BCH contribution at order "<<i<<" : "<<z.norm()*factor<<std::endl;
result+= z*factor; result+= z*factor;
} }
return result; return result;
@ -169,10 +171,10 @@ return y;
inline int nextpow2(const double n) inline int nextpow2(const double n)
{ {
const double log2=log(2.); const double log2=std::log(2.);
if(n<=.75) return 0; //try to keep the taylor expansion short if(n<=.75) return 0; //try to keep the taylor expansion short
if(n<=1.) return 1; if(n<=1.) return 1;
return int(ceil(log(n)/log2-log(.75))); return int(std::ceil(std::log(n)/log2-std::log(.75)));
} }
//should better be computed by mathematica to have accurate last digits, perhaps chebyshev instead, see exp in glibc //should better be computed by mathematica to have accurate last digits, perhaps chebyshev instead, see exp in glibc
@ -208,10 +210,10 @@ template<class T, class C, class S>
NRVec<C> exp_aux(const T &x, int &power, int maxpower, int maxtaylor, S prescale) NRVec<C> exp_aux(const T &x, int &power, int maxpower, int maxtaylor, S prescale)
{ {
double mnorm= x.norm() * abs(prescale); double mnorm= x.norm() * std::abs(prescale);
power=nextpow2(mnorm); power=nextpow2(mnorm);
if(maxpower>=0 && power>maxpower) power=maxpower; if(maxpower>=0 && power>maxpower) power=maxpower;
double scale=exp(-log(2.)*power); double scale=std::exp(-std::log(2.)*power);
//find how long taylor expansion will be necessary //find how long taylor expansion will be necessary
@ -236,8 +238,8 @@ for(i=0,t=1.;i<=n;i++)
taylor2[i]=exptaylor[i]*t; taylor2[i]=exptaylor[i]*t;
t*=scale; t*=scale;
} }
//cout <<"TEST power, scale "<<power<<" "<<scale<<endl; //std::cout <<"TEST power, scale "<<power<<" "<<scale<<std::endl;
//cout <<"TEST taylor2 "<<taylor2<<endl; //std::cout <<"TEST taylor2 "<<taylor2<<std::endl;
return taylor2; return taylor2;
} }
@ -246,10 +248,10 @@ return taylor2;
template<class T, class C, class S> template<class T, class C, class S>
void sincos_aux(NRVec<C> &si, NRVec<C> &co, const T &x, int &power,int maxpower, int maxtaylor, const S prescale) void sincos_aux(NRVec<C> &si, NRVec<C> &co, const T &x, int &power,int maxpower, int maxtaylor, const S prescale)
{ {
double mnorm= x.norm() * abs(prescale); double mnorm= x.norm() * std::abs(prescale);
power=nextpow2(mnorm); power=nextpow2(mnorm);
if(maxpower>=0 && power>maxpower) power=maxpower; if(maxpower>=0 && power>maxpower) power=maxpower;
double scale=exp(-log(2.)*power); double scale=std::exp(-std::log(2.)*power);
//find how long taylor expansion will be necessary //find how long taylor expansion will be necessary
const double precision=1e-14; //further decreasing brings nothing const double precision=1e-14; //further decreasing brings nothing
@ -275,8 +277,8 @@ for(i=0,t=1.;i<=n;i++)
else co[i>>1] = exptaylor[i]* (i&2?-t:t); else co[i>>1] = exptaylor[i]* (i&2?-t:t);
t*=scale; t*=scale;
} }
//cout <<"TEST sin "<<si<<endl; //std::cout <<"TEST sin "<<si<<std::endl;
//cout <<"TEST cos "<<co<<endl; //std::cout <<"TEST cos "<<co<<std::endl;
} }
@ -290,6 +292,7 @@ int power;
//prepare the polynom of and effectively scale T //prepare the polynom of and effectively scale T
NRVec<typename LA_traits<T>::normtype> taylor2=exp_aux<T,typename LA_traits<T>::normtype,double>(x,power,maxpower,maxtaylor,1.); NRVec<typename LA_traits<T>::normtype> taylor2=exp_aux<T,typename LA_traits<T>::normtype,double>(x,power,maxpower,maxtaylor,1.);
//std::cerr <<"TEST power "<<power<<std::endl;
T r= horner?polynom0(x,taylor2):polynom(x,taylor2); T r= horner?polynom0(x,taylor2):polynom(x,taylor2);
//for accuracy summing from the smallest terms up would be better, but this is more efficient for matrices //for accuracy summing from the smallest terms up would be better, but this is more efficient for matrices
@ -309,6 +312,7 @@ int power;
NRVec<typename LA_traits<T>::normtype> taylors,taylorc; NRVec<typename LA_traits<T>::normtype> taylors,taylorc;
sincos_aux<T,typename LA_traits<T>::normtype>(taylors,taylorc,x,power,maxpower,maxtaylor,1.); sincos_aux<T,typename LA_traits<T>::normtype>(taylors,taylorc,x,power,maxpower,maxtaylor,1.);
//could we save something by computing both polynoms simultaneously? //could we save something by computing both polynoms simultaneously?
{ {
T x2 = x*x; T x2 = x*x;
@ -446,16 +450,14 @@ NRVec<typename LA_traits<V>::normtype> taylors,taylorc;
sincos_aux<M,typename LA_traits<V>::normtype>(taylors,taylorc,mat,power,maxpower,maxtaylor,scale); sincos_aux<M,typename LA_traits<V>::normtype>(taylors,taylorc,mat,power,maxpower,maxtaylor,scale);
if(taylors.size()!=taylorc.size()) laerror("internal error - same size of sin and cos expansions assumed"); if(taylors.size()!=taylorc.size()) laerror("internal error - same size of sin and cos expansions assumed");
//the actual computation and resursive "squaring" //the actual computation and resursive "squaring"
cout <<"TEST power "<<power<<endl; //std::cout <<"TEST power "<<power<<std::endl;
sincostimes_aux(mat,si,co,rhs,taylors,taylorc,transpose,scale,power); sincostimes_aux(mat,si,co,rhs,taylors,taylorc,transpose,scale,power);
return; return;
} }
//@@@ power series matrix logarithm? //@@@ power series matrix logarithm?
}//namespace
#endif #endif

4
noncblas.cc
View File

@ -21,6 +21,8 @@
#include "laerror.h" #include "laerror.h"
#include "mat.h" #include "mat.h"
namespace LA {
#ifdef FORTRAN_ #ifdef FORTRAN_
#define FORNAME(x) x##_ #define FORNAME(x) x##_
#else #else
@ -255,3 +257,5 @@ return INFO;
#endif #endif
}//namespace

63
nonclass.cc
View File

@ -24,6 +24,7 @@
#include "nonclass.h" #include "nonclass.h"
#include "qsort.h" #include "qsort.h"
namespace LA {
#ifdef FORTRAN_ #ifdef FORTRAN_
#define FORNAME(x) x##_ #define FORNAME(x) x##_
@ -39,9 +40,13 @@ INSTANTIZE(complex<double>)
INSTANTIZE(int) INSTANTIZE(int)
INSTANTIZE(short) INSTANTIZE(short)
INSTANTIZE(char) INSTANTIZE(char)
INSTANTIZE(long)
INSTANTIZE(long long)
INSTANTIZE(unsigned char) INSTANTIZE(unsigned char)
INSTANTIZE(unsigned long) INSTANTIZE(unsigned short)
INSTANTIZE(unsigned int) INSTANTIZE(unsigned int)
INSTANTIZE(unsigned long)
INSTANTIZE(unsigned long long)
#define EPSDET 1e-300 #define EPSDET 1e-300
@ -146,7 +151,7 @@ static void linear_solve_do(NRMat<double> &A, double *B, const int nrhs, const i
if (det && r==0) { if (det && r==0) {
*det = 1.; *det = 1.;
//take into account some numerical instabilities in dgesv for singular matrices //take into account some numerical instabilities in dgesv for singular matrices
for (int i=0; i<n; ++i) {double t=A[i][i]; if(!finite(t) || abs(t) < EPSDET ) {*det=0.; break;} else *det *=t;} for (int i=0; i<n; ++i) {double t=A[i][i]; if(!finite(t) || std::abs(t) < EPSDET ) {*det=0.; break;} else *det *=t;}
//change sign of det by parity of ipiv permutation //change sign of det by parity of ipiv permutation
if(*det) for (int i=0; i<n; ++i) if( if(*det) for (int i=0; i<n; ++i) if(
#ifdef NONCBLAS #ifdef NONCBLAS
@ -158,9 +163,9 @@ static void linear_solve_do(NRMat<double> &A, double *B, const int nrhs, const i
} }
if(det && r>0) *det = 0; if(det && r>0) *det = 0;
/* /*
cout <<"ipiv = "; std::cout <<"ipiv = ";
for (int i=0; i<n; ++i) cout <<ipiv[i]<<" "; for (int i=0; i<n; ++i) std::cout <<ipiv[i]<<" ";
cout <<endl; std::cout <<std::endl;
*/ */
delete [] ipiv; delete [] ipiv;
if (r>0 && B) laerror("singular matrix in lapack_gesv"); if (r>0 && B) laerror("singular matrix in lapack_gesv");
@ -202,7 +207,7 @@ static void linear_solve_do(NRSMat<double> &a, double *b, const int nrhs, const
} }
if (det && r == 0) { if (det && r == 0) {
*det = 1.; *det = 1.;
for (int i=1; i<n; i++) {double t=a(i,i); if(!finite(t) || abs(t) < EPSDET ) {*det=0.; break;} else *det *= t;} for (int i=1; i<n; i++) {double t=a(i,i); if(!finite(t) || std::abs(t) < EPSDET ) {*det=0.; break;} else *det *= t;}
//do not use ipiv, since the permutation matrix occurs twice in the decomposition and signs thus cancel (man dspsv) //do not use ipiv, since the permutation matrix occurs twice in the decomposition and signs thus cancel (man dspsv)
} }
if (det && r>0) *det = 0; if (det && r>0) *det = 0;
@ -452,6 +457,9 @@ void singular_decomposition(NRMat<double> &a, NRMat<double> *u, NRVec<double> &s
if (r > 0) laerror("convergence problem in gesvd() of ingular_decomposition()"); if (r > 0) laerror("convergence problem in gesvd() of ingular_decomposition()");
} }
//QR decomposition
//extern "C" void FORNAME(dgeqrf)(const int *M, const int *N, double *A, const int *LDA, double *TAU, double *WORK, int *LWORK, int *INFO);
extern "C" void FORNAME(dgeev)(const char *JOBVL, const char *JOBVR, const int *N, extern "C" void FORNAME(dgeev)(const char *JOBVL, const char *JOBVR, const int *N,
double *A, const int *LDA, double *WR, double *WI, double *VL, const int *LDVL, double *A, const int *LDA, double *WR, double *WI, double *VL, const int *LDVL,
@ -555,7 +563,7 @@ void gdiagonalize(NRMat<double> &a, NRVec<double> &wr, NRVec<double> &wi,
&ldvr, vl?vl[0]:(double *)0, &ldvl, work, &lwork, &r); &ldvr, vl?vl[0]:(double *)0, &ldvl, work, &lwork, &r);
delete[] work; delete[] work;
//cout <<"TEST dgeev\n"<<wr<<wi<<*vr<<*vl<<endl; //std::cout <<"TEST dgeev\n"<<wr<<wi<<*vr<<*vl<<std::endl;
if (r < 0) laerror("illegal argument in ggev/geev in gdiagonalize()"); if (r < 0) laerror("illegal argument in ggev/geev in gdiagonalize()");
if (r > 0) laerror("convergence problem in ggev/geev in gdiagonalize()"); if (r > 0) laerror("convergence problem in ggev/geev in gdiagonalize()");
@ -589,9 +597,9 @@ void gdiagonalize(NRMat<double> &a, NRVec<double> &wr, NRVec<double> &wi,
x=.5*t*s11; x=.5*t*s11;
y=.5*t*s12; y=.5*t*s12;
double alp,bet; double alp,bet;
t=.5*sqrt(t); t=.5*std::sqrt(t);
alp=sqrt(.5*(t+x)); alp=std::sqrt(.5*(t+x));
bet=sqrt(.5*(t-x)); bet=std::sqrt(.5*(t-x));
if(y<0.) bet= -bet; if(y<0.) bet= -bet;
//rotate left ev //rotate left ev
@ -722,6 +730,16 @@ const NRMat< complex<double> > imagmatrix (const NRMat<double> &a)
return result; return result;
} }
const NRMat< complex<double> > complexmatrix (const NRMat<double> &re, const NRMat<double> &im)
{
if(re.nrows()!=im.nrows() || re.ncols() != im.ncols()) laerror("incompatible sizes of real and imaginary parts");
NRMat< complex<double> > result(re.nrows(), re.ncols());
cblas_dcopy(re.nrows()*re.ncols(), re, 1, (double *)result[0], 2);
cblas_dcopy(re.nrows()*re.ncols(), im, 1, (double *)result[0]+1, 2);
return result;
}
NRMat<double> matrixfunction(NRMat<double> a, complex<double> NRMat<double> matrixfunction(NRMat<double> a, complex<double>
(*f)(const complex<double> &), const bool adjust) (*f)(const complex<double> &), const bool adjust)
@ -735,13 +753,13 @@ NRMat<complex<double> > a0=complexify(a);
gdiagonalize(a, w, &u, &v);//a gets destroyed, eigenvectors are rows gdiagonalize(a, w, &u, &v);//a gets destroyed, eigenvectors are rows
NRVec< complex<double> > z = diagofproduct(u, v, 1, 1); NRVec< complex<double> > z = diagofproduct(u, v, 1, 1);
/* /*
cout <<"TEST matrixfunction\n"<<w<<u<<v<<z; std::cout <<"TEST matrixfunction\n"<<w<<u<<v<<z;
cout <<"TEST matrixfunction1 "<< u*a0 - diagonalmatrix(w)*u<<endl; std::cout <<"TEST matrixfunction1 "<< u*a0 - diagonalmatrix(w)*u<<std::endl;
cout <<"TEST matrixfunction2 "<< a0*v.transpose(1) - v.transpose(1)*diagonalmatrix(w)<<endl; std::cout <<"TEST matrixfunction2 "<< a0*v.transpose(1) - v.transpose(1)*diagonalmatrix(w)<<std::endl;
cout <<"TEST matrixfunction3 "<< u*v.transpose(1)<<diagonalmatrix(z)<<endl; std::cout <<"TEST matrixfunction3 "<< u*v.transpose(1)<<diagonalmatrix(z)<<std::endl;
NRVec< complex<double> > wz(n); NRVec< complex<double> > wz(n);
for (int i=0; i<a.nrows(); i++) wz[i] = w[i]/z[i]; for (int i=0; i<a.nrows(); i++) wz[i] = w[i]/z[i];
cout <<"TEST matrixfunction4 "<< a0<< v.transpose(true)*diagonalmatrix(wz)*u<<endl; std::cout <<"TEST matrixfunction4 "<< a0<< v.transpose(true)*diagonalmatrix(wz)*u<<std::endl;
*/ */
for (int i=0; i<a.nrows(); i++) w[i] = (*f)(w[i])/z[i]; for (int i=0; i<a.nrows(); i++) w[i] = (*f)(w[i])/z[i];
@ -751,7 +769,7 @@ cout <<"TEST matrixfunction4 "<< a0<< v.transpose(true)*diagonalmatrix(wz)*u<<en
r.gemm(0.0, v, 'c', u, 'n', 1.0); r.gemm(0.0, v, 'c', u, 'n', 1.0);
double inorm = cblas_dnrm2(n*n, (double *)r[0]+1, 2); double inorm = cblas_dnrm2(n*n, (double *)r[0]+1, 2);
if (inorm > 1e-10) { if (inorm > 1e-10) {
cout << "norm = " << inorm << endl; std::cout << "norm = " << inorm << std::endl;
laerror("nonzero norm of imaginary part of real matrixfunction"); laerror("nonzero norm of imaginary part of real matrixfunction");
} }
return realpart(r); return realpart(r);
@ -817,18 +835,18 @@ complex<double> myclog (const complex<double> &x)
complex<double> mycexp (const complex<double> &x) complex<double> mycexp (const complex<double> &x)
{ {
return exp(x); return std::exp(x);
} }
complex<double> sqrtinv (const complex<double> &x) complex<double> sqrtinv (const complex<double> &x)
{ {
return 1./sqrt(x); return 1./std::sqrt(x);
} }
double sqrtinv (const double x) double sqrtinv (const double x)
{ {
return 1./sqrt(x); return 1./std::sqrt(x);
} }
@ -959,7 +977,7 @@ for(j=i; j<n; ++j)
if(i==j) if(i==j)
{ {
if(x<=0) laerror("not positive definite metric in gendiagonalize"); if(x<=0) laerror("not positive definite metric in gendiagonalize");
dl[i] = sqrt(x); dl[i] = std::sqrt(x);
} }
else else
b(j,i) = x / dl[i]; b(j,i) = x / dl[i];
@ -1023,12 +1041,13 @@ if(nchange&1)//still adjust to get determinant=1
{ {
int imin=-1; double min=1e200; int imin=-1; double min=1e200;
for(int i=0; i<n;++i) for(int i=0; i<n;++i)
if(abs(v[i][i])<min) if(std::abs(v[i][i])<min)
{ {
imin=i; imin=i;
min=abs(v[i][i]); min=std::abs(v[i][i]);
} }
cblas_dscal(n,-1.,v[imin],1); cblas_dscal(n,-1.,v[imin],1);
} }
} }
}//namespace

21
nonclass.h
View File

@ -17,12 +17,12 @@
*/ */
#ifndef _LA_NONCLASS_H_ #ifndef _LA_NONCLASS_H_
#define _LA_NONCLASS_H_ #define _LA_NONCLASS_H_
#include "vec.h" #include "vec.h"
#include "smat.h" #include "smat.h"
#include "mat.h" #include "mat.h"
#include "la_traits.h" #include "la_traits.h"
namespace LA {
//MISC //MISC
template <class T> template <class T>
@ -120,11 +120,11 @@ extern complex<double> sqrtinv(const complex<double> &);
extern double sqrtinv(const double); extern double sqrtinv(const double);
//functions on matrices //functions on matrices
inline NRMat<double> sqrt(const NRSMat<double> &a) { return matrixfunction(a,&sqrt); } inline NRMat<double> sqrt(const NRSMat<double> &a) { return matrixfunction(a,&std::sqrt); }
inline NRMat<double> sqrtinv(const NRSMat<double> &a) { return matrixfunction(a,&sqrtinv); } inline NRMat<double> sqrtinv(const NRSMat<double> &a) { return matrixfunction(a,&sqrtinv); }
inline NRMat<double> realsqrt(const NRMat<double> &a) { return realmatrixfunction(a,&sqrt); } inline NRMat<double> realsqrt(const NRMat<double> &a) { return realmatrixfunction(a,&std::sqrt); }
inline NRMat<double> realsqrtinv(const NRMat<double> &a) { return realmatrixfunction(a,&sqrtinv); } inline NRMat<double> realsqrtinv(const NRMat<double> &a) { return realmatrixfunction(a,&sqrtinv); }
inline NRMat<double> log(const NRSMat<double> &a) { return matrixfunction(a,&log); } inline NRMat<double> log(const NRSMat<double> &a) { return matrixfunction(a,&std::log); }
extern NRMat<double> log(const NRMat<double> &a); extern NRMat<double> log(const NRMat<double> &a);
extern NRMat<double> exp0(const NRMat<double> &a); extern NRMat<double> exp0(const NRMat<double> &a);
@ -133,6 +133,7 @@ extern const NRMat<double> realpart(const NRMat< complex<double> >&);
extern const NRMat<double> imagpart(const NRMat< complex<double> >&); extern const NRMat<double> imagpart(const NRMat< complex<double> >&);
extern const NRMat< complex<double> > realmatrix (const NRMat<double>&); extern const NRMat< complex<double> > realmatrix (const NRMat<double>&);
extern const NRMat< complex<double> > imagmatrix (const NRMat<double>&); extern const NRMat< complex<double> > imagmatrix (const NRMat<double>&);
extern const NRMat< complex<double> > complexmatrix (const NRMat<double>&, const NRMat<double>&);
//inverse by means of linear solve, preserving rhs intact //inverse by means of linear solve, preserving rhs intact
template<typename T> template<typename T>
@ -186,7 +187,7 @@ if(ignoresign)
{ {
for(int i=0; i<nrows; ++i) for(int i=0; i<nrows; ++i)
for(int j=0; j<ncols; ++j) for(int j=0; j<ncols; ++j)
r(i,j) = rows[i]*cols[j]<0?0.:a(abs(rows[i])+indexshift,abs(cols[j])+indexshift); r(i,j) = rows[i]*cols[j]<0?0.:a(std::abs(rows[i])+indexshift,std::abs(cols[j])+indexshift);
} }
else else
{ {
@ -201,7 +202,7 @@ if(ignoresign)
{ {
for(int i=0; i<nrows; ++i) for(int i=0; i<nrows; ++i)
for(int j=0; j<ncols; ++j) for(int j=0; j<ncols; ++j)
r(i,j) = a(abs(rows[i])+indexshift,abs(cols[j])+indexshift); r(i,j) = a(std::abs(rows[i])+indexshift,std::abs(cols[j])+indexshift);
} }
else else
{ {
@ -218,11 +219,6 @@ return r;
//auxiliary routine to adjust eigenvectors to guarantee real logarithm //auxiliary routine to adjust eigenvectors to guarantee real logarithm
extern void adjustphases(NRMat<double> &v); extern void adjustphases(NRMat<double> &v);
template<class T>
T abssqr(const complex<T> &x)
{
return x.real()*x.real()+x.imag()*x.imag();
}
//declaration of template interface to cblas routines with full options available //declaration of template interface to cblas routines with full options available
//just to facilitate easy change between float, double, complex in a user code //just to facilitate easy change between float, double, complex in a user code
@ -272,6 +268,5 @@ return r;
} }
}//namespace
#endif #endif

5
permutation.h
View File

@ -15,6 +15,9 @@
You should have received a copy of the GNU General Public License You should have received a copy of the GNU General Public License
along with this program. If not, see <http://www.gnu.org/licenses/>. along with this program. If not, see <http://www.gnu.org/licenses/>.
*/ */
#ifndef _PERMUTATION_H
#define _PERMUTATION_H
namespace LA {
template<typename T> template<typename T>
const NRVec<T> inversepermutation(const NRVec<T> &p, const T offset=0) const NRVec<T> inversepermutation(const NRVec<T> &p, const T offset=0)
{ {
@ -24,3 +27,5 @@ if(!offset) for(int i=0; i<n; ++i) q[p[i]]=i;
else for(int i=0; i<n; ++i) q[p[i]-offset]=i+offset; else for(int i=0; i<n; ++i) q[p[i]-offset]=i+offset;
return q; return q;
} }
}//namespace
#endif

3
qsort.h
View File

@ -18,6 +18,8 @@
#ifndef _QSORT_H #ifndef _QSORT_H
#define _QSORT_H #define _QSORT_H
//quicksort, returns parity of the permutation //quicksort, returns parity of the permutation
//
namespace LA {
template<typename INDEX, typename COMPAR> template<typename INDEX, typename COMPAR>
int genqsort(INDEX l, INDEX r,COMPAR (*cmp)(const INDEX, const INDEX), void (*swap)(const INDEX,const INDEX)) int genqsort(INDEX l, INDEX r,COMPAR (*cmp)(const INDEX, const INDEX), void (*swap)(const INDEX,const INDEX))
@ -181,4 +183,5 @@ else
return parity; return parity;
} }
}//namespace
#endif #endif

21
smat.cc
View File

@ -31,7 +31,7 @@ extern ssize_t write(int, const void *, size_t);
// TODO // TODO
// specialize unary minus // specialize unary minus
namespace LA {
/* /*
@ -298,7 +298,7 @@ NRSMat< complex<double> >::dot(const NRSMat< complex<double> > &rhs) const
if (nn != rhs.nn) laerror("dot of incompatible SMat's"); if (nn != rhs.nn) laerror("dot of incompatible SMat's");
#endif #endif
complex<double> dot; complex<double> dot;
cblas_zdotc_sub(NN2, (void *)v, 1, (void *)rhs.v, 1, (void *)(&dot)); cblas_zdotc_sub(NN2, v, 1, rhs.v, 1, &dot);
return dot; return dot;
} }
@ -322,7 +322,7 @@ NRSMat< complex<double> >::dot(const NRVec< complex<double> > &rhs) const
if (NN2 != rhs.nn) laerror("dot of incompatible SMat's"); if (NN2 != rhs.nn) laerror("dot of incompatible SMat's");
#endif #endif
complex<double> dot; complex<double> dot;
cblas_zdotc_sub(NN2, (void *)v, 1, (void *)rhs.v, 1, (void *)(&dot)); cblas_zdotc_sub(NN2, v, 1, rhs.v, 1, &dot);
return dot; return dot;
} }
@ -341,7 +341,7 @@ const double NRSMat<double>::norm(const double scalar) const
if (i == j) tmp -= scalar; if (i == j) tmp -= scalar;
sum += tmp*tmp; sum += tmp*tmp;
} }
return sqrt(sum); return std::sqrt(sum);
} }
@ -350,7 +350,7 @@ template<>
const double NRSMat< complex<double> >::norm(const complex<double> scalar) const const double NRSMat< complex<double> >::norm(const complex<double> scalar) const
{ {
if (!(scalar.real()) && !(scalar.imag())) if (!(scalar.real()) && !(scalar.imag()))
return cblas_dznrm2(NN2, (void *)v, 1); return cblas_dznrm2(NN2, v, 1);
double sum = 0; double sum = 0;
complex<double> tmp; complex<double> tmp;
int k = 0; int k = 0;
@ -360,7 +360,7 @@ const double NRSMat< complex<double> >::norm(const complex<double> scalar) const
if (i == j) tmp -= scalar; if (i == j) tmp -= scalar;
sum += tmp.real()*tmp.real() + tmp.imag()*tmp.imag(); sum += tmp.real()*tmp.real() + tmp.imag()*tmp.imag();
} }
return sqrt(sum); return std::sqrt(sum);
} }
@ -388,7 +388,7 @@ void NRSMat< complex<double> >::axpy(const complex<double> alpha,
if (nn != x.nn) laerror("axpy of incompatible SMats"); if (nn != x.nn) laerror("axpy of incompatible SMats");
#endif #endif
copyonwrite(); copyonwrite();
cblas_zaxpy(nn, (void *)(&alpha), (void *)x.v, 1, (void *)v, 1); cblas_zaxpy(nn, &alpha, x.v, 1, v, 1);
} }
//complex from real //complex from real
@ -407,9 +407,16 @@ cblas_dcopy(nn*(nn+1)/2,&rhs(0,0),1,((double *)v) + (imagpart?1:0),2);
////// forced instantization in the corresponding object file ////// forced instantization in the corresponding object file
template class NRSMat<double>; template class NRSMat<double>;
template class NRSMat< complex<double> >; template class NRSMat< complex<double> >;
template class NRSMat<long long>;
template class NRSMat<long>;
template class NRSMat<int>; template class NRSMat<int>;
template class NRSMat<short>; template class NRSMat<short>;
template class NRSMat<char>; template class NRSMat<char>;
template class NRSMat<unsigned char>;
template class NRSMat<unsigned short>;
template class NRSMat<unsigned int>; template class NRSMat<unsigned int>;
template class NRSMat<unsigned long>; template class NRSMat<unsigned long>;
template class NRSMat<unsigned long long>;
}//namespace

20
smat.h
View File

@ -19,9 +19,9 @@
*/ */
#ifndef _LA_SMAT_H_ #ifndef _LA_SMAT_H_
#define _LA_SMAT_H_ #define _LA_SMAT_H_
#include "la_traits.h" #include "la_traits.h"
namespace LA {
#define NN2 (nn*(nn+1)/2) #define NN2 (nn*(nn+1)/2)
template <class T> template <class T>
class NRSMat { // symmetric or complex hermitean matrix in packed form class NRSMat { // symmetric or complex hermitean matrix in packed form
@ -94,15 +94,17 @@ public:
void fscanf(FILE *f, const char *format); void fscanf(FILE *f, const char *format);
//members concerning sparse matrix //members concerning sparse matrix
explicit NRSMat(const SparseMat<T> &rhs); // dense from sparse explicit NRSMat(const SparseMat<T> &rhs); // dense from sparse
explicit NRSMat(const SparseSMat<T> &rhs); // dense from sparse
inline void simplify() {}; //just for compatibility with sparse ones inline void simplify() {}; //just for compatibility with sparse ones
bool issymmetric() const {return 1;} bool issymmetric() const {return 1;}
}; };
}//namespace
//due to mutual includes this has to be after full class declaration //due to mutual includes this has to be after full class declaration
#include "vec.h" #include "vec.h"
#include "mat.h" #include "mat.h"
#include "sparsemat.h"
namespace LA {
// ctors // ctors
template <typename T> template <typename T>
@ -162,7 +164,7 @@ inline NRSMat< complex<double> > &
NRSMat< complex<double> >::operator*=(const complex<double> & a) NRSMat< complex<double> >::operator*=(const complex<double> & a)
{ {
copyonwrite(); copyonwrite();
cblas_zscal(NN2, (void *)(&a), (void *)v, 1); cblas_zscal(NN2, &a, v, 1);
return *this; return *this;
} }
template <typename T> template <typename T>
@ -214,7 +216,7 @@ NRSMat< complex<double> >::operator+=(const NRSMat< complex<double> > & rhs)
if (nn != rhs.nn) laerror("incompatible SMats in SMat::operator+="); if (nn != rhs.nn) laerror("incompatible SMats in SMat::operator+=");
#endif #endif
copyonwrite(); copyonwrite();
cblas_zaxpy(NN2, (void *)(&CONE), (void *)(&rhs.v), 1, (void *)(&v), 1); cblas_zaxpy(NN2, &CONE, rhs.v, 1, v, 1);
return *this; return *this;
} }
@ -251,7 +253,7 @@ NRSMat< complex<double> >::operator-=(const NRSMat< complex<double> > & rhs)
if (nn != rhs.nn) laerror("incompatible SMats in SMat::operator-="); if (nn != rhs.nn) laerror("incompatible SMats in SMat::operator-=");
#endif #endif
copyonwrite(); copyonwrite();
cblas_zaxpy(NN2, (void *)(&CMONE), (void *)(&rhs.v), 1, (void *)(&v), 1); cblas_zaxpy(NN2, &CMONE, rhs.v, 1, v, 1);
return *this; return *this;
} }
@ -406,7 +408,7 @@ inline const double NRSMat<double>::amax() const
template<> template<>
inline const complex<double> NRSMat< complex<double> >::amax() const inline const complex<double> NRSMat< complex<double> >::amax() const
{ {
return v[cblas_izamax(NN2, (void *)v, 1)]; return v[cblas_izamax(NN2, v, 1)];
} }
// reference pointer to Smat // reference pointer to Smat
@ -554,7 +556,7 @@ return r;
// I/O // I/O
template <typename T> template <typename T>
ostream& operator<<(ostream &s, const NRSMat<T> &x) std::ostream& operator<<(std::ostream &s, const NRSMat<T> &x)
{ {
int i,j,n; int i,j,n;
n=x.nrows(); n=x.nrows();
@ -568,7 +570,7 @@ ostream& operator<<(ostream &s, const NRSMat<T> &x)
template <typename T> template <typename T>
istream& operator>>(istream &s, NRSMat<T> &x) std::istream& operator>>(std::istream &s, NRSMat<T> &x)
{ {
int i,j,n,m; int i,j,n,m;
s >> n >> m; s >> n >> m;
@ -618,5 +620,5 @@ public:
} }
}; };
}//namespace
#endif /* _LA_SMAT_H_ */ #endif /* _LA_SMAT_H_ */

124
sparsemat.cc
View File

@ -25,6 +25,8 @@
#include <errno.h> #include <errno.h>
#include "sparsemat.h" #include "sparsemat.h"
namespace LA {
template<typename T> template<typename T>
static inline const T MAX(const T &a, const T &b) static inline const T MAX(const T &a, const T &b)
{return b > a ? (b) : (a);} {return b > a ? (b) : (a);}
@ -333,11 +335,7 @@ for(i=1; i<n;i++)
//check if summed to zero //check if summed to zero
for(i=0; i<n;i++) if(rowsorted[i] && for(i=0; i<n;i++) if(rowsorted[i] &&
#ifdef SPARSEEPSILON std::abs(rowsorted[i]->elem)<=SPARSEEPSILON
abs(rowsorted[i]->elem)<SPARSEEPSILON
#else
! rowsorted[i]->elem
#endif
) {delete rowsorted[i]; rowsorted[i]=NULL;} ) {delete rowsorted[i]; rowsorted[i]=NULL;}
//restore connectivity //restore connectivity
@ -397,11 +395,7 @@ void SparseMat<T>::addsafe(const SPMatindex n, const SPMatindex m, const T elem)
#ifdef debug #ifdef debug
if(n<0||n>=nn||m<0||m>=mm) laerror("SparseMat out of range"); if(n<0||n>=nn||m<0||m>=mm) laerror("SparseMat out of range");
#endif #endif
#ifdef SPARSEEPSILON if(std::abs(elem)<=SPARSEEPSILON) return;
if(abs(elem)<SPARSEEPSILON) return;
#else
if(!elem) return;
#endif
if(!count) {count=new int; *count=1; list=NULL;} //blank new matrix if(!count) {count=new int; *count=1; list=NULL;} //blank new matrix
else copyonwrite(); //makes also unsort else copyonwrite(); //makes also unsort
add(n,m,elem); add(n,m,elem);
@ -452,9 +446,7 @@ register matel<T> *l=rhs.list;
while(l) while(l)
{ {
if(rhs.symmetric || lower && l->row <=l->col || !lower && l->row >=l->col) if(rhs.symmetric || lower && l->row <=l->col || !lower && l->row >=l->col)
#ifdef SPARSEEPSILON if(std::abs(l->elem)>SPARSEEPSILON)
if(abs(l->elem)>SPARSEEPSILON)
#endif
add( l->row,l->col,sign=='+'?l->elem:- l->elem); add( l->row,l->col,sign=='+'?l->elem:- l->elem);
l=l->next; l=l->next;
} }
@ -473,21 +465,18 @@ register matel<T> *l=rhs.list;
if(symmetrize) if(symmetrize)
while(l) while(l)
{ {
#ifdef SPARSEEPSILON if(std::abs(l->elem)>SPARSEEPSILON)
if(abs(l->elem)>SPARSEEPSILON)
#endif
{add( l->row,l->col,l->elem); if( l->row!=l->col) add( l->col,l->row,l->elem);} {add( l->row,l->col,l->elem); if( l->row!=l->col) add( l->col,l->row,l->elem);}
l=l->next; l=l->next;
} }
else else
while(l) while(l)
{ {
#ifdef SPARSEEPSILON if(std::abs(l->elem)>SPARSEEPSILON)
if(abs(l->elem)>SPARSEEPSILON)
#endif
add( l->row,l->col,l->elem); add( l->row,l->col,l->elem);
l=l->next; l=l->next;
} }
simplify(); //maybe leave up to the user
return *this; return *this;
} }
@ -503,21 +492,18 @@ register matel<T> *l=rhs.list;
if(symmetrize) if(symmetrize)
while(l) while(l)
{ {
#ifdef SPARSEEPSILON if(std::abs(l->elem)>SPARSEEPSILON)
if(abs(l->elem)>SPARSEEPSILON)
#endif
{add( l->row,l->col,- l->elem); if( l->row!=l->col) add( l->col,l->row,- l->elem);} {add( l->row,l->col,- l->elem); if( l->row!=l->col) add( l->col,l->row,- l->elem);}
l=l->next; l=l->next;
} }
else else
while(l) while(l)
{ {
#ifdef SPARSEEPSILON if(std::abs(l->elem)>SPARSEEPSILON)
if(abs(l->elem)>SPARSEEPSILON)
#endif
add( l->row,l->col,- l->elem); add( l->row,l->col,- l->elem);
l=l->next; l=l->next;
} }
simplify(); //maybe leave up to the user
return *this; return *this;
} }
@ -537,11 +523,7 @@ SPMatindex i,j;
for(i=0;i<nn;++i) for(i=0;i<nn;++i)
for(j=0; j<mm;++j) for(j=0; j<mm;++j)
{register T t(rhs(i,j)); {register T t(rhs(i,j));
#ifdef SPARSEEPSILON if( std::abs(t)>SPARSEEPSILON)
if( abs(t)>SPARSEEPSILON)
#else
if(t)
#endif
add(i,j,t); add(i,j,t);
} }
} }
@ -660,40 +642,38 @@ else while(l)
//assignment of a scalar matrix //assignment of a scalar matrix
template <class T> template <class T>
SparseMat<T> & SparseMat<T>::operator=(const T a) SparseMat<T> & SparseMat<T>::operator=(const T &a)
{ {
if(!count ||nn<=0||mm<=0) laerror("assignment of scalar to undefined sparse matrix"); if(!count ||nn<=0||mm<=0) laerror("assignment of scalar to undefined sparse matrix");
if(nn!=mm) laerror("assignment of scalar to non-square sparse matrix"); if(nn!=mm) laerror("assignment of scalar to non-square sparse matrix");
resize(nn,mm);//clear resize(nn,mm);//clear
#ifdef SPARSEEPSILON if(std::abs(a)<=SPARSEEPSILON) return *this;
if(abs(a)<SPARSEEPSILON) return *this;
#else
if(a==(T)0) return *this;
#endif
SPMatindex i; SPMatindex i;
for(i=0;i<nn;++i) add(i,i,a); for(i=0;i<nn;++i) add(i,i,a);
return *this; return *this;
} }
template <class T> template <class T>
SparseMat<T> & SparseMat<T>::operator+=(const T a) SparseMat<T> & SparseMat<T>::operator+=(const T &a)
{ {
if(!count ||nn<=0||mm<=0) laerror("assignment of scalar to undefined sparse matrix"); if(!count ||nn<=0||mm<=0) laerror("assignment of scalar to undefined sparse matrix");
if(nn!=mm) laerror("assignment of scalar to non-square sparse matrix"); if(nn!=mm) laerror("assignment of scalar to non-square sparse matrix");
if(a==(T)0) return *this; if(a==(T)0) return *this;
SPMatindex i; SPMatindex i;
for(i=0;i<nn;++i) add(i,i,a); for(i=0;i<nn;++i) add(i,i,a);
simplify(); //maybe leave up to the user
return *this; return *this;
} }
template <class T> template <class T>
SparseMat<T> & SparseMat<T>::operator-=(const T a) SparseMat<T> & SparseMat<T>::operator-=(const T &a)
{ {
if(!count ||nn<=0||mm<=0) laerror("assignment of scalar to undefined sparse matrix"); if(!count ||nn<=0||mm<=0) laerror("assignment of scalar to undefined sparse matrix");
if(nn!=mm) laerror("assignment of scalar to non-square sparse matrix"); if(nn!=mm) laerror("assignment of scalar to non-square sparse matrix");
if(a==(T)0) return *this; if(a==(T)0) return *this;
SPMatindex i; SPMatindex i;
for(i=0;i<nn;++i) add(i,i,-a); for(i=0;i<nn;++i) add(i,i,-a);
simplify(); //maybe leave up to the user
return *this; return *this;
} }
@ -715,11 +695,7 @@ for(i=0;i<nn;++i)
for(j=0; j<=i;++j) for(j=0; j<=i;++j)
{register T t; {register T t;
if( if(
#ifdef SPARSEEPSILON std::abs(t=rhs(i,j))>SPARSEEPSILON
abs(t=rhs(i,j))>SPARSEEPSILON
#else
t=rhs(i,j)
#endif
) add(i,j,t); ) add(i,j,t);
} }
} }
@ -751,9 +727,7 @@ matel<T> *l=list;
while(l) //include the mirror picture of elements into the list while(l) //include the mirror picture of elements into the list
{ {
if( if(
#ifdef SPARSEEPSILON std::abs(l->elem)>SPARSEEPSILON &&
abs(l->elem)>SPARSEEPSILON &&
#endif
l->row!=l->col) add(l->col,l->row,l->elem); //not OK for complex-hermitean l->row!=l->col) add(l->col,l->row,l->elem); //not OK for complex-hermitean
l=l->next; l=l->next;
} }
@ -761,12 +735,12 @@ while(l) //include the mirror picture of elements into the list
template <class T> template <class T>
SparseMat<T> & SparseMat<T>::operator*=(const T a) SparseMat<T> & SparseMat<T>::operator*=(const T &a)
{ {
if(!count) laerror("operator*= on undefined lhs"); if(!count) laerror("operator*= on undefined lhs");
if(!list||a==(T)1) return *this; if(!list||a==(T)1) return *this;
if(a==(T)0) resize(nn,mm); if(a==(T)0) {clear(); return *this;}
else copyonwrite(); copyonwrite();
register matel<T> *l=list; register matel<T> *l=list;
while(l) while(l)
@ -936,15 +910,14 @@ return sum;
} }
//not OK for complex hermitean
template<class T> template<class T>
const T SparseMat<T>::norm(const T scalar) const const typename LA_traits<T>::normtype SparseMat<T>::norm(const T scalar) const
{ {
if(!list) return T(0); if(!list) return typename LA_traits<T>::normtype(0);
const_cast<SparseMat<T> *>(this)->simplify(); const_cast<SparseMat<T> *>(this)->simplify();
matel<T> *l=list; matel<T> *l=list;
T sum(0); typename LA_traits<T>::normtype sum(0);
if(scalar!=(T)0) if(scalar!=(T)0)
{ {
if(symmetric) if(symmetric)
@ -953,7 +926,7 @@ if(scalar!=(T)0)
T hlp=l->elem; T hlp=l->elem;
bool b=l->row==l->col; bool b=l->row==l->col;
if(b) hlp-=scalar; if(b) hlp-=scalar;
T tmp=hlp*hlp; typename LA_traits<T>::normtype tmp=LA_traits<T>::sqrabs(hlp);
sum+= tmp; sum+= tmp;
if(!b) sum+=tmp; if(!b) sum+=tmp;
l=l->next; l=l->next;
@ -963,7 +936,7 @@ if(scalar!=(T)0)
{ {
T hlp=l->elem; T hlp=l->elem;
if(l->row==l->col) hlp-=scalar; if(l->row==l->col) hlp-=scalar;
sum+= hlp*hlp; sum+= LA_traits<T>::sqrabs(hlp);
l=l->next; l=l->next;
} }
} }
@ -972,7 +945,7 @@ else
if(symmetric) if(symmetric)
while(l) while(l)
{ {
T tmp=l->elem*l->elem; typename LA_traits<T>::normtype tmp=LA_traits<T>::sqrabs(l->elem);
sum+= tmp; sum+= tmp;
if(l->row!=l->col) sum+=tmp; if(l->row!=l->col) sum+=tmp;
l=l->next; l=l->next;
@ -980,11 +953,11 @@ else
else else
while(l) while(l)
{ {
sum+= l->elem*l->elem; sum+= LA_traits<T>::sqrabs(l->elem);
l=l->next; l=l->next;
} }
} }
return sqrt(sum); //not OK for int, would need traits technique return (typename LA_traits<T>::normtype) std::sqrt(sum); //not OK for int, would need traits technique
} }
@ -1009,9 +982,7 @@ if(transp)
while(l) while(l)
{ {
register T t=alpha*l->elem; register T t=alpha*l->elem;
#ifdef SPARSEEPSILON if(std::abs(t)>SPARSEEPSILON)
if(abs(t)>SPARSEEPSILON)
#endif
add( l->col,l->row,t); add( l->col,l->row,t);
l=l->next; l=l->next;
} }
@ -1019,9 +990,7 @@ else
while(l) while(l)
{ {
register T t=alpha*l->elem; register T t=alpha*l->elem;
#ifdef SPARSEEPSILON if(std::abs(t)>SPARSEEPSILON)
if(abs(t)>SPARSEEPSILON)
#endif
add( l->row,l->col,t); add( l->row,l->col,t);
l=l->next; l=l->next;
} }
@ -1227,11 +1196,7 @@ incsize(rhs.nrows(), rhs.ncols());
T tmp; T tmp;
for(SPMatindex i=0; i<(SPMatindex)rhs.nrows(); ++i) for(SPMatindex i=0; i<(SPMatindex)rhs.nrows(); ++i)
for(SPMatindex j=0; j<(SPMatindex)rhs.ncols(); ++j) for(SPMatindex j=0; j<(SPMatindex)rhs.ncols(); ++j)
#ifdef SPARSEEPSILON if(std::abs(tmp=rhs(i,j))>SPARSEEPSILON) add(n0+i,m0+j,tmp);
if(abs(tmp=rhs(i,j))>SPARSEEPSILON) add(n0+i,m0+j,tmp);
#else
if((tmp=rhs(i,j))!=(T)0) add(n0+i,m0+j,tmp);
#endif
return *this; return *this;
} }
@ -1247,11 +1212,7 @@ T tmp;
incsize(rhs.nrows(), rhs.ncols()); incsize(rhs.nrows(), rhs.ncols());
for(SPMatindex i=0; i<(SPMatindex)rhs.nrows(); ++i) for(SPMatindex i=0; i<(SPMatindex)rhs.nrows(); ++i)
for(SPMatindex j=0; j<(SPMatindex)rhs.ncols(); ++j) for(SPMatindex j=0; j<(SPMatindex)rhs.ncols(); ++j)
#ifdef SPARSEEPSILON if(std::abs(tmp=rhs(i,j))>SPARSEEPSILON) add(n0+i,m0+j,tmp);
if(abs(tmp=rhs(i,j))>SPARSEEPSILON) add(n0+i,m0+j,tmp);
#else
if((tmp=rhs(i,j))!=(T)0) add(n0+i,m0+j,tmp);
#endif
return *this; return *this;
} }
@ -1270,9 +1231,7 @@ incsize(rhs.nrows(), rhs.ncols());
register matel<T> *l=rhs.list; register matel<T> *l=rhs.list;
while(l) while(l)
{ {
#ifdef SPARSEEPSILON if(std::abs(l->elem)>SPARSEEPSILON)
if(abs(l->elem)>SPARSEEPSILON)
#endif
add(n0+l->row,m0+l->col,l->elem); add(n0+l->row,m0+l->col,l->elem);
l=l->next; l=l->next;
} }
@ -1303,18 +1262,18 @@ template SparseMat<T>::SparseMat(const NRMat<T> &rhs); \
template SparseMat<T>::SparseMat(const NRSMat<T> &rhs); \ template SparseMat<T>::SparseMat(const NRSMat<T> &rhs); \
template void SparseMat<T>::transposeme(); \ template void SparseMat<T>::transposeme(); \
template const T* SparseMat<T>::diagonalof(NRVec<T> &r, const bool divide, bool cache) const; \ template const T* SparseMat<T>::diagonalof(NRVec<T> &r, const bool divide, bool cache) const; \
template SparseMat<T> & SparseMat<T>::operator*=(const T a); \ template SparseMat<T> & SparseMat<T>::operator*=(const T &a); \
template void SparseMat<T>::setunsymmetric(); \ template void SparseMat<T>::setunsymmetric(); \
template SparseMat<T> & SparseMat<T>::operator=(const T a); \ template SparseMat<T> & SparseMat<T>::operator=(const T &a); \
template SparseMat<T> & SparseMat<T>::operator+=(const T a); \ template SparseMat<T> & SparseMat<T>::operator+=(const T &a); \
template SparseMat<T> & SparseMat<T>::operator-=(const T a); \ template SparseMat<T> & SparseMat<T>::operator-=(const T &a); \
template NRMat<T>::NRMat(const SparseMat<T> &rhs); \ template NRMat<T>::NRMat(const SparseMat<T> &rhs); \
template NRSMat<T>::NRSMat(const SparseMat<T> &rhs); \ template NRSMat<T>::NRSMat(const SparseMat<T> &rhs); \
template NRVec<T>::NRVec(const SparseMat<T> &rhs); \ template NRVec<T>::NRVec(const SparseMat<T> &rhs); \
template const NRVec<T> NRVec<T>::operator*(const SparseMat<T> &mat) const; \ template const NRVec<T> NRVec<T>::operator*(const SparseMat<T> &mat) const; \
template SparseMat<T> & SparseMat<T>::join(SparseMat<T> &rhs); \ template SparseMat<T> & SparseMat<T>::join(SparseMat<T> &rhs); \
template const T SparseMat<T>::trace() const; \ template const T SparseMat<T>::trace() const; \
template const T SparseMat<T>::norm(const T scalar) const; \ template const LA_traits<T>::normtype SparseMat<T>::norm(const T scalar) const; \
template void SparseMat<T>::axpy(const T alpha, const SparseMat<T> &x, const bool transp); \ template void SparseMat<T>::axpy(const T alpha, const SparseMat<T> &x, const bool transp); \
template const SparseMat<T> SparseMat<T>::operator*(const SparseMat<T> &rhs) const; \ template const SparseMat<T> SparseMat<T>::operator*(const SparseMat<T> &rhs) const; \
template const T SparseMat<T>::dot(const SparseMat<T> &rhs) const; \ template const T SparseMat<T>::dot(const SparseMat<T> &rhs) const; \
@ -1336,3 +1295,4 @@ INSTANTIZE(complex<double>) //some functions are not OK for hermitean matrices,
template class SparseMat<double>; template class SparseMat<double>;
template class SparseMat<complex<double> >; template class SparseMat<complex<double> >;
}//namespace

23
sparsemat.h
View File

@ -17,13 +17,13 @@
*/ */
#ifndef _SPARSEMAT_H_ #ifndef _SPARSEMAT_H_
#define _SPARSEMAT_H_ #define _SPARSEMAT_H_
#include "la_traits.h" #include "la_traits.h"
namespace LA {
//threshold for neglecting elements, if not defined, no tests are done except exact zero test in simplify - might be even faster //threshold for neglecting elements, if not defined, no tests are done except exact zero test in simplify - might be even faster
//seems to perform better with a threshold, in spite of abs() tests //seems to perform better with a threshold, in spite of abs() tests
#define SPARSEEPSILON 1e-13 const double SPARSEEPSILON=1e-35;
typedef unsigned int SPMatindex; typedef unsigned int SPMatindex;
typedef int SPMatindexdiff; //more clear would be to use traits typedef int SPMatindexdiff; //more clear would be to use traits
@ -81,10 +81,10 @@ public:
explicit SparseMat(const NRMat<T> &rhs); //construct from a dense one explicit SparseMat(const NRMat<T> &rhs); //construct from a dense one
explicit SparseMat(const NRSMat<T> &rhs); //construct from a dense symmetric one explicit SparseMat(const NRSMat<T> &rhs); //construct from a dense symmetric one
SparseMat & operator=(const SparseMat &rhs); SparseMat & operator=(const SparseMat &rhs);
SparseMat & operator=(const T a); //assign a to diagonal SparseMat & operator=(const T &a); //assign a to diagonal
SparseMat & operator+=(const T a); //assign a to diagonal SparseMat & operator+=(const T &a); //assign a to diagonal
SparseMat & operator-=(const T a); //assign a to diagonal SparseMat & operator-=(const T &a); //assign a to diagonal
SparseMat & operator*=(const T a); //multiply by a scalar SparseMat & operator*=(const T &a); //multiply by a scalar
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
@ -135,17 +135,21 @@ public:
void clear() {copyonwrite();unsort();deletelist();} 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 typename LA_traits<T>::normtype norm(const T scalar=(T)0) const; //is const only mathematically, not in internal implementation - we have to simplify first
unsigned int sort(int type) const; unsigned int sort(int type) const;
inline void add(const SPMatindex n, const SPMatindex m, const T elem) {matel<T> *ltmp= new matel<T>; ltmp->next=list; list=ltmp; list->row=n; list->col=m; list->elem=elem;} inline void add(const SPMatindex n, const SPMatindex m, const T elem) {matel<T> *ltmp= new matel<T>; ltmp->next=list; list=ltmp; list->row=n; list->col=m; list->elem=elem;}
void addsafe(const SPMatindex n, const SPMatindex m, const T elem); void addsafe(const SPMatindex n, const SPMatindex m, const T elem);
}; };
}//namespace
//due to mutual includes this has to be after full class declaration //due to mutual includes this has to be after full class declaration
#include "vec.h" #include "vec.h"
#include "smat.h" #include "smat.h"
#include "mat.h" #include "mat.h"
namespace LA {
template <typename T> template <typename T>
inline const NRVec<T> SparseMat<T>::operator*(const NRVec<T> &rhs) const inline const NRVec<T> SparseMat<T>::operator*(const NRVec<T> &rhs) const
{NRVec<T> result(nn); result.gemv((T)0,*this,'n',(T)1,rhs); return result;}; {NRVec<T> result(nn); result.gemv((T)0,*this,'n',(T)1,rhs); return result;};
@ -155,7 +159,7 @@ inline const NRMat<T> SparseMat<T>::operator*(const NRMat<T> &rhs) const
{NRMat<T> result(nn,rhs.ncols()); result.gemm((T)0,*this,'n',rhs,'n',(T)1); return result;}; {NRMat<T> result(nn,rhs.ncols()); result.gemm((T)0,*this,'n',rhs,'n',(T)1); return result;};
template <class T> template <class T>
ostream& operator<<(ostream &s, const SparseMat<T> &x) std::ostream& operator<<(std::ostream &s, const SparseMat<T> &x)
{ {
SPMatindex n,m; SPMatindex n,m;
n=x.nrows(); n=x.nrows();
@ -172,7 +176,7 @@ ostream& operator<<(ostream &s, const SparseMat<T> &x)
} }
template <class T> template <class T>
istream& operator>>(istream &s, SparseMat<T> &x) std::istream& operator>>(std::istream &s, SparseMat<T> &x)
{ {
int i,j; int i,j;
int n,m; int n,m;
@ -298,4 +302,5 @@ while(l)
return *this; return *this;
} }
}//namespace
#endif #endif

13
sparsesmat.cc
View File

@ -25,10 +25,11 @@
#include <errno.h> #include <errno.h>
#include "sparsesmat.h" #include "sparsesmat.h"
namespace LA {
template <typename T> template <typename T>
void SparseSMat<T>::gemm(const T beta, const SparseSMat &a, const char transa, const SparseSMat &b, const char transb, const T alpha) void SparseSMat<T>::gemm(const T beta, const SparseSMat &a, const char transa, const SparseSMat &b, const char transb, const T alpha)
{ {
std::cerr << "enter gemm\n";
(*this) *= beta; (*this) *= beta;
if(alpha==(T)0) return; if(alpha==(T)0) return;
if(a.nn!=b.nn || a.nn!=nn) laerror("incompatible sizes in SparseSMat::gemm"); if(a.nn!=b.nn || a.nn!=nn) laerror("incompatible sizes in SparseSMat::gemm");
@ -64,11 +65,7 @@ for(SPMatindex k=0; k<nn; ++k) //summation loop
add(ai[i],bi[j],prod(i,j),false); add(ai[i],bi[j],prod(i,j),false);
} }
std::cerr << "before simplify in gemm\n"; simplify();
std::cerr << " nterms = "<<
simplify()
<<std::endl; //erase elements below threshold
std::cerr << "regular exit gemm\n";
} }
@ -183,7 +180,7 @@ for(SPMatindex i=0; i<nn; ++i)
} }
else sum += LA_traits<T>::sqrabs(scalar); //missing diagonal element else sum += LA_traits<T>::sqrabs(scalar); //missing diagonal element
return sqrt(sum); return std::sqrt(sum);
} }
@ -208,4 +205,4 @@ INSTANTIZE(complex<double>)
template class SparseSMat<double>; template class SparseSMat<double>;
template class SparseSMat<complex<double> >; template class SparseSMat<complex<double> >;
}//namespace

4
sparsesmat.h
View File

@ -19,7 +19,6 @@
#ifndef _SPARSESMAT_H_ #ifndef _SPARSESMAT_H_
#define _SPARSESMAT_H_ #define _SPARSESMAT_H_
#include <string> #include <string>
#include <cmath> #include <cmath>
#include <stdlib.h> #include <stdlib.h>
@ -36,6 +35,7 @@
#include <map> #include <map>
#include <list> #include <list>
namespace LA {
//symmetric sparse matrix class with a representation for efficient exponentiatiation //symmetric sparse matrix class with a representation for efficient exponentiatiation
//in particular we need a unitary symmetric complex matrix as exp(iH) with H real symmetric //in particular we need a unitary symmetric complex matrix as exp(iH) with H real symmetric
@ -360,5 +360,5 @@ std::istream& operator>>(std::istream &s, SparseSMat<T> &x)
} }
}//namespace
#endif //_SPARSESMAT_H_ #endif //_SPARSESMAT_H_

2
strassen.cc
View File

@ -20,6 +20,7 @@
#include "mat.h" #include "mat.h"
#ifndef NO_STRASSEN #ifndef NO_STRASSEN
namespace LA {
/*Strassen algorithm*/ /*Strassen algorithm*/
// called routine is fortran-compatible // called routine is fortran-compatible
extern "C" void fmm(const char c_transa,const char c_transb,const int m,const int n,const int k,const double alpha, extern "C" void fmm(const char c_transa,const char c_transb,const int m,const int n,const int k,const double alpha,
@ -45,4 +46,5 @@ copyonwrite();
//swap transpositions and order of matrices //swap transpositions and order of matrices
fmm(transb,transa,mm,nn,k,alpha,b,b.mm, a, a.mm, beta,*this, mm,NULL,0); fmm(transb,transa,mm,nn,k,alpha,b,b.mm, a, a.mm, beta,*this, mm,NULL,0);
} }
}//namespace
#endif #endif

66
t.cc
View File

@ -20,18 +20,9 @@
#include <time.h> #include <time.h>
#include "la.h" #include "la.h"
#include "sparsemat.h"
#include "matexp.h"
#include "fourindex.h"
#include "davidson.h"
#include "gmres.h"
#include "conjgrad.h"
#include "diis.h"
#include "bitvector.h"
#ifdef USE_TRACEBACK
#include "traceback.h"
#endif
using namespace std;
using namespace LA;
extern void test(const NRVec<double> &x); extern void test(const NRVec<double> &x);
@ -1335,7 +1326,7 @@ cout << "Unitary\n"<<U;
cout <<"Error = "<<(U*U.transpose() -1.).norm()<<endl; cout <<"Error = "<<(U*U.transpose() -1.).norm()<<endl;
} }
if(1) if(0)
{ {
double t; double t;
NRMat<complex<double> > ham; NRMat<complex<double> > ham;
@ -1397,7 +1388,58 @@ cout <<"sparse sincostimes "<<clock()/((double) (CLOCKS_PER_SEC))-t <<"\n";
NRVec<complex<double> > r4 = cov2 + siv2*complex<double>(0,1); NRVec<complex<double> > r4 = cov2 + siv2*complex<double>(0,1);
cout <<"sincostimes error 2 = "<<(r1-r4).norm()<<endl; cout <<"sincostimes error 2 = "<<(r1-r4).norm()<<endl;
*/ */
}
if(0)
{
double t;
SparseMat<complex<double> > h;
cin >> h;
h *= complex<double>(0,1);
h.simplify();
cout <<"length of hamiltonian "<<h.length()<<endl;
t=clock()/((double) (CLOCKS_PER_SEC));
SparseMat<complex<double> > hexp = exp(h);
cout <<"sparse exp time "<<clock()/((double) (CLOCKS_PER_SEC))-t <<"\n";
cout <<"length of exp hamiltonian "<<hexp.length()<<endl;
NRMat<complex<double> > he(hexp);
NRMat<complex<double> > hec(he);
NRMat<complex<double> > het(he);
hec.conjugateme();
het.transposeme();
//cout <<he;
cout <<"unitarity error "<<(hec*he).norm(1)<<endl;
cout <<"symmetry error "<<(het-he).norm()<<endl;
}
if(0)
{
NRSMat<double> hd(100);
hd.randomize(1);
SparseSMat<double> h(hd);
NRMat<double> rd = hd*hd;
SparseSMat<double> r = h*h;
NRSMat<double>rx(r);
NRMat<double> r2(rx);
cout <<"Error = "<<(r2-rd).norm()<<endl;
}
if(1)
{
SparseSMat<complex<double> > h;
cin>>h;
h *= complex<double>(0,1);
double t=clock()/((double) (CLOCKS_PER_SEC));
SparseSMat<complex<double> > r = exp(h);
cout <<"SparseSMat time "<<clock()/((double) (CLOCKS_PER_SEC))-t <<"\n";
if(h.nrows()<=1024)
{
NRSMat<complex<double> > h3(h);
NRMat<complex<double> > h2(h3);
NRMat<complex<double> >r2 = exp(h2);
cout <<"error = "<<(r2-NRMat<complex<double> >(NRSMat<complex<double> >(r))).norm()<<endl;
cout <<"errorx = "<<(r2-NRSMat<complex<double> >(r)).norm()<<endl;
}
} }

6
t2.cc
View File

@ -28,13 +28,15 @@
#include "traceback.h" #include "traceback.h"
#endif #endif
using namespace LA;
void test(const NRVec<double> &x) void test(const NRVec<double> &x)
{ {
NRMat<double> aa(0.,2,2); NRMat<double> aa(0.,2,2);
NRMat<double> cc(aa); NRMat<double> cc(aa);
cc.copyonwrite(); cc[0][0]=1.; cc.copyonwrite(); cc[0][0]=1.;
cout << aa << cc <<"\n"; std::cout << aa << cc <<"\n";
cout <<"test x= "<<x<<"\n"; std::cout <<"test x= "<<x<<"\n";
} }

2
test.cc
View File

@ -25,6 +25,8 @@
#include "traceback.h" #include "traceback.h"
#endif #endif
using namespace std;
void test2(char *msg) void test2(char *msg)
{ {
laerror(msg); laerror(msg);

76
vec.cc
View File

@ -18,19 +18,22 @@
along with this program. If not, see <http://www.gnu.org/licenses/>. along with this program. If not, see <http://www.gnu.org/licenses/>.
*/ */
#include <iostream> #include <iostream>
#include "vec.h"
#include <stdlib.h> #include <stdlib.h>
#include <sys/types.h> #include <sys/types.h>
#include <sys/stat.h> #include <sys/stat.h>
#include <fcntl.h> #include <fcntl.h>
#include <errno.h> #include <errno.h>
#include "vec.h"
#include "qsort.h" #include "qsort.h"
extern "C" { extern "C" {
extern ssize_t read(int, void *, size_t); extern ssize_t read(int, void *, size_t);
extern ssize_t write(int, const void *, size_t); extern ssize_t write(int, const void *, size_t);
} }
namespace LA {
////////////////////////////////////////////////////////////////////////////// //////////////////////////////////////////////////////////////////////////////
//// forced instantization in the corespoding object file //// forced instantization in the corespoding object file
@ -42,13 +45,16 @@ template void NRVec<T>::get(int fd, bool dim, bool transp); \
INSTANTIZE(double) INSTANTIZE(double)
INSTANTIZE(complex<double>) INSTANTIZE(complex<double>)
INSTANTIZE(int)
INSTANTIZE(short)
INSTANTIZE(unsigned short)
INSTANTIZE(char) INSTANTIZE(char)
INSTANTIZE(short)
INSTANTIZE(int)
INSTANTIZE(long)
INSTANTIZE(long long)
INSTANTIZE(unsigned char) INSTANTIZE(unsigned char)
INSTANTIZE(unsigned short)
INSTANTIZE(unsigned int) INSTANTIZE(unsigned int)
INSTANTIZE(unsigned long) INSTANTIZE(unsigned long)
INSTANTIZE(unsigned long long)
@ -206,7 +212,7 @@ void NRVec< complex<double> >::axpy(const complex<double> alpha,
if (nn != x.nn) laerror("axpy of incompatible vectors"); if (nn != x.nn) laerror("axpy of incompatible vectors");
#endif #endif
copyonwrite(); copyonwrite();
cblas_zaxpy(nn, (void *)(&alpha), (void *)(x.v), 1, (void *)v, 1); cblas_zaxpy(nn, &alpha, x.v, 1, v, 1);
} }
// axpy call for T = double (strided) // axpy call for T = double (strided)
@ -223,7 +229,7 @@ void NRVec< complex<double> >::axpy(const complex<double> alpha,
const complex<double> *x, const int stride) const complex<double> *x, const int stride)
{ {
copyonwrite(); copyonwrite();
cblas_zaxpy(nn, (void *)(&alpha), (void *)x, stride, v, 1); cblas_zaxpy(nn, &alpha, x, stride, v, 1);
} }
// unary minus // unary minus
@ -242,7 +248,7 @@ NRVec< complex<double> >::operator-() const
{ {
NRVec< complex<double> > result(*this); NRVec< complex<double> > result(*this);
result.copyonwrite(); result.copyonwrite();
cblas_zdscal(nn, -1.0, (void *)(result.v), 1); cblas_zdscal(nn, -1.0, result.v, 1);
return result; return result;
} }
@ -279,37 +285,26 @@ template<>
NRVec< complex<double> > & NRVec< complex<double> >::normalize() NRVec< complex<double> > & NRVec< complex<double> >::normalize()
{ {
complex<double> tmp; complex<double> tmp;
tmp = cblas_dznrm2(nn, (void *)v, 1); tmp = cblas_dznrm2(nn, v, 1);
#ifdef DEBUG #ifdef DEBUG
if(!(tmp.real()) && !(tmp.imag())) laerror("normalization of zero vector"); if(!(tmp.real()) && !(tmp.imag())) laerror("normalization of zero vector");
#endif #endif
copyonwrite(); copyonwrite();
tmp = 1.0/tmp; tmp = 1.0/tmp;
cblas_zscal(nn, (void *)(&tmp), (void *)v, 1); cblas_zscal(nn, &tmp, v, 1);
return *this; return *this;
} }
//stubs for linkage //stubs for linkage
template<>
NRVec<unsigned char> & NRVec<unsigned char>::normalize() {laerror("normalize() impossible for integer types"); return *this;}
template<>
NRVec<int> & NRVec<int>::normalize() {laerror("normalize() impossible for integer types"); return *this;}
template<>
NRVec<short> & NRVec<short>::normalize() {laerror("normalize() impossible for integer types"); return *this;}
template<>
NRVec<char> & NRVec<char>::normalize() {laerror("normalize() impossible for integer types"); return *this;}
template<>
NRVec<unsigned long> & NRVec<unsigned long>::normalize() {laerror("normalize() impossible for integer types"); return *this;}
template<>
NRVec<unsigned int> & NRVec<unsigned int>::normalize() {laerror("normalize() impossible for integer types"); return *this;}
#define INSTANTIZE_DUMMY(T) \ #define INSTANTIZE_DUMMY(T) \
template<> void NRVec<T>::gemv(const T beta, const NRMat<T> &a, const char trans, const T alpha, const NRVec<T> &x) { laerror("gemv on unsupported types"); } \ template<> void NRVec<T>::gemv(const T beta, const NRMat<T> &a, const char trans, const T alpha, const NRVec<T> &x) { laerror("gemv on unsupported types"); } \
template<> void NRVec<T>::gemv(const T beta, const NRSMat<T> &a, const char trans, const T alpha, const NRVec<T> &x) { laerror("gemv on unsupported types"); } \ template<> void NRVec<T>::gemv(const T beta, const NRSMat<T> &a, const char trans, const T alpha, const NRVec<T> &x) { laerror("gemv on unsupported types"); } \
template<> void NRVec<T>::gemv(const T beta, const SparseMat<T> &a, const char trans, const T alpha, const NRVec<T> &x, bool s) { laerror("gemv on unsupported types"); } \ template<> void NRVec<T>::gemv(const T beta, const SparseMat<T> &a, const char trans, const T alpha, const NRVec<T> &x, bool s) { laerror("gemv on unsupported types"); } \
template<> void NRVec<T>::gemv(const LA_traits_complex<T>::Component_type beta, const LA_traits_complex<T>::NRMat_Noncomplex_type &a, const char trans, const LA_traits_complex<T>::Component_type alpha, const NRVec<T> &x) { laerror("gemv on unsupported types"); } \ template<> void NRVec<T>::gemv(const LA_traits_complex<T>::Component_type beta, const LA_traits_complex<T>::NRMat_Noncomplex_type &a, const char trans, const LA_traits_complex<T>::Component_type alpha, const NRVec<T> &x) { laerror("gemv on unsupported types"); } \
template<> void NRVec<T>::gemv(const LA_traits_complex<T>::Component_type beta, const LA_traits_complex<T>::NRSMat_Noncomplex_type &a, const char trans, const LA_traits_complex<T>::Component_type alpha, const NRVec<T> &x) { laerror("gemv on unsupported types"); } template<> void NRVec<T>::gemv(const LA_traits_complex<T>::Component_type beta, const LA_traits_complex<T>::NRSMat_Noncomplex_type &a, const char trans, const LA_traits_complex<T>::Component_type alpha, const NRVec<T> &x) { laerror("gemv on unsupported types"); } \
template<> NRVec<T> & NRVec<T>::normalize() {laerror("normalize() impossible for integer types"); return *this;} \
template<> const NRMat<T> NRVec<T>::otimes(const NRVec<T> &b,const bool conj, const T &scale) const {laerror("otimes presently implemented only for double and complex double"); return NRMat<T> ();}
@ -401,20 +396,22 @@ void NRVec< complex<double> >::gemv(const complex<double> beta,
// Direc product Mat = Vec | Vec // Direct product Mat = Vec | Vec
template<> template<>
const NRMat<double> NRVec<double>::operator|(const NRVec<double> &b) const const NRMat<double> NRVec<double>::otimes(const NRVec<double> &b,const bool conj, const double &scale) const
{ {
NRMat<double> result(0.,nn,b.nn); NRMat<double> result(0.,nn,b.nn);
cblas_dger(CblasRowMajor, nn, b.nn, 1., v, 1, b.v, 1, result, b.nn); cblas_dger(CblasRowMajor, nn, b.nn, scale, v, 1, b.v, 1, result, b.nn);
return result; return result;
} }
template<> template<>
const NRMat< complex<double> > const NRMat< complex<double> >
NRVec< complex<double> >::operator|(const NRVec< complex<double> > &b) const NRVec<complex<double> >::otimes(const NRVec< complex<double> > &b, const bool conj, const complex<double> &scale) const
{ {
NRMat< complex<double> > result(0.,nn,b.nn); NRMat< complex<double> > result(0.,nn,b.nn);
cblas_zgerc(CblasRowMajor, nn, b.nn, &CONE, v, 1, b.v, 1, result, b.nn); if(conj) cblas_zgerc(CblasRowMajor, nn, b.nn, &scale, v, 1, b.v, 1, result, b.nn);
else cblas_zgeru(CblasRowMajor, nn, b.nn, &scale, v, 1, b.v, 1, result, b.nn);
return result; return result;
} }
@ -436,24 +433,39 @@ else return memqsort<0,NRVec<T>,int,int>(*this,perm,from,to);
template class NRVec<double>; template class NRVec<double>;
template class NRVec<complex<double> >; template class NRVec<complex<double> >;
template class NRVec<char>; template class NRVec<char>;
template class NRVec<unsigned char>;
template class NRVec<short>; template class NRVec<short>;
template class NRVec<int>; template class NRVec<int>;
template class NRVec<long>;
template class NRVec<long long>;
template class NRVec<unsigned char>;
template class NRVec<unsigned short>;
template class NRVec<unsigned int>; template class NRVec<unsigned int>;
template class NRVec<unsigned long>; template class NRVec<unsigned long>;
template class NRVec<unsigned long long>;
INSTANTIZE_DUMMY(char) INSTANTIZE_DUMMY(char)
INSTANTIZE_DUMMY(unsigned char)
INSTANTIZE_DUMMY(short) INSTANTIZE_DUMMY(short)
INSTANTIZE_DUMMY(int) INSTANTIZE_DUMMY(int)
INSTANTIZE_DUMMY(long)
INSTANTIZE_DUMMY(long long)
INSTANTIZE_DUMMY(unsigned char)
INSTANTIZE_DUMMY(unsigned short)
INSTANTIZE_DUMMY(unsigned int) INSTANTIZE_DUMMY(unsigned int)
INSTANTIZE_DUMMY(unsigned long) INSTANTIZE_DUMMY(unsigned long)
INSTANTIZE_DUMMY(unsigned long long)
INSTANTIZE_DUMMY(complex<char>) INSTANTIZE_DUMMY(complex<char>)
INSTANTIZE_DUMMY(complex<unsigned char>)
INSTANTIZE_DUMMY(complex<short>) INSTANTIZE_DUMMY(complex<short>)
INSTANTIZE_DUMMY(complex<int>) INSTANTIZE_DUMMY(complex<int>)
INSTANTIZE_DUMMY(complex<long>)
INSTANTIZE_DUMMY(complex<long long>)
INSTANTIZE_DUMMY(complex<unsigned char>)
INSTANTIZE_DUMMY(complex<unsigned short>)
INSTANTIZE_DUMMY(complex<unsigned int>) INSTANTIZE_DUMMY(complex<unsigned int>)
INSTANTIZE_DUMMY(complex<unsigned long>)
INSTANTIZE_DUMMY(complex<unsigned long long>)
INSTANTIZE_DUMMY(complex<complex<double> >) INSTANTIZE_DUMMY(complex<complex<double> >)
INSTANTIZE_DUMMY(complex<complex<float> >) INSTANTIZE_DUMMY(complex<complex<float> >)
INSTANTIZE_DUMMY(complex<unsigned long>)
}//namespace

33
vec.h
View File

@ -21,6 +21,8 @@
#include "la_traits.h" #include "la_traits.h"
namespace LA {
////////////////////////////////////////////////////////////////////////////// //////////////////////////////////////////////////////////////////////////////
// Forward declarations // Forward declarations
template <typename T> void lawritemat(FILE *file,const T *a,int r,int c, template <typename T> void lawritemat(FILE *file,const T *a,int r,int c,
@ -104,7 +106,8 @@ public:
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;};
const NRMat<T> operator|(const NRVec<T> &rhs) const; const NRMat<T> otimes(const NRVec<T> &rhs, const bool conjugate=false, const T &scale=1) const; //outer product
inline const NRMat<T> operator|(const NRVec<T> &rhs) const {return otimes(rhs,true);};
inline const T sum() const {T sum=0; for(int i=0; i<nn; i++) sum += v[i]; return sum;}; //sum of its elements inline const T sum() const {T sum=0; for(int i=0; i<nn; i++) sum += v[i]; return sum;}; //sum of its elements
inline const T asum() const; //sum of its elements absolute values inline const T asum() const; //sum of its elements absolute values
inline const T dot(const T *a, const int stride=1) const; // ddot with a stride-vector inline const T dot(const T *a, const int stride=1) const; // ddot with a stride-vector
@ -136,26 +139,28 @@ public:
int sort(int direction=0, int from=0, int to= -1, int *perm=NULL); //sort, ascending by default, returns parity of permutation int sort(int direction=0, int from=0, int to= -1, int *perm=NULL); //sort, ascending by default, returns parity of permutation
}; };
}//namespace
//due to mutual includes this has to be after full class declaration //due to mutual includes this has to be after full class declaration
#include "mat.h" #include "mat.h"
#include "smat.h" #include "smat.h"
#include "sparsemat.h" #include "sparsemat.h"
namespace LA {
// formatted I/O // formatted I/O
template <typename T> template <typename T>
ostream & operator<<(ostream &s, const NRVec<T> &x) std::ostream & operator<<(std::ostream &s, const NRVec<T> &x)
{ {
int i, n; int i, n;
n = x.size(); n = x.size();
s << n << endl; s << n << std::endl;
for(i=0; i<n; i++) s << (typename LA_traits_io<T>::IOtype)x[i] << (i == n-1 ? '\n' : ' '); for(i=0; i<n; i++) s << (typename LA_traits_io<T>::IOtype)x[i] << (i == n-1 ? '\n' : ' ');
return s; return s;
} }
template <typename T> template <typename T>
istream & operator>>(istream &s, NRVec<T> &x) std::istream & operator>>(std::istream &s, NRVec<T> &x)
{ {
int i,n; int i,n;
@ -238,7 +243,7 @@ inline NRVec< complex<double> > &
NRVec< complex<double> >::operator+=(const complex<double> &a) NRVec< complex<double> >::operator+=(const complex<double> &a)
{ {
copyonwrite(); copyonwrite();
cblas_zaxpy(nn, (void *)(&CONE), (void *)(&a), 0, (void *)v, 1); cblas_zaxpy(nn, &CONE, &a, 0, v, 1);
return *this; return *this;
} }
@ -267,7 +272,7 @@ inline NRVec< complex<double> > &
NRVec< complex<double> >::operator-=(const complex<double> &a) NRVec< complex<double> >::operator-=(const complex<double> &a)
{ {
copyonwrite(); copyonwrite();
cblas_zaxpy(nn, (void *)(&CMONE), (void *)(&a), 0, (void *)v, 1); cblas_zaxpy(nn, &CMONE, &a, 0, v, 1);
return *this; return *this;
} }
@ -302,7 +307,7 @@ NRVec< complex<double> >::operator+=(const NRVec< complex<double> > &rhs)
if (nn != rhs.nn) laerror("daxpy of incompatible vectors"); if (nn != rhs.nn) laerror("daxpy of incompatible vectors");
#endif #endif
copyonwrite(); copyonwrite();
cblas_zaxpy(nn, (void *)(&CONE), rhs.v, 1, v, 1); cblas_zaxpy(nn, &CONE, rhs.v, 1, v, 1);
return *this; return *this;
} }
@ -366,7 +371,7 @@ NRVec< complex<double> >::operator-=(const NRVec< complex<double> > &rhs)
if (nn != rhs.nn) laerror("daxpy of incompatible vectors"); if (nn != rhs.nn) laerror("daxpy of incompatible vectors");
#endif #endif
copyonwrite(); copyonwrite();
cblas_zaxpy(nn, (void *)(&CMONE), (void *)rhs.v, 1, (void *)v, 1); cblas_zaxpy(nn, &CMONE, rhs.v, 1, v, 1);
return *this; return *this;
} }
@ -398,7 +403,7 @@ inline NRVec< complex<double> > &
NRVec< complex<double> >::operator*=(const complex<double> &a) NRVec< complex<double> >::operator*=(const complex<double> &a)
{ {
copyonwrite(); copyonwrite();
cblas_zscal(nn, (void *)(&a), (void *)v, 1); cblas_zscal(nn, &a, v, 1);
return *this; return *this;
} }
@ -432,7 +437,7 @@ NRVec< complex<double> >::operator*(const NRVec< complex<double> > &rhs) const
if (nn != rhs.nn) laerror("dot of incompatible vectors"); if (nn != rhs.nn) laerror("dot of incompatible vectors");
#endif #endif
complex<double> dot; complex<double> dot;
cblas_zdotc_sub(nn, (void *)v, 1, (void *)rhs.v, 1, (void *)(&dot)); cblas_zdotc_sub(nn, v, 1, rhs.v, 1, &dot);
return dot; return dot;
} }
@ -468,7 +473,7 @@ inline const complex<double>
NRVec< complex<double> >::dot(const complex<double> *y, const int stride) const NRVec< complex<double> >::dot(const complex<double> *y, const int stride) const
{ {
complex<double> dot; complex<double> dot;
cblas_zdotc_sub(nn, y, stride, v, 1, (void *)(&dot)); cblas_zdotc_sub(nn, y, stride, v, 1, &dot);
return dot; return dot;
} }
@ -527,7 +532,7 @@ inline const double NRVec<double>::norm() const
template<> template<>
inline const double NRVec< complex<double> >::norm() const inline const double NRVec< complex<double> >::norm() const
{ {
return cblas_dznrm2(nn, (void *)v, 1); return cblas_dznrm2(nn, v, 1);
} }
// Max element of the array // Max element of the array
@ -539,7 +544,7 @@ inline const double NRVec<double>::amax() const
template<> template<>
inline const complex<double> NRVec< complex<double> >::amax() const inline const complex<double> NRVec< complex<double> >::amax() const
{ {
return v[cblas_izamax(nn, (void *)v, 1)]; return v[cblas_izamax(nn, v, 1)];
} }
@ -687,6 +692,6 @@ for(int i=0; i<rhs.size(); ++i) r[i]=rhs[i];
return r; return r;
} }
}//namespace
#endif /* _LA_VEC_H_ */ #endif /* _LA_VEC_H_ */