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

*** empty log message ***

This commit is contained in:
jiri 2006-08-16 21:43:45 +00:00
parent 24c048d210
commit 9b69ed529f
9 changed files with 58 additions and 71 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

13
fourindex.h
View File

@ -619,6 +619,11 @@ fourindex_dense<twoelectronrealmullikan,T,I>::fourindex_dense<twoelectronrealmul
{
if(rhs.getsymmetry() != twoelectronrealmullikan ) laerror("fourindex_dense symmetry mismatch");
typename fourindex<I,T>::iterator p;
#ifdef DEBUG
unsigned long I = SMat_index_1(p->index.indiv.i,p->index.indiv.j);
unsigned long J = SMat_index_1(p->index.indiv.k,p->index.indiv.l);
if (I<0 || I>=(unsigned long)NRSMat<T>::nn || J<0 || J>=(unsigned long)NRSMat<T>::nn) laerror("fourindex_dense index out of range in constructor");
#endif
for(p=rhs.begin(); p!= rhs.end(); ++p) (*this)(p->index.indiv.i,p->index.indiv.j,p->index.indiv.k,p->index.indiv.l) = p->elem;
}
@ -627,8 +632,12 @@ fourindex_dense<twoelectronrealmullikan,T,I>::fourindex_dense<twoelectronrealmul
{
if(rhs.getsymmetry() != twoelectronrealmullikan ) laerror("fourindex_dense symmetry mismatch");
typename fourindex_ext<I,T>::iterator p;
for(p=rhs.begin(); p!= rhs.end(); ++p) (*this)(p->index.indiv.i,p->index.indiv.j
,p->index.indiv.k,p->index.indiv.l) = p->elem;
#ifdef DEBUG
unsigned long I = SMat_index_1(p->index.indiv.i,p->index.indiv.j);
unsigned long J = SMat_index_1(p->index.indiv.k,p->index.indiv.l);
if (I<0 || I>=(unsigned long)NRSMat<T>::nn || J<0 || J>=(unsigned long)NRSMat<T>::nn) laerror("fourindex_dense index out of range in constructor");
#endif
for(p=rhs.begin(); p!= rhs.end(); ++p) (*this)(p->index.indiv.i,p->index.indiv.j ,p->index.indiv.k,p->index.indiv.l) = p->elem;
}

3
la_traits.h
View File

@ -21,9 +21,6 @@ extern "C" {
}
#endif
#ifdef _GLIBCPP_NO_TEMPLATE_EXPORT
# define export
#endif
//forward declarations
template<typename C> class NRVec;

31
laerror.cc
View File

@ -3,6 +3,8 @@
#include "laerror.h"
#include <stdio.h>
#include <errno.h>
#include <stdarg.h>
#ifdef USE_TRACEBACK
#include "traceback.h"
@ -30,7 +32,7 @@ void laerror(const char *s1)
//stub for f77 blas called from strassen routine
extern "C" void xerbla_(const char name[6], int *n)
{
char msg[128];
char msg[1024];
strcpy(msg,"LAPACK or BLAS error in routine ");
strncat(msg,name,6);
sprintf(msg+strlen(msg),": illegal value of parameter #%d",*n);
@ -38,3 +40,30 @@ laerror(msg);
}
//with atlas-cblas another error routine is necessary
extern "C" void ATL_xerbla(int p, char *rout, char *form, ...)
{
char msg0[1024], *msg;
va_list argptr;
va_start(argptr, form);
strcpy(msg0,"ATLAS error\n");
msg=msg0+strlen(msg0);
if (p) {sprintf(msg, "Parameter %d to routine %s was incorrect\n", p, rout); msg+=strlen(msg);}
vsprintf(msg, form, argptr);
va_end(argptr);
laerror(msg0);
}
int cblas_errprn(int ierr, int info, char *form, ...)
{
char msg0[1024], *msg;
va_list argptr;
va_start(argptr, form);
sprintf(msg0,"CBLAS error %d %d\n",ierr,info);
msg=msg0+strlen(msg0);
vsprintf(msg, form, argptr);
va_end(argptr);
laerror(msg0);
return 0;
}

32
mat.cc
View File

