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

*** empty log message ***

Browse Source
This commit is contained in:
jiri 2009-10-08 14:01:15 +00:00
parent c5309ee47b
commit 07c12d6896
15 changed files with 272 additions and 218 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

5
Makefile.am
View File

@ -12,11 +12,12 @@ EXTRA_DIST = LICENSE
CXXFLAGS += -fno-omit-frame-pointer -O3 -g -fPIC -finline-limit=1000 CXXFLAGS += -fno-omit-frame-pointer -O3 -g -fPIC -finline-limit=1000
CXXFLAGS += -DNO_STRASSEN -DFORTRAN_ CXXFLAGS += -DNO_STRASSEN -DFORTRAN_
#remove this for production code efficiency
CXXFLAGS += -DDEBUG CXXFLAGS += -DDEBUG
LDFLAGS += $(CBLASOPT) $(CLAPACKOPT) CXXFLAGS += $(CBLASOPT) $(CLAPACKOPT)
LDFLAGS += $(CBLASLIB) LDFLAGS += $(CBLASLIB)
LDFLAGS += $(TRACEBACKOPT) CXXFLAGS += $(TRACEBACKOPT)
LDFLAGS += $(TRACEBACKLIB) LDFLAGS += $(TRACEBACKLIB)

1
configure.ac
View File

@ -10,6 +10,7 @@ AM_INIT_AUTOMAKE([-Wall -Werror foreign])
AC_PROG_CXX AC_PROG_CXX
AC_PROG_CC AC_PROG_CC
AC_PROG_RANLIB AC_PROG_RANLIB
AC_LANG(C++)
# Checks for libraries. # Checks for libraries.
AC_CHECK_LIB([blas], [dgemm_]) AC_CHECK_LIB([blas], [dgemm_])

2
davidson.h
View File

@ -31,7 +31,7 @@
//therefore the whole implementation must be a template in a header //therefore the whole implementation must be a template in a header
//Note that for efficiency in a direct CI case the diagonalof() should cache its result //Note that for efficiency in a direct CI case the diagonalof() should cache its result
export template <typename T, typename Matrix> template <typename T, typename Matrix>
extern void davidson(const Matrix &bigmat, NRVec<T> &eivals, NRVec<T> *eivecs, const char *eivecsfile, extern void davidson(const Matrix &bigmat, NRVec<T> &eivals, NRVec<T> *eivecs, const char *eivecsfile,
int nroots=1, const bool verbose=0, const double eps=1e-6, int nroots=1, const bool verbose=0, const double eps=1e-6,
const bool incore=1, int maxit=100, const int maxkrylov = 500); const bool incore=1, int maxit=100, const int maxkrylov = 500);

1
la.h
View File

@ -19,6 +19,7 @@
#define _LA_H_ #define _LA_H_
//this should be the single include file for the end user //this should be the single include file for the end user
//
#include "vec.h" #include "vec.h"
#include "smat.h" #include "smat.h"

78
la_traits.h
View File

