LA_library/nonclass.h

396 lines
14 KiB
C
Raw Normal View History

2008-02-26 14:55:23 +01:00
/*
LA: linear algebra C++ interface library
Copyright (C) 2008 Jiri Pittner <jiri.pittner@jh-inst.cas.cz> or <jiri@pittnerovi.com>
This program is free software: you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
You should have received a copy of the GNU General Public License
along with this program. If not, see <http://www.gnu.org/licenses/>.
*/
2005-02-01 00:08:03 +01:00
#ifndef _LA_NONCLASS_H_
#define _LA_NONCLASS_H_
2004-03-17 04:07:21 +01:00
#include "vec.h"
#include "smat.h"
#include "mat.h"
2005-02-06 15:01:27 +01:00
#include "la_traits.h"
2004-03-17 04:07:21 +01:00
2009-11-12 22:01:19 +01:00
namespace LA {
2005-02-04 15:31:42 +01:00
2004-03-17 04:07:21 +01:00
//MISC
2009-10-08 16:01:15 +02:00
template <class T>
2006-08-15 22:29:01 +02:00
const NRSMat<T> twoside_transform(const NRSMat<T> &S, const NRMat<T> &C, bool transp=0) //calculate C^dagger S C
2006-08-15 22:23:32 +02:00
{
2006-08-15 22:29:01 +02:00
if(transp)
{
NRMat<T> tmp = C * S;
2006-08-15 22:29:55 +02:00
NRMat<T> result(C.nrows(),C.nrows());
2006-08-15 22:29:01 +02:00
result.gemm((T)0,tmp,'n',C,'t',(T)1);
return NRSMat<T>(result);
}
2006-08-15 22:23:32 +02:00
NRMat<T> tmp = S * C;
NRMat<T> result(C.ncols(),C.ncols());
result.gemm((T)0,C,'t',tmp,'n',(T)1);
return NRSMat<T>(result);
}
2009-10-08 16:01:15 +02:00
template <class T>
2004-03-17 17:39:07 +01:00
const NRMat<T> diagonalmatrix(const NRVec<T> &x)
{
int n=x.size();
NRMat<T> result((T)0,n,n);
T *p = result[0];
for(int j=0; j<n; j++) {*p = x[j]; p+=(n+1);}
return result;
}
2004-03-17 04:07:21 +01:00
//more efficient commutator for a special case of full matrices
template<class T>
inline const NRMat<T> commutator ( const NRMat<T> &x, const NRMat<T> &y, const bool trx=0, const bool tryy=0)
{
NRMat<T> r(trx?x.ncols():x.nrows(), tryy?y.nrows():y.ncols());
r.gemm((T)0,x,trx?'t':'n',y,tryy?'t':'n',(T)1);
r.gemm((T)1,y,tryy?'t':'n',x,trx?'t':'n',(T)-1);
return r;
}
//more efficient commutator for a special case of full matrices
template<class T>
inline const NRMat<T> anticommutator ( const NRMat<T> &x, const NRMat<T> &y, const bool trx=0, const bool tryy=0)
{
NRMat<T> r(trx?x.ncols():x.nrows(), tryy?y.nrows():y.ncols());
r.gemm((T)0,x,trx?'t':'n',y,tryy?'t':'n',(T)1);
r.gemm((T)1,y,tryy?'t':'n',x,trx?'t':'n',(T)1);
return r;
}
//////////////////////
// LAPACK interface //
//////////////////////
#define declare_la(T) \
extern const NRVec<T> diagofproduct(const NRMat<T> &a, const NRMat<T> &b,\
bool trb=0, bool conjb=0); \
extern T trace2(const NRMat<T> &a, const NRMat<T> &b, bool trb=0); \
extern T trace2(const NRSMat<T> &a, const NRSMat<T> &b, const bool diagscaled=0);\
2006-09-10 22:06:44 +02:00
extern T trace2(const NRSMat<T> &a, const NRMat<T> &b, const bool diagscaled=0);\
2019-11-13 23:22:25 +01:00
extern T trace2(const NRMat<T> &a, const NRSMat<T> &b, const bool diagscaled=0);\
2011-01-18 15:37:05 +01:00
extern void linear_solve(NRMat<T> &a, NRMat<T> *b, T *det=0,int n=0); /*solve Ax^T=b^T (b is nrhs x n) */ \
extern void linear_solve(NRSMat<T> &a, NRMat<T> *b, T *det=0, int n=0); /*solve Ax^T=b^T (b is nrhs x n) */\
2005-02-17 23:54:27 +01:00
extern void linear_solve(NRMat<T> &a, NRVec<T> &b, double *det=0, int n=0); \
extern void linear_solve(NRSMat<T> &a, NRVec<T> &b, double *det=0, int n=0); \
2009-05-28 15:02:01 +02:00
extern void diagonalize(NRMat<T> &a, NRVec<LA_traits<T>::normtype> &w, const bool eivec=1, const bool corder=1, int n=0, NRMat<T> *b=NULL, const int itype=1); \
extern void diagonalize(NRSMat<T> &a, NRVec<LA_traits<T>::normtype> &w, NRMat<T> *v, const bool corder=1, int n=0, NRSMat<T> *b=NULL, const int itype=1);\
2019-11-13 23:22:25 +01:00
extern void singular_decomposition(NRMat<T> &a, NRMat<T> *u, NRVec<LA_traits<T>::normtype> &s, NRMat<T> *v, const bool vnotdagger=0, int m=0, int n=0);
2004-03-17 04:07:21 +01:00
2006-05-28 16:12:01 +02:00
/*NOTE!!! all versions of diagonalize DESTROY A and generalized diagonalize also B matrix */
2004-03-17 04:07:21 +01:00
declare_la(double)
2020-01-06 21:50:34 +01:00
declare_la(std::complex<double>)
2004-03-17 04:07:21 +01:00
// Separate declarations
2005-02-17 23:54:27 +01:00
//general nonsymmetric matrix and generalized diagonalization
2011-01-18 15:37:05 +01:00
//corder =0 ... C rows are eigenvectors, =1 ... C columns are eigenvectors
2004-03-17 04:07:21 +01:00
extern void gdiagonalize(NRMat<double> &a, NRVec<double> &wr, NRVec<double> &wi,
NRMat<double> *vl=NULL, NRMat<double> *vr=NULL, const bool corder=1, int n=0, const int sorttype=0, const int biorthonormalize=0,
2011-01-18 15:37:05 +01:00
NRMat<double> *b=NULL, NRVec<double> *beta=NULL); //this used real storage of eigenvectors like dgeev
template<typename T>
2020-01-06 21:50:34 +01:00
extern void gdiagonalize(NRMat<T> &a, NRVec< std::complex<double> > &w,
NRMat< std::complex<double> >*vl=NULL, NRMat< std::complex<double> > *vr=NULL,
2008-03-02 17:40:22 +01:00
const bool corder=1, int n=0, const int sorttype=0, const int biorthonormalize=0,
2011-01-18 15:37:05 +01:00
NRMat<T> *b=NULL, NRVec<T> *beta=NULL); //eigenvectors are stored in complex matrices for T both double and complex
2004-03-17 04:07:21 +01:00
//for compatibility in davidson
extern void gdiagonalize(NRMat<std::complex<double> > &a, NRVec<double> &wr, NRVec<double> &wi,
NRMat<std::complex<double> > *vl=NULL, NRMat<std::complex<double> > *vr=NULL, const bool corder=1, int n=0, const int sorttype=0, const int biorthonormalize=0,
NRMat<std::complex<double> > *b=NULL, NRVec<std::complex<double> > *beta=NULL);
2011-01-18 15:37:05 +01:00
//complex,real,imaginary parts of various entities
template<typename T>
extern const typename LA_traits<T>::realtype realpart(const T&);
template<typename T>
extern const typename LA_traits<T>::realtype imagpart(const T&);
template<typename T>
extern const typename LA_traits<T>::complextype realmatrix (const T&);
template<typename T>
extern const typename LA_traits<T>::complextype imagmatrix (const T&);
template<typename T>
extern const typename LA_traits<T>::complextype complexmatrix (const T&, const T&);
2004-03-17 04:07:21 +01:00
2010-01-07 17:10:12 +01:00
//Cholesky decomposition
extern void cholesky(NRMat<double> &a, bool upper=1);
2020-01-06 21:50:34 +01:00
extern void cholesky(NRMat<std::complex<double> > &a, bool upper=1);
2010-01-07 17:10:12 +01:00
2004-03-17 04:07:21 +01:00
//inverse by means of linear solve, preserving rhs intact
template<typename T>
const NRMat<T> inverse(NRMat<T> a, T *det=0)
{
#ifdef DEBUG
if(a.nrows()!=a.ncols()) laerror("inverse() for non-square matrix");
#endif
NRMat<T> result(a.nrows(),a.nrows());
result = (T)1.;
2008-03-02 17:40:22 +01:00
a.copyonwrite();
2004-03-17 04:07:21 +01:00
linear_solve(a, &result, det);
2008-03-02 17:40:22 +01:00
result.transposeme(); //tested with noncblas
2004-03-17 04:07:21 +01:00
return result;
}
2010-01-17 21:28:38 +01:00
//several matrix norms
template<class MAT>
typename LA_traits<MAT>::normtype MatrixNorm(const MAT &A, const char norm);
//condition number
template<class MAT>
typename LA_traits<MAT>::normtype CondNumber(const MAT &A, const char norm);
2005-02-06 15:01:27 +01:00
//general determinant
template<class MAT>
2005-09-04 21:34:10 +02:00
const typename LA_traits<MAT>::elementtype determinant(MAT a)//passed by value
{
typename LA_traits<MAT>::elementtype det;
if(a.nrows()!=a.ncols()) laerror("determinant of non-square matrix");
linear_solve(a,NULL,&det);
return det;
}
//general determinant destructive on input
template<class MAT>
const typename LA_traits<MAT>::elementtype determinant_destroy(MAT &a) //passed by reference
2005-02-06 15:01:27 +01:00
{
2005-02-14 01:10:07 +01:00
typename LA_traits<MAT>::elementtype det;
2005-02-06 15:01:27 +01:00
if(a.nrows()!=a.ncols()) laerror("determinant of non-square matrix");
linear_solve(a,NULL,&det);
return det;
}
2005-09-11 22:04:24 +02:00
2010-02-25 21:47:01 +01:00
//------------------------------------------------------------------------------
// solves set of linear equations using gesvx
// input:
// A double precision matrix of dimension nn x mm, where min(nn, mm) >= n
// B double prec. array dimensioned as nrhs x n
// rhsCount nrhs - count of right hand sides
// eqCount n - count of equations
// eq use equilibration of matrix A before solving
// saveA if set, do no overwrite A if equilibration in effect
// rcond if not NULL, store the returned rcond value from dgesvx
// output:
// solution is stored in B
// the info parameter of gesvx is returned (see man dgesvx)
//------------------------------------------------------------------------------
2010-01-17 21:28:38 +01:00
template<class T>
2010-02-25 21:47:01 +01:00
int linear_solve_x(NRMat<T> &A, T *B, const int rhsCount, const int eqCount, const bool eq, const bool saveA, double *rcond);
//------------------------------------------------------------------------------
// for given square matrices A, B computes X = AB^{-1} as follows
// XB = A => B^TX^T = A^T
// input:
// _A double precision matrix of dimension nn x nn
// _B double prec. matrix of dimension nn x nn
// _useEq use equilibration suitable for badly conditioned matrices
// _rcond if not NULL, store the returned value of rcond fromd dgesvx
// output:
// solution is stored in _B
// the info parameter of dgesvx is returned (see man dgesvx)
//------------------------------------------------------------------------------
2010-01-17 21:28:38 +01:00
template<class T>
2020-01-06 21:50:34 +01:00
int multiply_by_inverse(NRMat<T> &A, NRMat<T> &B, bool useEq=false, double *rcond=NULL);
2010-01-17 21:28:38 +01:00
2005-09-05 16:19:44 +02:00
//general submatrix, INDEX will typically be NRVec<int> or even int*
//NOTE: in order to check consistency between nrows and rows in rows is a NRVec
//some advanced metaprogramming would be necessary
2005-09-06 17:55:07 +02:00
//application: e.g. ignoresign=true, equalsigns=true, indexshift= -1 ... elements of Slater overlaps for RHF
2005-09-05 16:19:44 +02:00
2005-09-04 21:34:10 +02:00
template<class MAT, class INDEX>
2005-09-06 17:55:07 +02:00
const NRMat<typename LA_traits<MAT>::elementtype> submatrix(const MAT a, const int nrows, const INDEX rows, const int ncols, const INDEX cols, int indexshift=0, bool ignoresign=false, bool equalsigns=false)
2005-09-04 21:34:10 +02:00
{
2005-09-05 16:19:44 +02:00
NRMat<typename LA_traits<MAT>::elementtype> r(nrows,ncols);
2005-09-04 21:34:10 +02:00
2005-09-06 17:55:07 +02:00
if(equalsigns) //make the element zero if signs of both indices are opposite
{
if(ignoresign)
{
for(int i=0; i<nrows; ++i)
for(int j=0; j<ncols; ++j)
2009-11-12 22:01:19 +01:00
r(i,j) = rows[i]*cols[j]<0?0.:a(std::abs(rows[i])+indexshift,std::abs(cols[j])+indexshift);
2005-09-06 17:55:07 +02:00
}
else
{
for(int i=0; i<nrows; ++i)
for(int j=0; j<ncols; ++j)
r(i,j) = rows[i]*cols[j]<0?0.:a(rows[i]+indexshift,cols[j]+indexshift);
}
}
else
{
2005-09-04 21:34:10 +02:00
if(ignoresign)
{
2005-09-05 16:19:44 +02:00
for(int i=0; i<nrows; ++i)
for(int j=0; j<ncols; ++j)
2009-11-12 22:01:19 +01:00
r(i,j) = a(std::abs(rows[i])+indexshift,std::abs(cols[j])+indexshift);
2005-09-04 21:34:10 +02:00
}
else
{
2005-09-05 16:19:44 +02:00
for(int i=0; i<nrows; ++i)
for(int j=0; j<ncols; ++j)
2005-09-04 21:34:10 +02:00
r(i,j) = a(rows[i]+indexshift,cols[j]+indexshift);
}
2005-09-06 17:55:07 +02:00
}
2005-09-04 21:34:10 +02:00
return r;
}
2005-02-06 15:01:27 +01:00
//auxiliary routine to adjust eigenvectors to guarantee real logarithm
extern void adjustphases(NRMat<double> &v);
2006-04-01 06:48:01 +02:00
2006-04-07 22:47:38 +02:00
//declaration of template interface to cblas routines with full options available
//just to facilitate easy change between float, double, complex in a user code
//very incomplete, add new ones as needed
template<class T> inline void xcopy(int n, const T *x, int incx, T *y, int incy);
template<class T> inline void xaxpy(int n, const T &a, const T *x, int incx, T *y, int incy);
template<class T> inline T xdot(int n, const T *x, int incx, const T *y, int incy);
//specialized definitions have to be in the header file to be inlineable, eliminating any runtime overhead
template<>
inline void xcopy<double> (int n, const double *x, int incx, double *y, int incy)
{
cblas_dcopy(n, x, incx, y, incy);
}
template<>
inline void xaxpy<double>(int n, const double &a, const double *x, int incx, double *y, int incy)
{
cblas_daxpy(n, a, x, incx, y, incy);
}
template<>
inline double xdot<double>(int n, const double *x, int incx, const double *y, int incy)
{
return cblas_ddot(n,x,incx,y,incy);
}
2006-04-01 06:48:01 +02:00
2006-09-12 22:27:41 +02:00
//debugging aid: reconstruct an explicit matrix from the implicit version
//which provides gemv only
template<typename M, typename T>
NRMat<T> reconstructmatrix(const M &implicitmat)
{
NRMat<T> r(implicitmat.nrows(),implicitmat.ncols());
NRVec<T> rhs(0.,implicitmat.ncols());
NRVec<T> tmp(implicitmat.nrows());
for(int i=0; i<implicitmat.ncols(); ++i)
{
rhs[i]=1.;
implicitmat.gemv(0.,tmp,'n',1.,rhs);
for(int j=0; j<implicitmat.nrows(); ++j) r(j,i)=tmp[j];
rhs[i]=0.;
}
return r;
}
2011-01-18 15:37:05 +01:00
//matrix functions via diagonalization
extern NRMat<double> realmatrixfunction(NRMat<double> a, double (*f) (double)); //a has to by in fact symmetric
2020-01-06 21:50:34 +01:00
extern NRMat<std::complex<double> > complexmatrixfunction(NRMat<double> a, double (*fre) (double), double (*fim) (double)); //a has to by in fact symmetric
2011-01-18 15:37:05 +01:00
template<typename T>
NRMat<T> matrixfunction(NRSMat<T> a, double (*f) (double)) //of symmetric/hermitian matrix
{
int n = a.nrows();
NRVec<double> w(n);
NRMat<T> v(n, n);
diagonalize(a, w, &v, 0);
for (int i=0; i<a.nrows(); i++) w[i] = (*f)(w[i]);
NRMat<T> u = v;
NRVec<T> ww=w; //diagmultl needs same type
v.diagmultl(ww);
NRMat<T> r(n, n);
r.gemm(0.0, u, 't', v, 'n', 1.0); //gemm will use 'c' for complex ones
return r;
}
template<typename T>
2020-01-06 21:50:34 +01:00
extern NRMat<T> matrixfunction(NRMat<T> a, std::complex<double> (*f)(const std::complex<double> &)) //of a general real/complex matrix
2011-01-18 15:37:05 +01:00
{
int n = a.nrows();
2020-01-06 21:50:34 +01:00
NRVec<std::complex<double> > w(n);
NRMat<std::complex<double> > u(n,n),v(n,n);
2011-01-18 15:37:05 +01:00
#ifdef debugmf
2020-01-06 21:50:34 +01:00
NRMat<std::complex<double> > a0=a;
2011-01-18 15:37:05 +01:00
#endif
gdiagonalize<T>(a, w, &u, &v, false,n,0,false,NULL,NULL);//a gets destroyed, eigenvectors are rows
2020-01-06 21:50:34 +01:00
NRVec< std::complex<double> > z = diagofproduct(u, v, 1, 1);
2011-01-18 15:37:05 +01:00
#ifdef debugmf
std::cout <<"TEST matrixfunction\n"<<w<<u<<v<<z;
std::cout <<"TEST matrixfunction1 "<< u*a0 - diagonalmatrix(w)*u<<std::endl;
std::cout <<"TEST matrixfunction2 "<< a0*v.transpose(1) - v.transpose(1)*diagonalmatrix(w)<<std::endl;
std::cout <<"TEST matrixfunction3 "<< u*v.transpose(1)<<diagonalmatrix(z)<<std::endl;
#endif
2020-01-06 21:50:34 +01:00
NRVec< std::complex<double> > wz(n);
2011-01-18 15:37:05 +01:00
for (int i=0; i<a.nrows(); i++) wz[i] = w[i]/z[i];
#ifdef debugmf
std::cout <<"TEST matrixfunction4 "<< a0<< v.transpose(true)*diagonalmatrix(wz)*u<<std::endl;
#endif
for (int i=0; i<a.nrows(); i++) w[i] = (*f)(w[i])/z[i];
u.diagmultl(w);
2020-01-06 21:50:34 +01:00
NRMat< std::complex<double> > r(n, n);
2011-01-18 15:37:05 +01:00
r.gemm(0.0, v, 'c', u, 'n', 1.0);
return (NRMat<T>) r; //convert back to real if applicable by the explicit decomplexifying constructor; it is NOT checked to which accuracy the imaginary part is actually zero
}
2020-01-06 21:50:34 +01:00
extern std::complex<double> sqrtinv(const std::complex<double> &);
2011-01-18 15:37:05 +01:00
extern double sqrtinv(const double);
//functions on matrices
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> 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> log(const NRSMat<double> &a) { return matrixfunction(a,&std::log); }
extern NRMat<double> log(const NRMat<double> &a);
2020-01-06 21:50:34 +01:00
extern NRMat<std::complex<double> > log(const NRMat<std::complex<double> > &a);
extern NRMat<std::complex<double> > exp0(const NRMat<std::complex<double> > &a);
extern NRMat<std::complex<double> > copytest(const NRMat<std::complex<double> > &a);
2011-01-18 15:37:05 +01:00
extern NRMat<double> copytest(const NRMat<double> &a);
extern NRMat<double> exp0(const NRMat<double> &a);
2009-11-12 22:01:19 +01:00
}//namespace
2005-02-01 00:08:03 +01:00
#endif