Added: Header files BLAS.h and LAPack.h with platform-dependent implementations.
Changed: Replaced potentially confusing macro symbol USE_BLAS with HAS_BLAS which is an enum value depending on which BLAS implementation (if any) is used. Added: Wrappers for the LAPack/BLAS functions when HAS_BLAS == 3 (Accelerate).
This commit is contained in:
Knut Morten Okstad
31
src/LinAlg/BLAS.h
Normal file
31
src/LinAlg/BLAS.h
Normal file
@@ -0,0 +1,31 @@
|
||||
// $Id$
|
||||
//==============================================================================
|
||||
//!
|
||||
//! \file BLAS.h
|
||||
//!
|
||||
//! \date Jan 10 2016
|
||||
//!
|
||||
//! \author Knut Morten Okstad / SINTEF
|
||||
//!
|
||||
//! \brief BLAS support for various platforms.
|
||||
//!
|
||||
//==============================================================================
|
||||
|
||||
#ifndef UTL_BLAS_H
|
||||
#define UTL_BLAS_H
|
||||
|
||||
#if defined(USE_MKL)
|
||||
#define HAS_BLAS 2
|
||||
#include <mkl_cblas.h>
|
||||
#elif defined(USE_ACCELERATE)
|
||||
#define HAS_BLAS 3
|
||||
#include <Accelerate/Accelerate.h>
|
||||
#elif defined(USE_CBLAS)
|
||||
#define HAS_BLAS 1
|
||||
extern "C"
|
||||
{
|
||||
#include <cblas.h>
|
||||
}
|
||||
#endif
|
||||
|
||||
#endif
|
||||
@@ -14,128 +14,47 @@
|
||||
#include "DenseMatrix.h"
|
||||
#include "SparseMatrix.h"
|
||||
#include "SAM.h"
|
||||
#include "LAPack.h"
|
||||
|
||||
#ifdef USE_F77SAM
|
||||
#if defined(_WIN32) && !defined(__MINGW32__) && !defined(__MINGW64__)
|
||||
#define addem2_ ADDEM2
|
||||
#define dgecon_ DGECON
|
||||
#define dgesv_ DGESV
|
||||
#define dgetrs_ DGETRS
|
||||
#define dlange_ DLANGE
|
||||
#define dposv_ DPOSV
|
||||
#define dpotrs_ DPOTRS
|
||||
#define dsyevx_ DSYEVX
|
||||
#define dsygvx_ DSYGVX
|
||||
#define dgeev_ DGEEV
|
||||
#elif defined(_AIX)
|
||||
#define addem2_ addem2
|
||||
#define dgecon_ dgecon
|
||||
#define dgesv_ dgesv
|
||||
#define dgetrs_ dgetrs
|
||||
#define dlange_ dlange
|
||||
#define dposv_ dposv
|
||||
#define dpotrs_ dpotrs
|
||||
#define dsyevx_ dsyevx
|
||||
#define dsygvx_ dsygvx
|
||||
#define dgeev_ dgeev
|
||||
#define addem2 ADDEM2
|
||||
#elif !defined(_AIX)
|
||||
#define addem2 addem2_
|
||||
#endif
|
||||
|
||||
#if defined(USE_CBLAS) || defined(USE_MKL)
|
||||
extern "C" {
|
||||
//! \brief Adds an element matrix \a eK into the system matrix \a sysK.
|
||||
//! \details This is a FORTRAN-77 subroutine in the SAM library.
|
||||
//! \sa SAM library documentation.
|
||||
void addem2_(const Real* eK, const Real* ttcc, const int* mpar,
|
||||
void addem2 (const Real* eK, const Real* ttcc, const int* mpar,
|
||||
const int* madof, const int* meqn, const int* mpmnpc,
|
||||
const int* mmnpc, const int* mpmceq, const int* mmceq,
|
||||
const int& iel, const int& nedof, const int& neq,
|
||||
const int& lpu, const int& nrhs, Real* sysK, Real* sysRHS,
|
||||
int* work, int& ierr);
|
||||
|
||||
//! \brief Estimates the reciprocal of the condition number of the matrix \b A.
|
||||
//! \details This is a FORTRAN-77 subroutine in the LAPack library.
|
||||
//! \sa LAPack library documentation.
|
||||
void dgecon_(const char* norm, const int* n, const Real* A, const int* lda,
|
||||
const Real* anorm, Real* rcond, Real* work, int* iwork, int* info);
|
||||
|
||||
//! \brief Solves the linear equation system \a A*x=b.
|
||||
//! \details This is a FORTRAN-77 subroutine in the LAPack library.
|
||||
//! \sa LAPack library documentation.
|
||||
void dgesv_(const int* n, const int* nrhs,
|
||||
Real* A, const int* lda, int* ipiv,
|
||||
Real* B, const int* ldb, int* info);
|
||||
|
||||
//! \brief Solves the equation system \a A*x=b when \a A is already factorized.
|
||||
//! \details This is a FORTRAN-77 subroutine in the LAPack library.
|
||||
//! \sa LAPack library documentation.
|
||||
void dgetrs_(const char* trans, const int* n, const int* nrhs,
|
||||
Real* A, const int* lda, int* ipiv,
|
||||
Real* B, const int* ldb, int* info);
|
||||
|
||||
//! \brief Returns a norm of the matrix \b A.
|
||||
//! \details This is a FORTRAN-77 subroutine in the LAPack library.
|
||||
//! \sa LAPack library documentation.
|
||||
double dlange_(const char* norm, const int* m, const int* n, const Real* A,
|
||||
const int* lda, Real* work);
|
||||
|
||||
//! \brief Solves the symmetric linear equation system \a A*x=b.
|
||||
//! \details This is a FORTRAN-77 subroutine in the LAPack library.
|
||||
//! \sa LAPack library documentation.
|
||||
void dposv_(const char* uplo, const int* n, const int* nrhs,
|
||||
Real* A, const int* lda, Real* B, const int* ldb, int* info);
|
||||
|
||||
//! \brief Solves the symmetric equation system \a A*x=b for prefactored \b A.
|
||||
//! \details This is a FORTRAN-77 subroutine in the LAPack library.
|
||||
//! \sa LAPack library documentation.
|
||||
void dpotrs_(const char* uplo, const int* n, const int* nrhs,
|
||||
Real* A, const int* lda, Real* B, const int* ldb, int* info);
|
||||
|
||||
//! \brief Solves the standard eigenproblem \a A*x=(lambda)*x.
|
||||
//! \details This is a FORTRAN-77 subroutine in the LAPack library.
|
||||
//! \sa LAPack library documentation.
|
||||
void dsyevx_(const char* jobz, const char* range, const char* uplo,
|
||||
const int* n, Real* a, const int* lda,
|
||||
const Real* vl, const Real* vu, const int* il, const int* iu,
|
||||
const Real* abstol, const int* m, Real* w, Real* z, const int* ldz,
|
||||
Real* work, const int* lwork, int* iwork, int* ifail, int* info);
|
||||
|
||||
//! \brief Solves the generalized eigenproblem \a A*x=(lambda)*B*x.
|
||||
//! \details This is a FORTRAN-77 subroutine in the LAPack library.
|
||||
//! \sa LAPack library documentation.
|
||||
void dsygvx_(const int* itype, const char* jobz, const char* range,
|
||||
const char* uplo, const int* n, Real* a, const int* lda,
|
||||
Real* b, const int* ldb, const Real* vl, const Real* vu,
|
||||
const int* il, const int* iu, const Real* abstol,
|
||||
const int* m, Real* w, Real* z, const int* ldz,
|
||||
Real* work, const int* lwork, int* iwork, int* ifail, int* info);
|
||||
|
||||
//! \brief Solves the non-symmetric eigenproblem \a A*x=(lambda)*x.
|
||||
//! \details This is a FORTRAN-77 subroutine in the LAPack library.
|
||||
//! \sa LAPack library documentation.
|
||||
void dgeev_(const char* jobvl, const char* jobvr,
|
||||
const int* n, Real* a, const int* lda,
|
||||
Real* wr, Real* wi, Real* vl, const int* ldvl,
|
||||
Real* vr, const int* ldvr,
|
||||
Real* work, const int* lwork, int* info);
|
||||
}
|
||||
#endif
|
||||
|
||||
|
||||
DenseMatrix::DenseMatrix (size_t m, size_t n, bool s) : myMat(m,n), ipiv(nullptr)
|
||||
DenseMatrix::DenseMatrix (size_t m, size_t n, bool s) : myMat(m,n)
|
||||
{
|
||||
ipiv = nullptr;
|
||||
symm = s && m == n;
|
||||
}
|
||||
|
||||
|
||||
DenseMatrix::DenseMatrix (const DenseMatrix& A) : ipiv(nullptr)
|
||||
DenseMatrix::DenseMatrix (const DenseMatrix& A)
|
||||
{
|
||||
myMat = A.myMat;
|
||||
ipiv = nullptr;
|
||||
symm = A.symm;
|
||||
if (A.ipiv)
|
||||
std::cerr <<"DenseMatrix constructor: Copying factored matrix"<< std::endl;
|
||||
}
|
||||
|
||||
|
||||
DenseMatrix::DenseMatrix (const RealArray& data, size_t nrows) : ipiv(nullptr)
|
||||
DenseMatrix::DenseMatrix (const RealArray& data, size_t nrows)
|
||||
{
|
||||
size_t ndata = data.size();
|
||||
if (nrows == 0) nrows = (size_t)sqrt((double)ndata);
|
||||
@@ -143,6 +62,7 @@ DenseMatrix::DenseMatrix (const RealArray& data, size_t nrows) : ipiv(nullptr)
|
||||
|
||||
myMat.resize(nrows,ncols);
|
||||
memcpy(myMat.ptr(),&data.front(),nrows*ncols*sizeof(Real));
|
||||
ipiv = nullptr;
|
||||
symm = false;
|
||||
}
|
||||
|
||||
@@ -281,7 +201,7 @@ bool DenseMatrix::assemble (const Matrix& eM, const SAM& sam, int e)
|
||||
#ifdef USE_F77SAM
|
||||
Real dummyRHS;
|
||||
int* work = new int[eM.rows()];
|
||||
addem2_(eM.ptr(), sam.ttcc, sam.mpar,
|
||||
addem2 (eM.ptr(), sam.ttcc, sam.mpar,
|
||||
sam.madof, sam.meqn, sam.mpmnpc, sam.mmnpc, sam.mpmceq, sam.mmceq,
|
||||
e, eM.rows(), sam.neq, 6, 0, myMat.ptr(), &dummyRHS, work, ierr);
|
||||
delete[] work;
|
||||
@@ -309,7 +229,7 @@ bool DenseMatrix::assemble (const Matrix& eM, const SAM& sam,
|
||||
int ierr = 0;
|
||||
#ifdef USE_F77SAM
|
||||
int* work = new int[eM.rows()];
|
||||
addem2_(eM.ptr(), sam.ttcc, sam.mpar,
|
||||
addem2 (eM.ptr(), sam.ttcc, sam.mpar,
|
||||
sam.madof, sam.meqn, sam.mpmnpc, sam.mmnpc, sam.mpmceq, sam.mmceq,
|
||||
e, eM.rows(), sam.neq, 6, 1, myMat.ptr(), Bptr->ptr(), work, ierr);
|
||||
delete[] work;
|
||||
@@ -471,45 +391,34 @@ bool DenseMatrix::solve (Real* B, size_t nrhs, Real* rcond)
|
||||
if (n > myMat.cols()) return false; // More equations than unknowns
|
||||
|
||||
const char* dsolv = symm ? "DGESV" : "DPOSV";
|
||||
#if USE_BLAS
|
||||
#ifdef HAS_BLAS
|
||||
int info = 0;
|
||||
int noncons = n;
|
||||
if (symm)
|
||||
{
|
||||
// Matrix is marked as symmetric - use Cholesky instead of full LU
|
||||
if (ipiv) {
|
||||
char uplo = 'U';
|
||||
int nr = nrhs;
|
||||
dpotrs_(&uplo,&noncons,&nr,myMat.ptr(),&noncons,B,&noncons,&info);
|
||||
} else
|
||||
if (ipiv)
|
||||
dpotrs ('U',n,nrhs,myMat.ptr(),n,B,n,info);
|
||||
else
|
||||
{
|
||||
ipiv = new int[1]; // dummy allocation to flag factorization reuse
|
||||
char uplo = 'U';
|
||||
int nr = nrhs;
|
||||
dposv_(&uplo,&noncons,&nr,myMat.ptr(),&noncons,B,&noncons,&info);
|
||||
dposv ('U',n,nrhs,myMat.ptr(),n,B,n,info);
|
||||
}
|
||||
}
|
||||
else if (ipiv) {
|
||||
char trans = 'N';
|
||||
int nr = nrhs;
|
||||
dgetrs_(&trans,&noncons,&nr,myMat.ptr(),&noncons,ipiv,B,&noncons,&info);
|
||||
} else
|
||||
else if (ipiv)
|
||||
dgetrs ('N',n,nrhs,myMat.ptr(),n,ipiv,B,n,info);
|
||||
else
|
||||
{
|
||||
Real anorm;
|
||||
if (rcond) { // Evaluate the 1-norm of the original LHS-matrix
|
||||
char norm = '1';
|
||||
anorm = dlange_(&norm,&noncons,&noncons,myMat.ptr(),&noncons,rcond);
|
||||
}
|
||||
if (rcond) // Evaluate the 1-norm of the original LHS-matrix
|
||||
anorm = dlange('1',n,n,myMat.ptr(),n,rcond);
|
||||
ipiv = new int[n];
|
||||
int nr = nrhs;
|
||||
dgesv_(&noncons,&nr,myMat.ptr(),&noncons,ipiv,B,&noncons,&info);
|
||||
dgesv (n,nrhs,myMat.ptr(),n,ipiv,B,n,info);
|
||||
if (rcond && info == 0)
|
||||
{
|
||||
// Estimate the condition number
|
||||
Real* work = new Real[4*n];
|
||||
int* iwork = new int[n];
|
||||
char norm='1';
|
||||
dgecon_(&norm,&noncons,myMat.ptr(),&noncons,&anorm,rcond,work,iwork,&info);
|
||||
dgecon ('1',n,myMat.ptr(),n,anorm,rcond,work,iwork,info);
|
||||
delete[] work;
|
||||
delete[] iwork;
|
||||
}
|
||||
@@ -535,20 +444,14 @@ bool DenseMatrix::solveEig (RealArray& val, Matrix& vec, int nv)
|
||||
if (n < 1 || nv < 1) return true; // No equations to solve
|
||||
if (n > myMat.cols()) return false;
|
||||
|
||||
#if USE_BLAS
|
||||
#ifdef HAS_BLAS
|
||||
std::cout <<" Solving dense eigenproblem using LAPACK::DSYEVX"<< std::endl;
|
||||
int m, info = 0;
|
||||
Real dummy = 0.0;
|
||||
Real abstol = 0.0;
|
||||
Real dummy = Real(0);
|
||||
Real abstol = Real(0);
|
||||
// Invoke with Lwork = -1 to estimate work space size
|
||||
char jobz='V';
|
||||
char range='I';
|
||||
char uplo='U';
|
||||
int il = 1;
|
||||
int lwork = -1;
|
||||
int noncons = n;
|
||||
dsyevx_(&jobz,&range,&uplo,&noncons,myMat.ptr(),&noncons,&dummy,&dummy,&il,&nv,
|
||||
&abstol,&m,&val.front(),vec.ptr(),&noncons,&dummy,&lwork,nullptr,nullptr,&info);
|
||||
dsyevx ('V','I','U',n,myMat.ptr(),n,dummy,dummy,1,nv,
|
||||
abstol,m,&val.front(),vec.ptr(),n,&dummy,-1,nullptr,nullptr,info);
|
||||
|
||||
if (info == 0)
|
||||
{
|
||||
@@ -559,12 +462,8 @@ bool DenseMatrix::solveEig (RealArray& val, Matrix& vec, int nv)
|
||||
val.resize(n);
|
||||
vec.resize(n,nv);
|
||||
// Solve the eigenproblem
|
||||
char jobz='V';
|
||||
char range='I';
|
||||
char uplo='U';
|
||||
int il = 1;
|
||||
dsyevx_(&jobz,&range,&uplo,&noncons,myMat.ptr(),&noncons,&dummy,&dummy,&il,&nv,
|
||||
&abstol,&m,&val.front(),vec.ptr(),&noncons,work,&Lwork,Iwork+n,Iwork,&info);
|
||||
dsyevx ('V','I','U',n,myMat.ptr(),n,dummy,dummy,1,nv,
|
||||
abstol,m,&val.front(),vec.ptr(),n,work,Lwork,Iwork+n,Iwork,info);
|
||||
delete[] work;
|
||||
delete[] Iwork;
|
||||
val.resize(nv);
|
||||
@@ -590,22 +489,15 @@ bool DenseMatrix::solveEig (DenseMatrix& B, RealArray& val, Matrix& vec, int nv,
|
||||
if (n < 1 || nv < 1) return true; // No equations to solve
|
||||
if (n > myMat.cols()) return false;
|
||||
|
||||
#if USE_BLAS
|
||||
#ifdef HAS_BLAS
|
||||
std::cout <<" Solving dense eigenproblem using LAPACK::DSYGVX"<< std::endl;
|
||||
int m, info = 0;
|
||||
Real dummy = 0.0;
|
||||
Real abstol = 0.0;
|
||||
int itype = 1;
|
||||
char jobz = 'V';
|
||||
char range = 'I';
|
||||
char uplo = 'U';
|
||||
int il = 1;
|
||||
int lwork = -1;
|
||||
int noncons = n;
|
||||
Real dummy = Real(0);
|
||||
Real abstol = Real(0);
|
||||
// Invoke with Lwork = -1 to estimate work space size
|
||||
dsygvx_(&itype,&jobz,&range,&uplo,&noncons,myMat.ptr(),&noncons,B.myMat.ptr(),&noncons,
|
||||
&dummy,&dummy,&il,&nv,&abstol,&m,&val.front(),vec.ptr(),&noncons,
|
||||
&dummy,&lwork,nullptr,nullptr,&info);
|
||||
dsygvx (1,'V','I','U',n,myMat.ptr(),n,B.myMat.ptr(),n,
|
||||
dummy,dummy,1,nv,abstol,m,&val.front(),vec.ptr(),n,
|
||||
&dummy,-1,nullptr,nullptr,info);
|
||||
|
||||
if (info == 0)
|
||||
{
|
||||
@@ -616,9 +508,9 @@ bool DenseMatrix::solveEig (DenseMatrix& B, RealArray& val, Matrix& vec, int nv,
|
||||
val.resize(n);
|
||||
vec.resize(n,nv);
|
||||
// Solve the eigenproblem
|
||||
dsygvx_(&itype,&jobz,&range,&uplo,&noncons,myMat.ptr(),&noncons,B.myMat.ptr(),&noncons,
|
||||
&dummy,&dummy,&il,&nv,&abstol,&m,&val.front(),vec.ptr(),&noncons,
|
||||
work,&Lwork,Iwork+n,Iwork,&info);
|
||||
dsygvx (1,'V','I','U',n,myMat.ptr(),n,B.myMat.ptr(),n,
|
||||
dummy,dummy,1,nv,abstol,m,&val.front(),vec.ptr(),n,
|
||||
work,Lwork,Iwork+n,Iwork,info);
|
||||
delete[] work;
|
||||
delete[] Iwork;
|
||||
val.resize(nv);
|
||||
@@ -642,19 +534,15 @@ bool DenseMatrix::solveEig (DenseMatrix& B, RealArray& val, Matrix& vec, int nv,
|
||||
|
||||
bool DenseMatrix::solveEigNon (RealArray& r_val, RealArray& c_val)
|
||||
{
|
||||
const int n = myMat.rows();
|
||||
const size_t n = myMat.rows();
|
||||
if (n < 1) return true; // No equations to solve
|
||||
|
||||
#if USE_BLAS
|
||||
#ifdef HAS_BLAS
|
||||
int info = 0;
|
||||
Real dummy = 0.0;
|
||||
char jobvl = 'N';
|
||||
int ldvl = 1;
|
||||
int lwork = -1;
|
||||
int noncons = n;
|
||||
Real dummy = Real(0);
|
||||
// Invoke with Lwork = -1 to estimate work space size
|
||||
dgeev_(&jobvl,&jobvl,&noncons,myMat.ptr(),&noncons,&r_val.front(),&c_val.front(),
|
||||
&dummy,&ldvl,&dummy,&ldvl,&dummy,&lwork,&info);
|
||||
dgeev ('N','N',n,myMat.ptr(),n,&r_val.front(),&c_val.front(),
|
||||
&dummy,1,&dummy,1,&dummy,-1,info);
|
||||
|
||||
if (info == 0)
|
||||
{
|
||||
@@ -664,8 +552,8 @@ bool DenseMatrix::solveEigNon (RealArray& r_val, RealArray& c_val)
|
||||
r_val.resize(n);
|
||||
c_val.resize(n);
|
||||
// Solve the eigenproblem
|
||||
dgeev_(&jobvl,&jobvl,&noncons,myMat.ptr(),&noncons,&r_val.front(),&c_val.front(),
|
||||
&dummy,&ldvl,&dummy,&ldvl,work,&Lwork,&info);
|
||||
dgeev ('N','N',n,myMat.ptr(),n,&r_val.front(),&c_val.front(),
|
||||
&dummy,1,&dummy,1,work,Lwork,info);
|
||||
delete[] work;
|
||||
if (info == 0) return true;
|
||||
}
|
||||
@@ -676,14 +564,16 @@ bool DenseMatrix::solveEigNon (RealArray& r_val, RealArray& c_val)
|
||||
return false;
|
||||
}
|
||||
|
||||
DenseMatrix operator*(Real alpha, const DenseMatrix& A)
|
||||
|
||||
DenseMatrix operator* (Real alpha, const DenseMatrix& A)
|
||||
{
|
||||
DenseMatrix B(A);
|
||||
B.getMat() *= alpha;
|
||||
return B;
|
||||
}
|
||||
|
||||
DenseMatrix operator*(const DenseMatrix& A, Real alpha)
|
||||
|
||||
DenseMatrix operator* (const DenseMatrix& A, Real alpha)
|
||||
{
|
||||
DenseMatrix B(A);
|
||||
B.getMat() *= alpha;
|
||||
|
||||
239
src/LinAlg/LAPack.h
Normal file
239
src/LinAlg/LAPack.h
Normal file
@@ -0,0 +1,239 @@
|
||||
// $Id$
|
||||
//==============================================================================
|
||||
//!
|
||||
//! \file LAPack.h
|
||||
//!
|
||||
//! \date Jan 10 2016
|
||||
//!
|
||||
//! \author Knut Morten Okstad / SINTEF
|
||||
//!
|
||||
//! \brief LAPack support for various platforms.
|
||||
//!
|
||||
//==============================================================================
|
||||
|
||||
#ifndef UTL_LAPACK_H
|
||||
#define UTL_LAPACK_H
|
||||
|
||||
#include "BLAS.h"
|
||||
|
||||
#if HAS_BLAS == 3
|
||||
#define Subroutine static inline void
|
||||
|
||||
//! \brief Estimates the reciprocal of the condition number of the matrix \b A.
|
||||
//! \details This is a FORTRAN-77 subroutine in the LAPack library.
|
||||
//! \sa LAPack library documentation.
|
||||
Subroutine dgecon (char norm, int n, const Real* A, int lda,
|
||||
Real anorm, Real* rcond, Real* work, int* iwork, int& info)
|
||||
{ dgecon_(&norm,&n,const_cast<Real*>(A),&lda,&anorm,rcond,work,iwork,&info); }
|
||||
|
||||
//! \brief Solves the linear equation system \a A*x=b.
|
||||
//! \details This is a FORTRAN-77 subroutine in the LAPack library.
|
||||
//! \sa LAPack library documentation.
|
||||
Subroutine dgesv (int n, int nrhs,
|
||||
Real* A, int lda, int* ipiv,
|
||||
Real* B, int ldb, int& info)
|
||||
{ dgesv_(&n,&nrhs,A,&lda,ipiv,B,&ldb,&info); }
|
||||
|
||||
//! \brief Computes an LU factorization of a general M-by-N matrix A.
|
||||
//! \details This is a FORTRAN-77 subroutine in the LAPack library.
|
||||
//! \sa LAPack library documentation.
|
||||
Subroutine dgetrf (int m, int n, Real* A, int lda,
|
||||
int* ipiv, int& info)
|
||||
{ dgetrf_(&m,&n,A,&lda,ipiv,&info); }
|
||||
|
||||
//! \brief Computes the inverse of a matrix based on its LU factorization.
|
||||
//! \details This is a FORTRAN-77 subroutine in the LAPack library.
|
||||
//! \sa LAPack library documentation.
|
||||
Subroutine dgetri (int n, Real* A, int lda,
|
||||
int* ipiv, Real* work, int lwork, int& info)
|
||||
{ dgetri_(&n,A,&lda,ipiv,work,&lwork,&info); }
|
||||
|
||||
//! \brief Solves the equation system \a A*x=b when \a A is already factorized.
|
||||
//! \details This is a FORTRAN-77 subroutine in the LAPack library.
|
||||
//! \sa LAPack library documentation.
|
||||
Subroutine dgetrs (char trans, int n, int nrhs,
|
||||
Real* A, int lda, int* ipiv,
|
||||
Real* B, int ldb, int& info)
|
||||
{ dgetrs_(&trans,&n,&nrhs,A,&lda,ipiv,B,&ldb,&info); }
|
||||
|
||||
//! \brief Returns a norm of the matrix \b A.
|
||||
//! \details This is a FORTRAN-77 function in the LAPack library.
|
||||
//! \sa LAPack library documentation.
|
||||
static inline double dlange (char norm, int m, int n, const Real* A,
|
||||
int lda, Real* work)
|
||||
{ return dlange_(&norm,&m,&n,const_cast<Real*>(A),&lda,work); }
|
||||
|
||||
//! \brief Solves the symmetric linear equation system \a A*x=b.
|
||||
//! \details This is a FORTRAN-77 subroutine in the LAPack library.
|
||||
//! \sa LAPack library documentation.
|
||||
Subroutine dposv (char uplo, int n, int nrhs,
|
||||
Real* A, int lda, Real* B, int ldb, int& info)
|
||||
{ dposv_(&uplo,&n,&nrhs,A,&lda,B,&ldb,&info); }
|
||||
|
||||
//! \brief Solves the symmetric equation system \a A*x=b for prefactored \b A.
|
||||
//! \details This is a FORTRAN-77 subroutine in the LAPack library.
|
||||
//! \sa LAPack library documentation.
|
||||
Subroutine dpotrs (char uplo, int n, int nrhs,
|
||||
Real* A, int lda, Real* B, int ldb, int& info)
|
||||
{ dpotrs_(&uplo,&n,&nrhs,A,&lda,B,&ldb,&info); }
|
||||
|
||||
//! \brief Solves the standard eigenproblem \a A*x=(lambda)*x.
|
||||
//! \details This is a FORTRAN-77 subroutine in the LAPack library.
|
||||
//! \sa LAPack library documentation.
|
||||
Subroutine dsyev (char jobz, char uplo,
|
||||
int n, double* A, int lda, double* w,
|
||||
double* work, int lwork, int& info)
|
||||
{ dsyev_(&jobz,&uplo,&n,A,&lda,w,work,&lwork,&info); }
|
||||
|
||||
//! \brief Solves the standard eigenproblem \a A*x=(lambda)*x.
|
||||
//! \details This is a FORTRAN-77 subroutine in the LAPack library.
|
||||
//! \sa LAPack library documentation.
|
||||
Subroutine dsyevx (char jobz, char range, char uplo,
|
||||
int n, Real* A, int lda,
|
||||
Real vl, Real vu, int il, int iu,
|
||||
Real abstol, int& m, Real* w, Real* z, int ldz,
|
||||
Real* work, int lwork, int* iwork, int* ifail, int& info)
|
||||
{ dsyevx_(&jobz,&range,&uplo,&n,A,&lda,&vl,&vu,&il,&iu,&abstol,
|
||||
&m,w,z,&ldz,work,&lwork,iwork,ifail,&info); }
|
||||
|
||||
//! \brief Solves the generalized eigenproblem \a A*x=(lambda)*B*x.
|
||||
//! \details This is a FORTRAN-77 subroutine in the LAPack library.
|
||||
//! \sa LAPack library documentation.
|
||||
Subroutine dsygvx (int itype, char jobz, char range,
|
||||
char uplo, int n, Real* A, int lda,
|
||||
Real* B, int ldb, Real vl, Real vu,
|
||||
int il, int iu, Real abstol,
|
||||
int& m, Real* w, Real* z, int ldz,
|
||||
Real* work, int lwork, int* iwork, int* ifail, int& info)
|
||||
{ dsygvx_(&itype,&jobz,&range,&uplo,&n,A,&lda,B,&ldb,&vl,&vu,&il,&iu,&abstol,
|
||||
&m,w,z,&ldz,work,&lwork,iwork,ifail,&info); }
|
||||
|
||||
//! \brief Solves the non-symmetric eigenproblem \a A*x=(lambda)*x.
|
||||
//! \details This is a FORTRAN-77 subroutine in the LAPack library.
|
||||
//! \sa LAPack library documentation.
|
||||
Subroutine dgeev (char jobvl, char jobvr,
|
||||
int n, Real* A, int lda,
|
||||
Real* wr, Real* wi, Real* vl, int ldvl,
|
||||
Real* vr, int ldvr,
|
||||
Real* work, int lwork, int& info)
|
||||
{ dgeev_(&jobvl,&jobvr,&n,A,&lda,wr,wi,vl,&ldvl,vr,&ldvr,work,&lwork,&info); }
|
||||
|
||||
#else
|
||||
|
||||
#if defined(_WIN32) && !defined(__MINGW32__) && !defined(__MINGW64__)
|
||||
#define dgecon DGECON
|
||||
#define dgesv DGESV
|
||||
#define dgetrf DGETRF
|
||||
#define dgetri DGETRI
|
||||
#define dgetrs DGETRS
|
||||
#define dlange DLANGE
|
||||
#define dposv DPOSV
|
||||
#define dpotrs DPOTRS
|
||||
#define dsyev DSYEV
|
||||
#define dsyevx DSYEVX
|
||||
#define dsygvx DSYGVX
|
||||
#define dgeev DGEEV
|
||||
#elif !defined(_AIX)
|
||||
#define dgecon dgecon_
|
||||
#define dgesv dgesv_
|
||||
#define dgetrf dgetrf_
|
||||
#define dgetri dgetri_
|
||||
#define dgetrs dgetrs_
|
||||
#define dlange dlange_
|
||||
#define dposv dposv_
|
||||
#define dpotrs dpotrs_
|
||||
#define dsyev dsyev_
|
||||
#define dsyevx dsyevx_
|
||||
#define dsygvx dsygvx_
|
||||
#define dgeev dgeev_
|
||||
#endif
|
||||
|
||||
extern "C" {
|
||||
//! \brief Estimates the reciprocal of the condition number of the matrix \b A.
|
||||
//! \details This is a FORTRAN-77 subroutine in the LAPack library.
|
||||
//! \sa LAPack library documentation.
|
||||
void dgecon (const char& norm, const int& n, const Real* A, const int& lda,
|
||||
const Real& anorm, Real* rcond, Real* work, int* iwork, int& info);
|
||||
|
||||
//! \brief Solves the linear equation system \a A*x=b.
|
||||
//! \details This is a FORTRAN-77 subroutine in the LAPack library.
|
||||
//! \sa LAPack library documentation.
|
||||
void dgesv (const int& n, const int& nrhs,
|
||||
Real* A, const int& lda, int* ipiv,
|
||||
Real* B, const int& ldb, int& info);
|
||||
|
||||
//! \brief Computes an LU factorization of a general M-by-N matrix A.
|
||||
//! \details This is a FORTRAN-77 subroutine in the LAPack library.
|
||||
//! \sa LAPack library documentation.
|
||||
void dgetrf (const int& m, const int& n, Real* A, const int& lda,
|
||||
int* ipiv, int& info);
|
||||
|
||||
//! \brief Computes the inverse of a matrix based on its LU factorization.
|
||||
//! \details This is a FORTRAN-77 subroutine in the LAPack library.
|
||||
//! \sa LAPack library documentation.
|
||||
void dgetri (const int& n, Real* A, const int& lda,
|
||||
int* ipiv, Real* work, const int& lwork, int& info);
|
||||
|
||||
//! \brief Solves the equation system \a A*x=b when \a A is already factorized.
|
||||
//! \details This is a FORTRAN-77 subroutine in the LAPack library.
|
||||
//! \sa LAPack library documentation.
|
||||
void dgetrs (const char& trans, const int& n, const int& nrhs,
|
||||
Real* A, const int& lda, int* ipiv,
|
||||
Real* B, const int& ldb, int& info);
|
||||
|
||||
//! \brief Returns a norm of the matrix \b A.
|
||||
//! \details This is a FORTRAN-77 function in the LAPack library.
|
||||
//! \sa LAPack library documentation.
|
||||
double dlange (const char& norm, const int& m, const int& n, const Real* A,
|
||||
const int& lda, Real* work);
|
||||
|
||||
//! \brief Solves the symmetric linear equation system \a A*x=b.
|
||||
//! \details This is a FORTRAN-77 subroutine in the LAPack library.
|
||||
//! \sa LAPack library documentation.
|
||||
void dposv (const char& uplo, const int& n, const int& nrhs,
|
||||
Real* A, const int& lda, Real* B, const int& ldb, int& info);
|
||||
|
||||
//! \brief Solves the symmetric equation system \a A*x=b for prefactored \b A.
|
||||
//! \details This is a FORTRAN-77 subroutine in the LAPack library.
|
||||
//! \sa LAPack library documentation.
|
||||
void dpotrs (const char& uplo, const int& n, const int& nrhs,
|
||||
Real* A, const int& lda, Real* B, const int& ldb, int& info);
|
||||
|
||||
//! \brief Solves the standard eigenproblem \a A*x=(lambda)*x.
|
||||
//! \details This is a FORTRAN-77 subroutine in the LAPack library.
|
||||
//! \sa LAPack library documentation.
|
||||
void dsyev (const char& jobz, const char& uplo,
|
||||
const int& n, double* A, const int& lda, double* w,
|
||||
double* work, const int& lwork, int& info);
|
||||
|
||||
//! \brief Solves the standard eigenproblem \a A*x=(lambda)*x.
|
||||
//! \details This is a FORTRAN-77 subroutine in the LAPack library.
|
||||
//! \sa LAPack library documentation.
|
||||
void dsyevx (const char& jobz, const char& range, const char& uplo,
|
||||
const int& n, Real* a, const int& lda,
|
||||
const Real& vl, const Real& vu, const int& il, const int& iu,
|
||||
const Real& abstol, int& m, Real* w, Real* z, const int& ldz,
|
||||
Real* work, const int& lwork, int* iwork, int* ifail, int& info);
|
||||
|
||||
//! \brief Solves the generalized eigenproblem \a A*x=(lambda)*B*x.
|
||||
//! \details This is a FORTRAN-77 subroutine in the LAPack library.
|
||||
//! \sa LAPack library documentation.
|
||||
void dsygvx (const int& itype, const char& jobz, const char& range,
|
||||
const char& uplo, const int& n, Real* a, const int& lda,
|
||||
Real* b, const int& ldb, const Real& vl, const Real& vu,
|
||||
const int& il, const int& iu, const Real& abstol,
|
||||
int& m, Real* w, Real* z, const int& ldz,
|
||||
Real* work, const int& lwork, int* iwork, int* ifail, int& info);
|
||||
|
||||
//! \brief Solves the non-symmetric eigenproblem \a A*x=(lambda)*x.
|
||||
//! \details This is a FORTRAN-77 subroutine in the LAPack library.
|
||||
//! \sa LAPack library documentation.
|
||||
void dgeev (const char& jobvl, const char& jobvr,
|
||||
const int& n, Real* a, const int& lda,
|
||||
Real* wr, Real* wi, Real* vl, const int& ldvl,
|
||||
Real* vr, const int& ldvr,
|
||||
Real* work, const int& lwork, int& info);
|
||||
}
|
||||
#endif
|
||||
|
||||
#endif
|
||||
@@ -15,29 +15,7 @@
|
||||
//==============================================================================
|
||||
|
||||
#include "MatVec.h"
|
||||
|
||||
#if defined(_WIN32) && !defined(__MINGW32__) && !defined(__MINGW64__)
|
||||
#define dgetrf_ DGETRF
|
||||
#define dgetri_ DGITRI
|
||||
#elif defined(_AIX)
|
||||
#define dgetrf_ dgetrf
|
||||
#define dgetri_ dgetri
|
||||
#endif
|
||||
|
||||
#if defined(USE_CBLAS) || defined(USE_MKL)
|
||||
extern "C" {
|
||||
//! \brief Computes an LU factorization of a general M-by-N matrix A.
|
||||
//! \details This is a FORTRAN-77 subroutine in the LAPack library.
|
||||
//! \sa LAPack library documentation.
|
||||
void dgetrf_(const int* M, const int* N, double* A, const int* LDA,
|
||||
int* IPIV, int* INFO);
|
||||
//! \brief Computes the inverse of a matrix based on its LU factorization.
|
||||
//! \details This is a FORTRAN-77 subroutine in the LAPack library.
|
||||
//! \sa LAPack library documentation.
|
||||
void dgetri_(const int* N, double* A, const int* LDA,
|
||||
int* IPIV, double* WORK, const int* LWORK, int* INFO);
|
||||
}
|
||||
#endif
|
||||
#include "LAPack.h"
|
||||
|
||||
|
||||
// Set default values on some global variables.
|
||||
@@ -133,12 +111,11 @@ bool utl::invert (Matrix& A)
|
||||
if (A.rows() <= 3)
|
||||
return A.inverse() > 0.0;
|
||||
|
||||
#if USE_BLAS
|
||||
#ifdef HAS_BLAS
|
||||
// Use LAPack/BLAS for larger matrices
|
||||
int INFO, N = A.rows() > A.cols() ? A.rows() : A.cols();
|
||||
int* IPIV = new int[N];
|
||||
int r = A.rows();
|
||||
dgetrf_(&N,&N,A.ptr(),&r,IPIV,&INFO);
|
||||
dgetrf (N,N,A.ptr(),A.rows(),IPIV,INFO);
|
||||
if (INFO != 0)
|
||||
{
|
||||
delete[] IPIV;
|
||||
@@ -146,11 +123,9 @@ bool utl::invert (Matrix& A)
|
||||
return false;
|
||||
}
|
||||
double NWORK;
|
||||
int m1 = -1;
|
||||
dgetri_(&N,A.ptr(),&r,IPIV,&NWORK,&m1,&INFO);
|
||||
dgetri (N,A.ptr(),A.rows(),IPIV,&NWORK,-1,INFO);
|
||||
double* WORK = new double[int(NWORK)];
|
||||
int nwork = NWORK;
|
||||
dgetri_(&N,A.ptr(),&r,IPIV,WORK,&nwork,&INFO);
|
||||
dgetri (N,A.ptr(),A.rows(),IPIV,WORK,int(NWORK),INFO);
|
||||
delete[] IPIV;
|
||||
delete[] WORK;
|
||||
if (INFO == 0) return true;
|
||||
|
||||
@@ -10,8 +10,8 @@
|
||||
//! \brief Simple template classes for dense rectangular matrices and vectors.
|
||||
//! \details The classes have some algebraic operators defined, such that the
|
||||
//! class type \a T has to be of a numerical type, i.e., \a float or \a double.
|
||||
//! The multiplication methods are implemented based on the CBLAS library if
|
||||
//! either of the macro symbols USE_CBLAS or USE_MKL are defined.
|
||||
//! The multiplication methods are implemented based on the BLAS library if
|
||||
//! an implementation is provided (HAS_BLAS is defined, see BLAS.h).
|
||||
//! Otherwise, inlined methods are used.
|
||||
//!
|
||||
//==============================================================================
|
||||
@@ -24,18 +24,7 @@
|
||||
#include <algorithm>
|
||||
#include <cstring>
|
||||
#include <cmath>
|
||||
#ifdef USE_MKL
|
||||
#include <mkl_cblas.h>
|
||||
#elif defined(USE_ACCELERATE)
|
||||
#include <Accelerate/Accelerate.h>
|
||||
#elif defined(USE_CBLAS)
|
||||
extern "C"
|
||||
{
|
||||
#include <cblas.h>
|
||||
}
|
||||
#endif
|
||||
|
||||
#define USE_BLAS defined(USE_CBLAS) || defined(USE_MKL) || defined(USE_ACCELERATE)
|
||||
#include "BLAS.h"
|
||||
|
||||
#ifdef INDEX_CHECK
|
||||
#if INDEX_CHECK > 1
|
||||
@@ -157,6 +146,17 @@ namespace utl //! General utility classes and functions.
|
||||
return this->operator=(Y).relax(alfa,X);
|
||||
}
|
||||
|
||||
//! \brief Dot product between \a *this and another vector.
|
||||
//! \param[in] v The vector to dot this vector with
|
||||
//! \param[in] nv Length of the vector \a v
|
||||
//! \param[in] off1 Offset for this vector
|
||||
//! \param[in] inc1 Increment for this vector
|
||||
//! \param[in] off2 Offset for vector \b v
|
||||
//! \param[in] inc2 Increment for vector \b v
|
||||
T dot(const T* v, size_t nv,
|
||||
size_t off1 = 0, int inc1 = 1,
|
||||
size_t off2 = 0, int inc2 = 1) const;
|
||||
|
||||
//! \brief Dot product between \a *this and another vector.
|
||||
//! \param[in] v The vector to dot this vector with
|
||||
//! \param[in] off1 Offset for this vector
|
||||
@@ -165,7 +165,10 @@ namespace utl //! General utility classes and functions.
|
||||
//! \param[in] inc2 Increment for vector \b v
|
||||
T dot(const std::vector<T>& v,
|
||||
size_t off1 = 0, int inc1 = 1,
|
||||
size_t off2 = 0, int inc2 = 1) const;
|
||||
size_t off2 = 0, int inc2 = 1) const
|
||||
{
|
||||
return this->dot(&v.front(),v.size(),off1,inc1,off2,inc2);
|
||||
}
|
||||
|
||||
//! \brief Return the Euclidean norm of the vector.
|
||||
//! \param[in] off Index offset relative to the first vector component
|
||||
@@ -382,7 +385,7 @@ namespace utl //! General utility classes and functions.
|
||||
size_t nc = transposed ? block.rows() : block.cols();
|
||||
for (size_t i = 1; i <= nr && i+row-1 <= nrow; i++)
|
||||
for (size_t j = 1; j <= nc && j+col-1 <= ncol; j++)
|
||||
(*this)(i+row-1,j+col-1) = transposed ? block(j,i) : block(i,j);
|
||||
elem[i+row-2+nrow*(j+col-2)] = transposed ? block(j,i) : block(i,j);
|
||||
}
|
||||
|
||||
//! \brief Create a diagonal matrix.
|
||||
@@ -815,7 +818,7 @@ namespace utl //! General utility classes and functions.
|
||||
};
|
||||
|
||||
|
||||
#if USE_BLAS
|
||||
#ifdef HAS_BLAS
|
||||
//============================================================================
|
||||
//=== BLAS-implementation of the matrix/vector multiplication methods ====
|
||||
//============================================================================
|
||||
@@ -835,23 +838,23 @@ namespace utl //! General utility classes and functions.
|
||||
}
|
||||
|
||||
template<> inline
|
||||
float vector<float>::dot(const std::vector<float>& v,
|
||||
float vector<float>::dot(const float* v, size_t nv,
|
||||
size_t o1, int i1, size_t o2, int i2) const
|
||||
{
|
||||
int n1 = i1 > 1 || i1 < -1 ? this->size()/abs(i1) : this->size()-o1;
|
||||
int n2 = i2 > 1 || i2 < -1 ? v.size()/abs(i2) : v.size()-o2;
|
||||
int n2 = i2 > 1 || i2 < -1 ? nv/abs(i2) : nv-o2;
|
||||
int n = n1 < n2 ? n1 : n2;
|
||||
return cblas_sdot(n,this->ptr()+o1,i1,&v.front()+o2,i2);
|
||||
return cblas_sdot(n,this->ptr()+o1,i1,v+o2,i2);
|
||||
}
|
||||
|
||||
template<> inline
|
||||
double vector<double>::dot(const std::vector<double>& v,
|
||||
double vector<double>::dot(const double* v, size_t nv,
|
||||
size_t o1, int i1, size_t o2, int i2) const
|
||||
{
|
||||
int n1 = i1 > 1 || i1 < -1 ? this->size()/abs(i1) : this->size()-o1;
|
||||
int n2 = i2 > 1 || i2 < -1 ? v.size()/abs(i2) : v.size()-o2;
|
||||
int n2 = i2 > 1 || i2 < -1 ? nv/abs(i2) : nv-o2;
|
||||
int n = n1 < n2 ? n1 : n2;
|
||||
return cblas_ddot(n,this->ptr()+o1,i1,&v.front()+o2,i2);
|
||||
return cblas_ddot(n,this->ptr()+o1,i1,v+o2,i2);
|
||||
}
|
||||
|
||||
template<> inline
|
||||
@@ -1176,12 +1179,12 @@ namespace utl //! General utility classes and functions.
|
||||
}
|
||||
|
||||
template<class T> inline
|
||||
T vector<T>::dot(const std::vector<T>& v,
|
||||
T vector<T>::dot(const T* v, size_t nv,
|
||||
size_t o1, int i1, size_t o2, int i2) const
|
||||
{
|
||||
size_t i, j;
|
||||
T dotprod = T(0);
|
||||
for (i = o1, j = o2; i < this->size() && j < v.size(); i += i1, j += i2)
|
||||
for (i = o1, j = o2; i < this->size() && j < nv; i += i1, j += i2)
|
||||
dotprod += this->operator[](i) * v[j];
|
||||
return dotprod;
|
||||
}
|
||||
|
||||
@@ -13,28 +13,11 @@
|
||||
|
||||
#include "Tensor.h"
|
||||
#include "Vec3.h"
|
||||
#include "matrix.h"
|
||||
#include "LAPack.h"
|
||||
#include <array>
|
||||
#include <algorithm>
|
||||
#include <cstring>
|
||||
|
||||
#if defined(_WIN32) && !defined(__MINGW32__) && !defined(__MINGW64__)
|
||||
#define dsyev_ DSYEV
|
||||
#elif defined(_AIX)
|
||||
#define dsyev_ dsyev
|
||||
#endif
|
||||
|
||||
#if defined(USE_CBLAS) || defined(USE_MKL)
|
||||
extern "C" {
|
||||
//! \brief Solves the standard eigenproblem \a A*x=(lambda)*x.
|
||||
//! \details This is a FORTRAN-77 subroutine in the LAPack library.
|
||||
//! \sa LAPack library documentation.
|
||||
void dsyev_(const char* jobz, const char* uplo,
|
||||
const int* n, double* a, const int* lda, double* w,
|
||||
double* work, const int* lwork, int* info);
|
||||
}
|
||||
#endif
|
||||
|
||||
#ifndef epsZ
|
||||
//! \brief Zero tolerance for the incremental rotations.
|
||||
#define epsZ 1.0e-16
|
||||
@@ -1132,7 +1115,7 @@ bool SymmTensor::principal (Vec3& p) const
|
||||
bool SymmTensor::principal (Vec3& p, Vec3* pdir, int ndir) const
|
||||
{
|
||||
p = 0.0;
|
||||
#ifdef USE_BLAS
|
||||
#ifdef HAS_BLAS
|
||||
// Set up the upper triangle of the symmetric tensor
|
||||
double A[9], W[3];
|
||||
if (n == 3)
|
||||
@@ -1157,11 +1140,7 @@ bool SymmTensor::principal (Vec3& p, Vec3* pdir, int ndir) const
|
||||
int info = 0;
|
||||
const int Lwork = 12;
|
||||
double Work[Lwork];
|
||||
char jobz='V';
|
||||
char uplo='U';
|
||||
int noncons = n;
|
||||
int Noncons = Lwork;
|
||||
dsyev_(&jobz,&uplo,&noncons,A,&noncons,W,Work,&Noncons,&info);
|
||||
dsyev ('V','U',n,A,n,W,Work,Lwork,info);
|
||||
if (info)
|
||||
{
|
||||
std::cerr <<" *** LAPack::dsyev: info ="<< info;
|
||||
|
||||
Reference in New Issue
Block a user