@ -15,6 +15,11 @@
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/>.
*/ */
//
//for autotools
//
#include "config.h"
//////////////////////////////////////////////////////////////////////////// ////////////////////////////////////////////////////////////////////////////
//LA traits classes and generally needed includes //LA traits classes and generally needed includes
@ -22,6 +27,8 @@
#define _LA_TRAITS_INCL #define _LA_TRAITS_INCL
extern bool _LA_count_check;
using namespace std; using namespace std;
#include <stdio.h> #include <stdio.h>
@ -38,6 +45,14 @@ extern "C" {
} }
#endif #endif
#ifdef NONCLAPACK
#include "noncblas.h"
#else
extern "C" {
#include "clapack.h"
}
#endif
//forward declarations //forward declarations
template<typename C> class NRVec; template<typename C> class NRVec;
@ -47,6 +62,44 @@ 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;
typedef class {} Dummy_type;
typedef class {} Dummy_type2;
//for components of complex numbers
//
template<typename C>
struct LA_traits_complex
{
typedef Dummy_type Component_type;
typedef Dummy_type NRVec_Noncomplex_type;
typedef Dummy_type NRMat_Noncomplex_type;
typedef Dummy_type2 NRSMat_Noncomplex_type;
};
#define SPECIALIZE_COMPLEX(T) \
template<> \
struct LA_traits_complex<complex<T> > \
{ \
typedef T Component_type; \
typedef NRVec<T> NRVec_Noncomplex_type; \
typedef NRMat<T> NRMat_Noncomplex_type; \
typedef NRSMat<T> NRSMat_Noncomplex_type; \
};
SPECIALIZE_COMPLEX(double)
SPECIALIZE_COMPLEX(complex<double>)
SPECIALIZE_COMPLEX(float)
SPECIALIZE_COMPLEX(complex<float>)
SPECIALIZE_COMPLEX(char)
SPECIALIZE_COMPLEX(unsigned char)
SPECIALIZE_COMPLEX(short)
SPECIALIZE_COMPLEX(int)
SPECIALIZE_COMPLEX(unsigned int)
SPECIALIZE_COMPLEX(unsigned long)
//for general sortable classes //for general sortable classes
template<typename T, typename I, int type> struct LA_sort_traits; template<typename T, typename I, int type> struct LA_sort_traits;
@ -96,7 +149,12 @@ class isscalar { public: typedef scalar_false scalar_type;};
//specializations //specializations
#define SCALAR(X) \ #define SCALAR(X) \
template<>\ template<>\
class isscalar<X> {public: typedef scalar_true scalar_type;}; class isscalar<X> {public: typedef scalar_true scalar_type;};\
template<>\
class isscalar<complex<X> > {public: typedef scalar_true scalar_type;};\
template<>\
class isscalar<complex<complex<X> > > {public: typedef scalar_true scalar_type;};\
//declare what is scalar //declare what is scalar
SCALAR(char) SCALAR(char)
@ -111,21 +169,24 @@ SCALAR(unsigned long)
SCALAR(unsigned long long) SCALAR(unsigned long long)
SCALAR(float) SCALAR(float)
SCALAR(double) SCALAR(double)
SCALAR(complex<float>)
SCALAR(complex<double>)
SCALAR(void *) SCALAR(void *)
#undef SCALAR #undef SCALAR
//now declare the traits for scalars and for composed classes //declare this generically as traits for any unknown class
//NOTE! methods in traits classes have to be declared static, template<typename C, typename Scalar> struct LA_traits_aux
//since the class itself is never instantiated. {
//for performance, it can be also inlined at the same time typedef Dummy_type normtype;
template<typename C, typename Scalar> struct LA_traits_aux {}; };
//TRAITS SPECIALIZATIONS //TRAITS SPECIALIZATIONS
////now declare the traits for scalars and for composed classes
////NOTE! methods in traits classes have to be declared static,
////since the class itself is never instantiated.
////for performance, it can be also inlined at the same time
//
//complex scalars //complex scalars
template<typename C> template<typename C>
@ -226,6 +287,7 @@ static void copyonwrite(X<C> &x) {x.copyonwrite();} \
generate_traits_smat(NRSMat) generate_traits_smat(NRSMat)
generate_traits_smat(NRSMat_from1) generate_traits_smat(NRSMat_from1)
//the final traits class //the final traits class
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> {};

3
laerror.cc
View File

@ -30,6 +30,9 @@
#include "traceback.h" #include "traceback.h"
#endif #endif
bool _LA_count_check=true;
extern "C" void _findme(void) {}; //for autoconf test we need a function with C linkage
void laerror(const char *s1) void laerror(const char *s1)
{ {

30
mat.cc
View File

@ -369,6 +369,25 @@ if(n==0) n=nn;
return *this; return *this;
} }
//complex from real
template<>
NRMat<complex<double> >::NRMat(const NRMat<double> &rhs, bool imagpart)
: nn(rhs.nrows()), mm(rhs.ncols()), count(new int(1))
{
#ifdef MATPTR
v = new complex<double>*[n];
v[0] = new complex<double>[mm*nn];
for (int i=1; i<n; i++) v[i] = v[i-1] + m;
memset(v[0], 0, nn*mm*sizeof(complex<double>));
cblas_dcopy(nn*mm,&rhs[0][0],1,((double *)v[0]) + (imagpart?1:0),2);
#else
v = new complex<double>[mm*nn];
memset(v, 0, nn*mm*sizeof(complex<double>));
cblas_dcopy(nn*mm,&rhs[0][0],1,((double *)v) + (imagpart?1:0),2);
#endif
}
// Output of Mat // Output of Mat
template <typename T> template <typename T>
void NRMat<T>::fprintf(FILE *file, const char *format, const int modulo) const void NRMat<T>::fprintf(FILE *file, const char *format, const int modulo) const
@ -501,6 +520,17 @@ for(int i=0; i<nn; ++i)
(*this)(i,j) = x*(2.*random()/(1.+RAND_MAX) -1.); (*this)(i,j) = x*(2.*random()/(1.+RAND_MAX) -1.);
} }
template<>
void NRMat<complex<double> >::randomize(const double &x)
{
for(int i=0; i<nn; ++i)
for(int j=0; j<mm; ++j)
{
(*this)(i,j).real() = x*(2.*random()/(1.+RAND_MAX) -1.);
(*this)(i,j).imag() = x*(2.*random()/(1.+RAND_MAX) -1.);
}
}
// Mat *= a // Mat *= a

18
matexp.h
View File