@ -779,8 +779,8 @@ void NRMat< complex<double> >::diagmultr(const NRVec< complex<double> > &rhs)
#ifdef oldversion
// Mat * Smat, decomposed to nn x Vec * Smat
//NOTE: dsymm is not appropriate as it works on UNPACKED symmetric matrix
template<>
const NRMat<double>
NRMat<double>::operator*(const NRSMat<double> &rhs) const
@ -811,36 +811,6 @@ NRMat< complex<double> >::operator*(const NRSMat< complex<double> > &rhs) const
return result;
}
#else
// Mat * Smat
template<>
const NRMat<double>
NRMat<double>::operator*(const NRSMat<double> &rhs) const
{
#ifdef DEBUG
if (mm != rhs.nrows()) laerror("incompatible dimension in Mat*SMat");
#endif
NRMat<double> result(nn, rhs.ncols());
cblas_dsymm(CblasRowMajor, CblasRight, CblasLower, nn, mm, 1., &rhs[0],mm,(*this)[0],mm,0.,result[0],mm);
return result;
}
template<>
const NRMat< complex<double> >
NRMat< complex<double> >::operator*(const NRSMat< complex<double> > &rhs) const
{
#ifdef DEBUG
if (mm != rhs.nrows()) laerror("incompatible dimension in Mat*SMat");
#endif
NRMat< complex<double> > result(nn, rhs.ncols());
cblas_zhemm(CblasRowMajor, CblasRight, CblasLower, nn, mm, (void *)&CONE, &rhs[0],mm,(*this)[0],mm,(void *)&CZERO,result[0],mm);
return result;
}
#endif
// sum of rows

1
noncblas.cc
View File

@ -222,4 +222,5 @@ for(int i=1; i<N; ++i) for(int j=0; j<i; ++j) {double t=A[j*lda+i]; A[j*lda+i]=
return INFO;
}
#endif

2
noncblas.h
View File

@ -16,7 +16,7 @@
#define CBLAS_INDEX int
int cblas_errprn(int ierr, int info, char *form, ...);
void ATL_xerbla(int p, char *rout, char *form, ...);
/*
* ===========================================================================
* Prototypes for level 1 BLAS functions (complex are recast as routines)

34
smat.cc
View File

@ -61,6 +61,7 @@ LA_traits<T>::multiget(NN2,fd,v,dim);
template <typename T>
NRSMat<T>::NRSMat(const NRMat<T> &rhs)
{
nn=rhs.nrows();
#ifdef DEBUG
if (nn != rhs.ncols()) laerror("attempt to convert non-square Mat to SMat");
#endif
@ -148,8 +149,10 @@ void NRSMat<T>::fscanf(FILE *f, const char *format)
* BLAS specializations for double and complex<double>
*/
#ifdef oldversion
// SMat * Mat
//NOTE: dsymm is not appropriate as it works on UNPACKED symmetric matrix
template<>
const NRMat<double> NRSMat<double>::operator*(const NRMat<double> &rhs) const
{
@ -178,35 +181,6 @@ NRSMat< complex<double> >::operator*(const NRMat< complex<double> > &rhs) const
return result;
}
#else
// SMat * Mat
template<>
const NRMat<double> NRSMat<double>::operator*(const NRMat<double> &rhs) const
{
#ifdef DEBUG
if (nn != rhs.nrows()) laerror("incompatible dimensions in SMat*Mat");
#endif
NRMat<double> result(nn, rhs.ncols());
cblas_dsymm(CblasRowMajor, CblasLeft, CblasLower, nn, rhs.ncols(), 1., (*this),nn, rhs[0],rhs.ncols(), 0.,result[0],rhs.ncols());
return result;
}
template<>
const NRMat< complex<double> >
NRSMat< complex<double> >::operator*(const NRMat< complex<double> > &rhs) const
{
#ifdef DEBUG
if (nn != rhs.nrows()) laerror("incompatible dimensions in SMat*Mat");
#endif
NRMat< complex<double> > result(nn, rhs.ncols());
cblas_zhemm(CblasRowMajor, CblasLeft, CblasLower, nn, rhs.ncols(), &CONE, (*this),nn, rhs[0],rhs.ncols(), &CZERO,result[0],rhs.ncols());
return result;
}
#endif
// SMat * SMat

2
sparsemat.cc
View File

@ -1309,7 +1309,6 @@ return *this;
}
#ifdef _GLIBCPP_NO_TEMPLATE_EXPORT
#define INSTANTIZE(T) \
template SparseMat<T> & SparseMat<T>::oplusequal(const SparseMat<T> &rhs);\
@ -1363,4 +1362,3 @@ INSTANTIZE(double)
INSTANTIZE(complex<double>) //some functions are not OK for hermitean matrices, needs a revision!!!
#endif

11
t.cc
View File

@ -1105,7 +1105,7 @@ for(int iter=1; iter<100 && norm>1e-8 ; ++iter)
}
}
if(1)
if(0)
{
NRMat<double> a,b;
cin >>a;
@ -1115,5 +1115,14 @@ cout <<b*b;
cout <<(b*b-a).norm();
}
if(1)
{
NRSMat<double> a;
NRMat<double> b;
cin >>a>>b;
cout <<a*b;
}
}