Commit 3074e5eb authored by thomas.forbriger's avatar thomas.forbriger Committed by thomas.forbriger
Browse files

all the new stuff for LAPACK access in C++

This is a legacy commit from before 2015-03-01.
It may be incomplete as well as inconsistent.
See COPYING.legacy and README.history for details.


SVN Path:     http://gpitrsvn.gpi.uni-karlsruhe.de/repos/TFSoftware/trunk
SVN Revision: 1477
SVN UUID:     67feda4a-a26e-11df-9d6e-31afc202ad0c
parent 819aab5a
This is a legacy version of the repository. It may be incomplete as well as
inconsistent. See README.history for details. For the old stock of the
repository copyright and licence conditions apply as specified for versions
commited after 2015-03-01. Use recent versions as a base for new development.
The legacy version is only stored to keep a record of history.
// LAPACK++ (V. 1.1)
// (C) 1992-1996 All Rights Reserved.
// Linkage names between C, C++, and Fortran (platform dependent)
#ifndef _ARCH_H_
#define _ARCH_H_
#if defined(RIOS) && !defined(CLAPACK)
#define F77NAME(x) x
#else
#include <generic.h>
#define F77NAME(x) name2(x,_)
#endif
#if defined(SGI) && !defined(SGI_DEC)
#define SGI_DEC
extern "C" {
void mkidxname() {}
void mkdatname() {}
}
#endif
#endif // _ARCH_H_
// LAPACK++ (V. 1.1)
// (C) 1992-1996 All Rights Reserved.
#ifndef _LA_BAND_FACT_DOUBLE_H
#define _LA_BAND_FACT_DOUBLE_H
#include "lafnames.h"
#include LA_VECTOR_LONG_INT_H
#include LA_BAND_MAT_DOUBLE_H
#include "lapack.h"
class LaBandFactDouble
{
LaBandMatDouble B_;
LaVectorLongInt pivot_;
int info_;
public:
// constructor
inline LaBandFactDouble();
inline LaBandFactDouble(int,int,int);
inline LaBandFactDouble(LaBandMatDouble &);
inline LaBandFactDouble(LaBandFactDouble &);
inline ~LaBandFactDouble();
// extraction functions for components
inline LaBandMatDouble& B();
inline LaVectorLongInt& pivot();
inline int& info();
// operators
inline LaBandFactDouble& ref(LaBandFactDouble &);
inline LaBandFactDouble& ref(LaBandMatDouble &);
};
// constructor/destructor functions
inline LaBandFactDouble::LaBandFactDouble():B_(),pivot_()
{
#ifdef BandFactDouble_DEBUG
cout << " called LaBandFactDouble::LaBandFactDouble() " << endl;
#endif
info_ = 0;
}
inline LaBandFactDouble::LaBandFactDouble(int N, int kl, int ku)
: B_(N,kl,ku),pivot_(kl+ku+1)
{
#ifdef BandFactDouble_DEBUG
cout << " called LaBandFactDouble::LaBandFactDouble(int,int,int) " << endl;
#endif
info_ = 0;
}
inline LaBandFactDouble::LaBandFactDouble(LaBandMatDouble &G):pivot_()
{
#ifdef BandFactDouble_DEBUG
cout << " called LaBandFactDouble::LaBandFactDouble(LaBandMatDouble &)"
<<endl;
#endif
B_.ref(G);
info_ = 0;
}
inline LaBandFactDouble::LaBandFactDouble(LaBandFactDouble &F)
{
#ifdef BandFactDouble_DEBUG
cout << " called LaBandFactDouble::LaBandFactDouble(LaBandFactDouble &) " << endl;
#endif
B_.ref(F.B_);
pivot_.ref(F.pivot_);
info_ = F.info_;
}
inline LaBandFactDouble::~LaBandFactDouble()
{
}
// member functions
inline LaBandMatDouble& LaBandFactDouble::B()
{
return B_;
}
inline LaVectorLongInt& LaBandFactDouble::pivot()
{
return pivot_;
}
inline int& LaBandFactDouble::info()
{
return info_;
}
// operators
inline LaBandFactDouble& LaBandFactDouble::ref(LaBandFactDouble& F)
{
B_.ref(F.B_);
pivot_.ref(F.pivot_);
info_ = F.info_;
return *this;
}
inline LaBandFactDouble& LaBandFactDouble::ref(LaBandMatDouble& F)
{
B_.ref(F);
return *this;
}
inline void LaBandMatFactorize(LaBandMatDouble &A, LaBandFactDouble &AF)
{
integer n = A.size(1), m = n, LDA = A.gdim(0);
integer KL = A.subdiags(), KU = A.superdiags(), info=0;
F77NAME(dgbtrf)(&m, &n, &KL, &KU, &A(0,0), &LDA, &(AF.pivot()(0)), &info);
}
inline void LaLinearSolve(LaBandFactDouble &AF, LaGenMatDouble &X,
LaGenMatDouble &B)
{
char trans = 'N';
integer n = AF.B().size(1), lda = AF.B().gdim(0), nrhs = X.size(1),
ldb = B.size(0), kl = AF.B().subdiags(),
ku = AF.B().superdiags(), info=0;
X.inject(B);
F77NAME(dgbtrs)(
&trans,
&n,
&kl, &ku, &nrhs, &(AF.B()(-kl,0)),
&lda, &(AF.pivot()(0)), &X(0,0), &ldb, &info);
}
#endif
// LAPACK++ (V. 1.1)
#ifndef _BLAS_PP_H_
#define _BLAS_PP_H_
// requires
//
#include "laexcp.h"
#include "blas1++.h"
#include "blas2++.h"
#include "blas3++.h"
double abs(double);
//-------------------------------------
// Vector/Vector operators
//-------------------------------------
#ifdef _LA_VECTOR_DOUBLE_H_
inline LaVectorDouble operator*(const LaVectorDouble &x, double a)
{
int N = x.size();
LaVectorDouble t(N);
for (int i=0; i<N; i++)
{
t(i) = a * x(i);
}
return t;
}
inline LaVectorDouble operator*(double a, const LaVectorDouble &x)
{
return operator*(x,a);
}
inline double operator*(const LaVectorDouble &dx,
const LaVectorDouble &dy)
{
assert(dx.size()==dy.size());
integer incx = dx.inc(), incy = dy.inc(), n = dx.size();
return F77NAME(ddot)(&n, &dx(0), &incx, &dy(0), &incy);
}
inline LaVectorDouble operator+(const LaVectorDouble &dx,
const LaVectorDouble &dy)
{
assert(dx.size()==dy.size());
integer incx = dx.inc(), incy = dx.inc(), n = dx.size();
double da = 1.0;
LaVectorDouble tmp((int) n);
tmp = dy;
F77NAME(daxpy)(&n, &da, &dx(0), &incx, &tmp(0), &incy);
return tmp;
}
inline LaVectorDouble operator-(const LaVectorDouble &dx,
const LaVectorDouble &dy)
{
assert(dx.size()==dy.size());
integer incx = dx.inc(), incy = dy.inc(), n = dx.size();
double da = -1.0;
LaVectorDouble tmp(n);
tmp = dx;
F77NAME(daxpy)(&n, &da, &dy(0), &incx, &tmp(0), &incy);
return tmp;
}
//-------------------------------------
// Matrix/Vector operators
//-------------------------------------
#ifdef _LA_GEN_MAT_DOUBLE_H_
inline LaVectorDouble operator*(const LaGenMatDouble &A,
const LaVectorDouble &dx)
{
char trans = 'N';
double alpha = 1.0, beta = 0.0;
integer M = A.size(0), N = A.size(1), lda = A.gdim(0);
LaVectorDouble dy(M);
integer incx = dx.inc();
integer incy = dy.inc();
dy = 0.0;
F77NAME(dgemv)(&trans, &M, &N, &alpha, &A(0,0), &lda, &dx(0), &incx,
&beta, &dy(0), &incy);
return dy;
}
#endif
#ifdef _LA_BAND_MAT_DOUBLE_H_
inline LaVectorDouble operator*(const LaBandMatDouble &A,
const LaVectorDouble &dx)
{
char trans = 'N';
double alpha = 1.0, beta = 0.0;
integer M = A.size(0), N = A.size(1), lda = A.gdim(0),
kl = A.subdiags(), ku = A.superdiags();
LaVectorDouble dy(N);
integer incx = dx.inc(), incy = dy.inc();
F77NAME(dgbmv)(&trans, &M, &N, &kl, &ku, &alpha, &A(0,0), &lda,
&dx(0), &incx, &beta, &dy(0), &incy);
return dy;
}
#endif
#ifdef _LA_SYMM_MAT_DOUBLE_H_
inline LaVectorDouble operator*(const LaSymmMatDouble &A,
const LaVectorDouble &dx)
{
char uplo = 'L';
double alpha = 1.0, beta = 0.0;
integer N = A.size(1), lda = A.gdim(0);
LaVectorDouble dy(N);
integer incx = dx.inc(), incy = dy.inc();
F77NAME(dsymv)(&uplo, &N, &alpha, &A(0,0), &lda, &dx(0), &incx,
&beta, &dy(0), &incy);
return dy;
}
#endif
#ifdef _LA_SYMM_BAND_MAT_DOUBLE_H_
inline LaVectorDouble operator*(const LaSymmBandMatDouble &A,
const LaVectorDouble &dx)
{
char uplo = 'L';
double alpha = 1.0, beta = 0.0;
integer N = A.size(1), lda = A.gdim(0), k = A.subdiags();
LaVectorDouble dy(N);
integer incx = dx.inc(), incy = dy.inc();
F77NAME(dsbmv)(&uplo, &N, &k, &alpha, &A(0,0), &lda, &dx(0), &incx,
&beta, &dy(0), &incy);
return dy;
}
#endif
#ifdef _LA_SPD_MAT_DOUBLE_H_
inline LaVectorDouble operator*(const LaSpdMatDouble &AP,
const LaVectorDouble &dx)
{
char uplo = 'L';
double alpha = 1.0, beta = 0.0;
integer N = AP.size(1), incx = dx.inc();
integer lda = AP.gdim(0);
LaVectorDouble dy(N);
integer incy = dy.inc();
F77NAME(dsymv)(&uplo, &N, &alpha, &AP(0,0), &lda, &dx(0), &incx, &beta,
&dy(0), &incy);
return dy;
}
#endif
#ifdef _LA_LOWER_TRIANG_MAT_DOUBLE_H_
inline LaVectorDouble operator*(const LaLowerTriangMatDouble &A,
const LaVectorDouble &dx)
{
char uplo = 'L', trans = 'N', diag = 'N';
integer N = A.size(1), lda = A.gdim(0),
incx = dx.inc();
LaVectorDouble dy(dx);
F77NAME(dtrmv)(&uplo, &trans, &diag, &N, &A(0,0), &lda, &dy(0), &incx);
return dy;
}
#endif
#ifdef _LA_UPPER_TRIANG_MAT_DOUBLE_H_
inline LaVectorDouble operator*(const LaUpperTriangMatDouble &A,
const LaVectorDouble &dx)
{
char uplo = 'U', trans = 'N', diag = 'N';
integer N = A.size(1), lda = A.gdim(0),
incx = dx.inc();
LaVectorDouble dy(dx);
F77NAME(dtrmv)(&uplo, &trans, &diag, &N, &A(0,0), &lda, &dy(0), &incx);
return dy;
}
#endif
#ifdef _LA_UNIT_LOWER_TRIANG_MAT_DOUBLE_H_
inline LaVectorDouble operator*(const LaUnitLowerTriangMatDouble &A,
const LaVectorDouble &dx)
{
char uplo = 'L', trans = 'N', diag = 'U';
integer N = A.size(1), lda = A.gdim(0),
incx = dx.inc();
LaVectorDouble dy(dx);
F77NAME(dtrmv)(&uplo, &trans, &diag, &N, &A(0,0), &lda, &dy(0), &incx);
return dy;
}
#endif
#ifdef _LA_UNIT_UPPER_TRIANG_MAT_DOUBLE_H_
inline LaVectorDouble operator*(const LaUnitUpperTriangMatDouble &A,
const LaVectorDouble &dx)
{
char uplo = 'U', trans = 'N', diag = 'U';
integer N = A.size(1), lda = A.gdim(0),
incx = dx.inc();
LaVectorDouble dy(dx);
F77NAME(dtrmv)(&uplo, &trans, &diag, &N, &A(0,0), &lda, &dy(0), &incx);
return dy;
}
#endif
//-------------------------------------
// Matrix/Matrix operators
//-------------------------------------
inline LaGenMatDouble operator*(const LaGenMatDouble &A,
const LaGenMatDouble &B)
{
char t = 'N';
integer m = A.size(0), k = A.size(1), n = B.size(1);
integer lda = A.gdim(0), ldb = B.gdim(0);
double alpha = 1.0, beta = 1.0;
LaGenMatDouble C(m,n);
integer ldc = A.gdim(0);
C = 0.0;
F77NAME(dgemm)(&t, &t, &m, &n, &k, &alpha, &A(0,0), &lda, &B(0,0), &ldb,
&beta, &C(0,0), &ldc);
return C;
}
#ifdef _LA_UNIT_LOWER_TRIANG_MAT_DOUBLE_H_
inline LaGenMatDouble operator*(const LaUnitLowerTriangMatDouble &A,
const LaGenMatDouble &B)
{
char side = 'L', uplo = 'L', transa = 'N', diag = 'U';
double alpha = 1.0;
integer m = B.size(0), n = B.size(1),
lda = A.gdim(0), ldb = B.gdim(0);
LaGenMatDouble C(B);
F77NAME(dtrmm)(&side, &uplo, &transa, &diag, &m, &n, &alpha,
&A(0,0), &lda, &C(0,0), &ldb);
return C;
}
#endif
#ifdef _LA_UNIT_UPPER_TRIANG_MAT_DOUBLE_H_
inline LaGenMatDouble operator*(const LaUnitUpperTriangMatDouble &A,
const LaGenMatDouble &B)
{
char side = 'L', uplo = 'U', transa = 'N', diag = 'U';
double alpha = 1.0;
integer m = B.size(0), n = B.size(1),
lda = A.gdim(0), ldb = B.gdim(0);
LaGenMatDouble C(B);
F77NAME(dtrmm)(&side, &uplo, &transa, &diag, &m, &n, &alpha,
&A(0,0), &lda, &C(0,0), &ldb);
return C;
}
#endif
#ifdef _LA_LOWER_TRIANG_MAT_DOUBLE_H_
inline LaGenMatDouble operator*(const LaLowerTriangMatDouble &A,
const LaGenMatDouble &B)
{
char side = 'L', uplo = 'L', transa = 'N', diag = 'N';
double alpha = 1.0;
integer m = B.size(0), n = B.size(1),
lda = A.gdim(0), ldb = B.gdim(0);
LaGenMatDouble C(B);
F77NAME(dtrmm)(&side, &uplo, &transa, &diag, &m, &n, &alpha,
&A(0,0), &lda, &C(0,0), &ldb);
return C;
}
#endif
#ifdef _LA_UPPER_TRIANG_MAT_DOUBLE_H_
inline LaGenMatDouble operator*(const LaUpperTriangMatDouble &A,
const LaGenMatDouble &B)
{
char side = 'L', uplo = 'U', transa = 'N', diag = 'N';
double alpha = 1.0;
integer m = B.size(0), n = B.size(1),
lda = A.gdim(0), ldb = B.gdim(0);
LaGenMatDouble C(B);
F77NAME(dtrmm)(&side, &uplo, &transa, &diag, &m, &n, &alpha,
&A(0,0), &lda, &C(0,0), &ldb);
return C;
}
#endif
#ifdef _LA_SYMM_MAT_DOUBLE_H_
inline LaGenMatDouble operator*(const LaSymmMatDouble &A,
const LaGenMatDouble &B)
{
char side = 'L', uplo = 'L';
double alpha = 1.0, beta = 1.0;
LaGenMatDouble C(B.size(1),A.size(1));
integer m = C.size(0), n = C.size(1), lda = A.gdim(0),
ldb = B.gdim(0), ldc = C.gdim(0);
F77NAME(dsymm)(&side, &uplo, &m, &n, &alpha, &A(0,0), &lda,
&B(0,0), &ldb, &beta, &C(0,0), &ldc);
return C;
}
#endif
#ifdef _LA_SYMM_TRIDIAG_MAT_DOUBLE_H_
inline LaVectorDouble operator*(const LaSymmTridiagMatDouble& A,
const LaVectorDouble& X)
{
integer M = A.size();
integer N = X.size();
LaVectorDouble R(M);
R(0) = ((A.diag(0)(0) * X(0)) + (A.diag(-1)(0) * X(1)));
for (integer i = 1; i < M-2; i++)
{
R(i) = ((A.diag(-1)(i-1) * X(i-1)) +
(A.diag(0)(i) * X(i)) +
(A.diag(-1)(i) * X(i+1)));
}
R(M-1) = ((A.diag(0)(M-1) * X(N-1)) + (A.diag(-1)(M-2) * X(N-2)));
return R;
}
#endif
#ifdef _LA_TRIDIAG_MAT_DOUBLE_H_
inline LaVectorDouble operator*(const LaTridiagMatDouble& A,
const LaVectorDouble& X)
{
integer M = A.size();
integer N = X.size();
LaVectorDouble R(M);
R(0) = ((A.diag(0)(0) * X(0)) + (A.diag(-1)(0) * X(1)));
for (integer i = 1; i < M-2; i++)
{
R(i) = ((A.diag(-1)(i-1) * X(i-1)) +
(A.diag(0)(i) * X(i)) +
(A.diag(1)(i) * X(i+1)));
}
R(M-1) = ((A.diag(0)(M-1) * X(N-1)) + (A.diag(1)(M-2) * X(N-2)));
return R;
}
#endif
inline LaGenMatDouble operator-(const LaGenMatDouble &A,
const LaGenMatDouble &B)
{
#ifndef HPPA
const char fname[] = "operator+(A,B)";
#else
char *fname = NULL;
#endif