@ -24,14 +24,6 @@
#include "la_traits.h" #include "la_traits.h"
#include "laerror.h" #include "laerror.h"
#ifndef NONCBLAS
extern "C" {
#include "cblas.h"
#include "clapack.h"
}
#else
#include "noncblas.h"
#endif
template<class T,class R> template<class T,class R>
@ -184,7 +176,7 @@ return int(ceil(log(n)/log2-log(.75)));
template<class T, class C> template<class T, class C>
NRVec<C> exp_aux(const T &x, int &power,int maxpower= -1, int maxtaylor= -1) NRVec<C> exp_aux(const T &x, int &power,int maxpower= -1, int maxtaylor= -1, C prescale=1.)
{ {
//should better be computed by mathematica to have accurate last digits, chebyshev instead, see exp in glibc //should better be computed by mathematica to have accurate last digits, chebyshev instead, see exp in glibc
static double exptaylor[]={ static double exptaylor[]={
@ -210,7 +202,7 @@ static double exptaylor[]={
8.2206352466243294955e-18, 8.2206352466243294955e-18,
4.1103176233121648441e-19, 4.1103176233121648441e-19,
0.}; 0.};
double mnorm= x.norm(); double mnorm= x.norm() * 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=exp(-log(2.)*power);
@ -268,14 +260,14 @@ return r;
//and probably not efficient either //and probably not efficient either
template<class M, class V> template<class M, class V>
void exptimesdestructive(const M &mat, V &result, V &rhs, bool transpose=false, const double scale=1., int maxpower= -1, int maxtaylor= -1, bool mat_is_0=false) //uses just matrix vector multiplication void exptimesdestructive(const M &mat, V &result, V &rhs, bool transpose=false, const typename LA_traits<V>::elementtype scale=1., int maxpower= -1, int maxtaylor= -1, bool mat_is_0=false) //uses just matrix vector multiplication
{ {
if(mat_is_0) {result = rhs; LA_traits<V>::copyonwrite(result); return;} //prevent returning a shallow copy of rhs if(mat_is_0) {result = rhs; LA_traits<V>::copyonwrite(result); return;} //prevent returning a shallow copy of rhs
if(mat.nrows()!=mat.ncols()||(unsigned int) mat.nrows() != (unsigned int)rhs.size()) laerror("inappropriate sizes in exptimes"); if(mat.nrows()!=mat.ncols()||(unsigned int) mat.nrows() != (unsigned int)rhs.size()) laerror("inappropriate sizes in exptimes");
int power; int power;
//prepare the polynom of and effectively scale the matrix //prepare the polynom of and effectively scale the matrix
NRVec<typename LA_traits<V>::elementtype> taylor2=exp_aux<M,typename LA_traits<V>::elementtype>(mat,power,maxpower,maxtaylor); NRVec<typename LA_traits<V>::elementtype> taylor2=exp_aux<M,typename LA_traits<V>::elementtype>(mat,power,maxpower,maxtaylor,scale);
V tmp; V tmp;
@ -298,7 +290,7 @@ return;
template<class M, class V> template<class M, class V>
const V exptimes(const M &mat, V rhs, bool transpose=false, const double scale=1., int maxpower= -1, int maxtaylor= -1, bool mat_is_0=false ) const V exptimes(const M &mat, V rhs, bool transpose=false, const typename LA_traits<V>::elementtype scale=1., int maxpower= -1, int maxtaylor= -1, bool mat_is_0=false )
{ {
V result; V result;
exptimesdestructive(mat,result,rhs,transpose,scale,maxpower,maxtaylor,mat_is_0); exptimesdestructive(mat,result,rhs,transpose,scale,maxpower,maxtaylor,mat_is_0);

9
nonclass.cc
View File

@ -16,15 +16,8 @@
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 NONCBLAS //this can be safely included since it contains ifdefs NONCBLAS and NONCLAPACK anyway
extern "C" {
#include "cblas.h"
#include "clapack.h"
}
#else
#include "noncblas.h" #include "noncblas.h"
#endif
#include "vec.h" #include "vec.h"
#include "smat.h" #include "smat.h"
#include "mat.h" #include "mat.h"

4
nonclass.h
View File

@ -25,7 +25,7 @@
//MISC //MISC
export template <class T> template <class T>
const NRSMat<T> twoside_transform(const NRSMat<T> &S, const NRMat<T> &C, bool transp=0) //calculate C^dagger S C const NRSMat<T> twoside_transform(const NRSMat<T> &S, const NRMat<T> &C, bool transp=0) //calculate C^dagger S C
{ {
if(transp) if(transp)
@ -44,7 +44,7 @@ return NRSMat<T>(result);
export template <class T> template <class T>
const NRMat<T> diagonalmatrix(const NRVec<T> &x) const NRMat<T> diagonalmatrix(const NRVec<T> &x)
{ {
int n=x.size(); int n=x.size();

18
smat.cc
View File

@ -138,6 +138,14 @@ void NRSMat<double>::randomize(const double &x)
for(int i=0; i<NN2; ++i) v[i] = x*(2.*random()/(1.+RAND_MAX) -1.); for(int i=0; i<NN2; ++i) v[i] = x*(2.*random()/(1.+RAND_MAX) -1.);
} }
template<>
void NRSMat<complex<double> >::randomize(const double &x)
{
for(int i=0; i<NN2; ++i) v[i].real() = x*(2.*random()/(1.+RAND_MAX) -1.);
for(int i=0; i<NN2; ++i) v[i].imag() = x*(2.*random()/(1.+RAND_MAX) -1.);
for(int i=0; i<nn; ++i) for(int j=0; j<=i; ++j) if(i==j) v[i*(i+1)/2+j].imag()=0; //hermitean
}
// write matrix to the file with specific format // write matrix to the file with specific format
@ -383,6 +391,16 @@ void NRSMat< complex<double> >::axpy(const complex<double> alpha,
cblas_zaxpy(nn, (void *)(&alpha), (void *)x.v, 1, (void *)v, 1); cblas_zaxpy(nn, (void *)(&alpha), (void *)x.v, 1, (void *)v, 1);
} }
//complex from real
template<>
NRSMat<complex<double> >::NRSMat(const NRSMat<double> &rhs, bool imagpart)
: nn(rhs.nrows()), v(new complex<double>[rhs.nrows()*(rhs.nrows()+1)/2]), count(new int(1))
{
memset(v,0,nn*(nn+1)/2*sizeof(complex<double>));
cblas_dcopy(nn*(nn+1)/2,&rhs(0,0),1,((double *)v) + (imagpart?1:0),2);
}
////////////////////////////////////////////////////////////////////////////// //////////////////////////////////////////////////////////////////////////////

62
sparsemat.cc
View File

@ -43,7 +43,7 @@ extern ssize_t write(int, const void *, size_t);
} }
export template <class T> template <class T>
void SparseMat<T>::get(int fd, bool dimen, bool transp) void SparseMat<T>::get(int fd, bool dimen, bool transp)
{ {
errno=0; errno=0;
@ -75,7 +75,7 @@ list=l;
} }
export template <class T> template <class T>
void SparseMat<T>::put(int fd,bool dimen, bool transp) const void SparseMat<T>::put(int fd,bool dimen, bool transp) const
{ {
errno=0; errno=0;
@ -103,7 +103,7 @@ if(2*sizeof(SPMatindex) != write(fd,&sentinel,2*sizeof(SPMatindex))) laerror("ca
//helpers to be used from different functions //helpers to be used from different functions
export template <class T> template <class T>
void SparseMat<T>::unsort() void SparseMat<T>::unsort()
{ {
if(symmetric) colsorted=NULL; if(symmetric) colsorted=NULL;
@ -113,7 +113,7 @@ colsorted=rowsorted=NULL;
nonzero=0; nonzero=0;
} }
export template <class T> template <class T>
void SparseMat<T>::deletelist() void SparseMat<T>::deletelist()
{ {
if(colsorted||rowsorted) unsort();//prevent obsolete pointers if(colsorted||rowsorted) unsort();//prevent obsolete pointers
@ -131,7 +131,7 @@ count=NULL;
} }
//no checks, not to be public //no checks, not to be public
export template <class T> template <class T>
void SparseMat<T>::copylist(const matel<T> *l) void SparseMat<T>::copylist(const matel<T> *l)
{ {
list=NULL; list=NULL;
@ -142,7 +142,7 @@ while(l)
} }
} }
export template <class T> template <class T>
void SparseMat<T>::copyonwrite() void SparseMat<T>::copyonwrite()
{ {
if(!count) laerror("probably an assignment to undefined sparse matrix"); if(!count) laerror("probably an assignment to undefined sparse matrix");
@ -241,7 +241,7 @@ else
} }
export template <class T> template <class T>
unsigned int SparseMat<T>::length() const unsigned int SparseMat<T>::length() const
{ {
if(nonzero) return nonzero; if(nonzero) return nonzero;
@ -258,7 +258,7 @@ return n;
} }
export template <class T> template <class T>
unsigned int SparseMat<T>::sort(int type) const //must be const since used from operator* which must be const to be compatible with other stuff, dirty casts here unsigned int SparseMat<T>::sort(int type) const //must be const since used from operator* which must be const to be compatible with other stuff, dirty casts here
{ {
if(type==0&&rowsorted || type==1&&colsorted) return nonzero; if(type==0&&rowsorted || type==1&&colsorted) return nonzero;
@ -295,7 +295,7 @@ return nonzero; //number of (in principle) nonzero elements
} }
export template <class T> template <class T>
void SparseMat<T>::simplify() void SparseMat<T>::simplify()
{ {
unsigned int n; unsigned int n;
@ -359,7 +359,7 @@ unsort(); //since there were NULLs introduced, rowsorted is not dense
} }
export template <class T> template <class T>
void SparseMat<T>::resize(const SPMatindex n, const SPMatindex m) void SparseMat<T>::resize(const SPMatindex n, const SPMatindex m)
{ {
unsort(); unsort();
@ -378,7 +378,7 @@ void SparseMat<T>::resize(const SPMatindex n, const SPMatindex m)
colsorted=rowsorted=NULL; colsorted=rowsorted=NULL;
} }
export template <class T> template <class T>
void SparseMat<T>::incsize(const SPMatindex n, const SPMatindex m) void SparseMat<T>::incsize(const SPMatindex n, const SPMatindex m)
{ {
if(symmetric && n!=m) laerror("unsymmetric size increment of a symmetric sparsemat"); if(symmetric && n!=m) laerror("unsymmetric size increment of a symmetric sparsemat");
@ -391,7 +391,7 @@ void SparseMat<T>::incsize(const SPMatindex n, const SPMatindex m)
export template <class T> template <class T>
void SparseMat<T>::addsafe(const SPMatindex n, const SPMatindex m, const T elem) void SparseMat<T>::addsafe(const SPMatindex n, const SPMatindex m, const T elem)
{ {
#ifdef debug #ifdef debug
@ -409,7 +409,7 @@ add(n,m,elem);
//assignment operator //assignment operator
export template <class T> template <class T>
SparseMat<T> & SparseMat<T>::operator=(const SparseMat<T> &rhs) SparseMat<T> & SparseMat<T>::operator=(const SparseMat<T> &rhs)
{ {
if (this != &rhs) if (this != &rhs)
@ -427,7 +427,7 @@ SparseMat<T> & SparseMat<T>::operator=(const SparseMat<T> &rhs)
return *this; return *this;
} }
export template <class T> template <class T>
SparseMat<T> & SparseMat<T>::join(SparseMat<T> &rhs) SparseMat<T> & SparseMat<T>::join(SparseMat<T> &rhs)
{ {
if(symmetric!=rhs.symmetric||nn!=rhs.nn||mm!=rhs.mm) laerror("incompatible matrices in join()"); if(symmetric!=rhs.symmetric||nn!=rhs.nn||mm!=rhs.mm) laerror("incompatible matrices in join()");
@ -442,7 +442,7 @@ return *this;
} }
export template <class T> template <class T>
SparseMat<T> & SparseMat<T>::addtriangle(const SparseMat &rhs, const bool lower, const char sign) SparseMat<T> & SparseMat<T>::addtriangle(const SparseMat &rhs, const bool lower, const char sign)
{ {
if(nn!=rhs.nn||mm!=rhs.mm) laerror("incompatible dimensions for +="); if(nn!=rhs.nn||mm!=rhs.mm) laerror("incompatible dimensions for +=");
@ -461,7 +461,7 @@ while(l)
return *this; return *this;
} }
export template <class T> template <class T>
SparseMat<T> & SparseMat<T>::operator+=(const SparseMat<T> &rhs) SparseMat<T> & SparseMat<T>::operator+=(const SparseMat<T> &rhs)
{ {
if(symmetric&&!rhs.symmetric) laerror("cannot add general to symmetric sparse"); if(symmetric&&!rhs.symmetric) laerror("cannot add general to symmetric sparse");
@ -491,7 +491,7 @@ while(l)
return *this; return *this;
} }
export template <class T> template <class T>
SparseMat<T> & SparseMat<T>::operator-=(const SparseMat<T> &rhs) SparseMat<T> & SparseMat<T>::operator-=(const SparseMat<T> &rhs)
{ {
if(symmetric&&!rhs.symmetric) laerror("cannot add general to symmetric sparse"); if(symmetric&&!rhs.symmetric) laerror("cannot add general to symmetric sparse");
@ -523,7 +523,7 @@ return *this;
//constructor from a dense matrix //constructor from a dense matrix
export template <class T> template <class T>
SparseMat<T>::SparseMat(const NRMat<T> &rhs) SparseMat<T>::SparseMat(const NRMat<T> &rhs)
{ {
nn=rhs.nrows(); nn=rhs.nrows();
@ -583,7 +583,7 @@ return divide?NULL:&r[0];
//constructor dense matrix from sparse //constructor dense matrix from sparse
export template <class T> template <class T>
NRMat<T>::NRMat(const SparseMat<T> &rhs) NRMat<T>::NRMat(const SparseMat<T> &rhs)
{ {
nn=rhs.nrows(); nn=rhs.nrows();
@ -618,7 +618,7 @@ while(l)
//constructor dense symmetric packed matrix from sparse //constructor dense symmetric packed matrix from sparse
#define nn2 (nn*(nn+1)/2) #define nn2 (nn*(nn+1)/2)
export template <class T> template <class T>
NRSMat<T>::NRSMat(const SparseMat<T> &rhs) NRSMat<T>::NRSMat(const SparseMat<T> &rhs)
{ {
if(!rhs.issymmetric()||rhs.nrows()!=rhs.ncols()) laerror("sparse matrix is not symmetric"); if(!rhs.issymmetric()||rhs.nrows()!=rhs.ncols()) laerror("sparse matrix is not symmetric");
@ -636,7 +636,7 @@ while(l)
#undef nn2 #undef nn2
//constructor dense vector from sparse //constructor dense vector from sparse
export template <class T> template <class T>
NRVec<T>::NRVec(const SparseMat<T> &rhs) NRVec<T>::NRVec(const SparseMat<T> &rhs)
{ {
if(rhs.nrows()>1 && rhs.ncols()>1) laerror("cannot construct a vector from a sparse matrix with more than one row/column"); if(rhs.nrows()>1 && rhs.ncols()>1) laerror("cannot construct a vector from a sparse matrix with more than one row/column");
@ -659,7 +659,7 @@ else while(l)
} }
//assignment of a scalar matrix //assignment of a scalar matrix
export 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");
@ -675,7 +675,7 @@ for(i=0;i<nn;++i) add(i,i,a);
return *this; return *this;
} }
export 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");
@ -686,7 +686,7 @@ for(i=0;i<nn;++i) add(i,i,a);
return *this; return *this;
} }
export 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");
@ -700,7 +700,7 @@ return *this;
//constructor from a dense symmetric matrix //constructor from a dense symmetric matrix
export template <class T> template <class T>
SparseMat<T>::SparseMat(const NRSMat<T> &rhs) SparseMat<T>::SparseMat(const NRSMat<T> &rhs)
{ {
nn=rhs.nrows(); nn=rhs.nrows();
@ -724,7 +724,7 @@ for(i=0;i<nn;++i)
} }
} }
export template <class T> template <class T>
void SparseMat<T>::transposeme() void SparseMat<T>::transposeme()
{ {
if(!count) laerror("transposeme on undefined lhs"); if(!count) laerror("transposeme on undefined lhs");
@ -739,7 +739,7 @@ while(l)
SWAP(nn,mm); SWAP(nn,mm);
} }
export template <class T> template <class T>
void SparseMat<T>::setunsymmetric() void SparseMat<T>::setunsymmetric()
{ {
if(!symmetric) return; if(!symmetric) return;
@ -760,7 +760,7 @@ while(l) //include the mirror picture of elements into the list
} }
export 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");
@ -1217,7 +1217,7 @@ simplify();
//direct sum and product -- only partly implemented at the moment //direct sum and product -- only partly implemented at the moment
export template<typename T> template<typename T>
SparseMat<T> & SparseMat<T>::oplusequal(const NRMat<T> &rhs) SparseMat<T> & SparseMat<T>::oplusequal(const NRMat<T> &rhs)
{ {
if(issymmetric()) laerror("oplusequal symmetric-unsymmetric"); if(issymmetric()) laerror("oplusequal symmetric-unsymmetric");
@ -1237,7 +1237,7 @@ return *this;
export template<typename T> template<typename T>
SparseMat<T> & SparseMat<T>::oplusequal(const NRSMat<T> &rhs) SparseMat<T> & SparseMat<T>::oplusequal(const NRSMat<T> &rhs)
{ {
if(!issymmetric()) laerror("oplusequal symmetric-unsymmetric"); if(!issymmetric()) laerror("oplusequal symmetric-unsymmetric");
@ -1257,7 +1257,7 @@ return *this;
export template <class T> template <class T>
SparseMat<T> & SparseMat<T>::oplusequal(const SparseMat<T> &rhs) SparseMat<T> & SparseMat<T>::oplusequal(const SparseMat<T> &rhs)
{ {
if(symmetric != rhs.symmetric) laerror("incompatible symmetry of sparsemats in oplusequal"); if(symmetric != rhs.symmetric) laerror("incompatible symmetry of sparsemats in oplusequal");

3
strassen.cc
View File

@ -27,8 +27,11 @@ extern "C" void fmm(const char c_transa,const char c_transb,const int m,const in
double *d_aux,int i_naux); double *d_aux,int i_naux);
extern "C" void strassen_cutoff(int c, int c1, int c2, int c3); extern "C" void strassen_cutoff(int c, int c1, int c2, int c3);
template<>
void NRMat<double>::s_cutoff(const int c, const int c1, const int c2, const int c3) const void NRMat<double>::s_cutoff(const int c, const int c1, const int c2, const int c3) const
{ strassen_cutoff(c,c1,c2,c3);} { strassen_cutoff(c,c1,c2,c3);}
template<>
void NRMat<double>::strassen(const double beta, const NRMat<double> &a, const char transa, const NRMat<double> &b, const char transb, const double alpha) void NRMat<double>::strassen(const double beta, const NRMat<double> &a, const char transa, const NRMat<double> &b, const char transb, const double alpha)
{ {
int l(transa=='n'?a.nn:a.mm); int l(transa=='n'?a.nn:a.mm);

38
t.cc
View File

@ -1172,7 +1172,7 @@ cout <<a.otimes(b);
} }
//test of general square matrix eigenvector derivatives //test of general square matrix eigenvector derivatives
if(1) if(0)
{ {
const bool biorthonormalize=1; const bool biorthonormalize=1;
@ -1270,9 +1270,45 @@ cout <<"right eigenvector analytic derivative\n"<<vrg<<endl;
cout <<"right eigenvector derivative error = "<<(vrd-vrg).norm()<<endl; cout <<"right eigenvector derivative error = "<<(vrd-vrg).norm()<<endl;
}
if(0)
{
try { laerror("test catch exception"); }
catch(LAerror x)
{
cout <<"caught exception: "<<x <<endl;
}
laerror("test exception 2");
}
if(0)
{
NRVec<double> v(3);
v[0]=1; v[1]=2; v[2]=3;
NRVec<complex<double> > vv = v;
NRVec<double>u(v);
vv += u;
cout <<vv;
}
if(1)
{
complex<double> scale; cin >> scale;
NRMat<complex<double> > h; cin >>h;
NRVec<complex<double> > x(h.nrows());
NRVec<complex<double> > y,z;
x.randomize(1.);
y=exptimes(h*scale,x,false,1.);
z=exptimes(h,x,false,scale);
cout <<x;
cout <<y;
cout <<z;
cout <<"Error "<<(y-z).norm()<<endl;
} }

218
vec.cc
View File

@ -168,6 +168,23 @@ void NRVec<double>::randomize(const double &x)
for(int i=0; i<nn; ++i) v[i] = x*(2.*random()/(1.+RAND_MAX) -1.); for(int i=0; i<nn; ++i) v[i] = x*(2.*random()/(1.+RAND_MAX) -1.);
} }
template<>
void NRVec<complex<double> >::randomize(const double &x)
{
for(int i=0; i<nn; ++i) v[i] = complex<double> (x*(2.*random()/(1.+RAND_MAX) -1.),x*(2.*random()/(1.+RAND_MAX) -1.));
}
//complex from real constructor
template<>
NRVec<complex<double> >::NRVec(const NRVec<double> &rhs, bool imagpart)
: nn(rhs.size()), v(new complex<double>[rhs.size()]), count(new int(1))
{
memset(v,0,nn*sizeof(complex<double>));
cblas_dcopy(nn,&rhs[0],1,((double *)v) + (imagpart?1:0),2);
}
// axpy call for T = double (not strided) // axpy call for T = double (not strided)
template<> template<>
@ -287,158 +304,12 @@ template<>
NRVec<unsigned int> & NRVec<unsigned int>::normalize() {laerror("normalize() impossible for integer types"); return *this;} NRVec<unsigned int> & NRVec<unsigned int>::normalize() {laerror("normalize() impossible for integer types"); return *this;}
template<> #define INSTANTIZE_DUMMY(T) \
void NRVec<int>::gemv(const int beta, 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"); } \
const NRSMat<int> &A, const char trans, 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"); } \
const int alpha, const NRVec &x) 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"); } \
laerror("not yet implemented"); 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<short>::gemv(const short beta,
const NRSMat<short> &A, const char trans,
const short alpha, const NRVec &x)
{
laerror("not yet implemented");
}
template<>
void NRVec<unsigned long>::gemv(const unsigned long beta,
const NRSMat<unsigned long> &A, const char trans,
const unsigned long alpha, const NRVec &x)
{
laerror("not yet implemented");
}
template<>
void NRVec<unsigned int>::gemv(const unsigned int beta,
const NRSMat<unsigned int> &A, const char trans,
const unsigned int alpha, const NRVec &x)
{
laerror("not yet implemented");
}
template<>
void NRVec<unsigned char>::gemv(const unsigned char beta,
const NRSMat<unsigned char> &A, const char trans,
const unsigned char alpha, const NRVec &x)
{
laerror("not yet implemented");
}
template<>
void NRVec<char>::gemv(const char beta,
const NRSMat<char> &A, const char trans,
const char alpha, const NRVec &x)
{
laerror("not yet implemented");
}
template<>
void NRVec<int>::gemv(const int beta,
const NRMat<int> &A, const char trans,
const int alpha, const NRVec &x)
{
laerror("not yet implemented");
}
template<>
void NRVec<short>::gemv(const short beta,
const NRMat<short> &A, const char trans,
const short alpha, const NRVec &x)
{
laerror("not yet implemented");
}
template<>
void NRVec<unsigned long>::gemv(const unsigned long beta,
const NRMat<unsigned long> &A, const char trans,
const unsigned long alpha, const NRVec &x)
{
laerror("not yet implemented");
}
template<>
void NRVec<unsigned int>::gemv(const unsigned int beta,
const NRMat<unsigned int> &A, const char trans,
const unsigned int alpha, const NRVec &x)
{
laerror("not yet implemented");
}
template<>
void NRVec<unsigned char>::gemv(const unsigned char beta,
const NRMat<unsigned char> &A, const char trans,
const unsigned char alpha, const NRVec &x)
{
laerror("not yet implemented");
}
template<>
void NRVec<char>::gemv(const char beta,
const NRMat<char> &A, const char trans,
const char alpha, const NRVec &x)
{
laerror("not yet implemented");
}
template<>
void NRVec<int>::gemv(const int beta,
const SparseMat<int> &A, const char trans,
const int alpha, const NRVec &x, bool s)
{
laerror("not yet implemented");
}
template<>
void NRVec<short>::gemv(const short beta,
const SparseMat<short> &A, const char trans,
const short alpha, const NRVec &x, bool s)
{
laerror("not yet implemented");
}
template<>
void NRVec<char>::gemv(const char beta,
const SparseMat<char> &A, const char trans,
const char alpha, const NRVec &x, bool s)
{
laerror("not yet implemented");
}
template<>
void NRVec<unsigned long>::gemv(const unsigned long beta,
const SparseMat<unsigned long> &A, const char trans,
const unsigned long alpha, const NRVec &x, bool s)
{
laerror("not yet implemented");
}
template<>
void NRVec<unsigned int>::gemv(const unsigned int beta,
const SparseMat<unsigned int> &A, const char trans,
const unsigned int alpha, const NRVec &x, bool s)
{
laerror("not yet implemented");
}
template<>
void NRVec<unsigned char>::gemv(const unsigned char beta,
const SparseMat<unsigned char> &A, const char trans,
const unsigned char alpha, const NRVec &x, bool s)
{
laerror("not yet implemented");
}
@ -457,6 +328,21 @@ void NRVec<double>::gemv(const double beta, const NRMat<double> &A,
A.nrows(), A.ncols(), alpha, A, A.ncols(), x.v, 1, beta, v, 1); A.nrows(), A.ncols(), alpha, A, A.ncols(), x.v, 1, beta, v, 1);
} }
template<>
void NRVec<complex<double> >::gemv(const double beta, const NRMat<double> &A,
const char trans, const double alpha, const NRVec<complex<double> > &x)
{
#ifdef DEBUG
if ((trans == 'n'?A.ncols():A.nrows()) != x.size())
laerror("incompatible sizes in gemv A*x");
#endif
copyonwrite();
cblas_dgemv(CblasRowMajor, (trans=='n' ? CblasNoTrans:CblasTrans), A.nrows(), A.ncols(), alpha, A, A.ncols(), (double *)x.v, 2, beta, (double *)v, 2);
cblas_dgemv(CblasRowMajor, (trans=='n' ? CblasNoTrans:CblasTrans), A.nrows(), A.ncols(), alpha, A, A.ncols(), ((double *)x.v) + 1, 2, beta, ((double *)v)+1, 2);
}
template<> template<>
void NRVec< complex<double> >::gemv(const complex<double> beta, void NRVec< complex<double> >::gemv(const complex<double> beta,
const NRMat< complex<double> > &A, const char trans, const NRMat< complex<double> > &A, const char trans,
@ -484,6 +370,19 @@ void NRVec<double>::gemv(const double beta, const NRSMat<double> &A,
cblas_dspmv(CblasRowMajor, CblasLower, A.ncols(), alpha, A, x.v, 1, beta, v, 1); cblas_dspmv(CblasRowMajor, CblasLower, A.ncols(), alpha, A, x.v, 1, beta, v, 1);
} }
template<>
void NRVec<complex<double> >::gemv(const double beta, const NRSMat<double> &A,
const char trans, const double alpha, const NRVec<complex<double> > &x)
{
#ifdef DEBUG
if (A.ncols()!=x.size()) laerror("incompatible dimension in gemv A*x");
#endif
copyonwrite();
cblas_dspmv(CblasRowMajor, CblasLower, A.ncols(), alpha, A, (double *)x.v, 2, beta, (double *)v, 2);
cblas_dspmv(CblasRowMajor, CblasLower, A.ncols(), alpha, A, ((double *)x.v)+1, 2, beta, ((double *)v)+1, 2);
}
template<> template<>
void NRVec< complex<double> >::gemv(const complex<double> beta, void NRVec< complex<double> >::gemv(const complex<double> beta,
@ -543,3 +442,18 @@ template class NRVec<int>;
template class NRVec<unsigned int>; template class NRVec<unsigned int>;
template class NRVec<unsigned long>; template class NRVec<unsigned long>;
INSTANTIZE_DUMMY(char)
INSTANTIZE_DUMMY(unsigned char)
INSTANTIZE_DUMMY(short)
INSTANTIZE_DUMMY(int)
INSTANTIZE_DUMMY(unsigned int)
INSTANTIZE_DUMMY(unsigned long)
INSTANTIZE_DUMMY(complex<char>)
INSTANTIZE_DUMMY(complex<unsigned char>)
INSTANTIZE_DUMMY(complex<short>)
INSTANTIZE_DUMMY(complex<int>)
INSTANTIZE_DUMMY(complex<unsigned int>)
INSTANTIZE_DUMMY(complex<complex<double> >)
INSTANTIZE_DUMMY(complex<complex<float> >)
INSTANTIZE_DUMMY(complex<unsigned long>)