Page Menu
Home
HEPForge
Search
Configure Global Search
Log In
Files
F19244355
lapack_interface.icc
No One
Temporary
Actions
View File
Edit File
Delete File
View Transforms
Subscribe
Award Token
Flag For Later
Size
38 KB
Referenced Files
None
Subscribers
None
lapack_interface.icc
View Options
#include <vector>
#include <cassert>
#include <stdexcept>
#include <sstream>
#include "npstat/nm/lapack.h"
namespace npstat {
namespace Private {
inline void check_lapack_status(const int INFO, const char* fcn1,
const char* fcnlap)
{
if (INFO)
{
std::ostringstream os;
os << "In npstat::" << fcn1 << ": LAPACK routine "
<< fcnlap << " failed with status " << INFO;
throw std::runtime_error(os.str());
}
}
}
template<typename Numeric>
inline void invert_posdef_sym_matrix(const Numeric* /* in */,
unsigned /* dim */,
Numeric* /* out */)
{
throw std::invalid_argument(
"In npstat::invert_posdef_sym_matrix: "
"matrix inversion is not supported for this template parameter");
}
template<typename Numeric>
inline void invert_sym_matrix(const Numeric* /* in */,
unsigned /* dim */,
Numeric* /* out */)
{
throw std::invalid_argument(
"In npstat::invert_sym_matrix: "
"matrix inversion is not supported for this template parameter");
}
template<typename Numeric>
inline void invert_general_matrix(const Numeric* /* in */,
unsigned /* dim */,
Numeric* /* out */)
{
throw std::invalid_argument(
"In npstat::invert_general_matrix: "
"matrix inversion is not supported for this template parameter");
}
template<typename Numeric>
inline void sym_matrix_eigenvalues(const Numeric* /* in */,
unsigned /* dim */,
Numeric* /* eigenvalues */)
{
throw std::invalid_argument(
"In npstat::sym_matrix_eigenvalues: calculation "
"of eigenvalues is not supported for this template parameter");
}
template<typename Numeric>
inline void sym_matrix_eigenvalues_dc(const Numeric* /* in */,
unsigned /* dim */,
Numeric* /* eigenvalues */)
{
throw std::invalid_argument(
"In npstat::sym_matrix_eigenvalues_dc: calculation "
"of eigenvalues is not supported for this template parameter");
}
template<typename Numeric>
inline void sym_matrix_eigenvalues_rrr(const Numeric* /* in */,
unsigned /* dim */,
Numeric* /* eigenvalues */)
{
throw std::invalid_argument(
"In npstat::sym_matrix_eigenvalues_rrr: calculation "
"of eigenvalues is not supported for this template parameter");
}
template<typename Numeric>
inline void sym_matrix_eigensystem(const Numeric* /* in */,
unsigned /* dim */,
Numeric* /* eigenvalues */,
Numeric* /* eigenvectors */)
{
throw std::invalid_argument(
"In npstat::sym_matrix_eigensystem: calculation "
"of eigenvectors is not supported for this template parameter");
}
template<typename Numeric>
inline void sym_matrix_eigensystem_dc(const Numeric* /* in */,
unsigned /* dim */,
Numeric* /* eigenvalues */,
Numeric* /* eigenvectors */)
{
throw std::invalid_argument(
"In npstat::sym_matrix_eigensystem_dc: calculation "
"of eigenvectors is not supported for this template parameter");
}
template<typename Numeric>
inline void sym_matrix_eigensystem_rrr(const Numeric* /* in */,
unsigned /* dim */,
Numeric* /* eigenvalues */,
Numeric* /* eigenvectors */)
{
throw std::invalid_argument(
"In npstat::sym_matrix_eigensystem_rrr: calculation "
"of eigenvectors is not supported for this template parameter");
}
template<typename Numeric>
inline Numeric lu_decomposition_matrix_det(const Numeric* /* in */,
unsigned /* dim */)
{
throw std::invalid_argument(
"In npstat::lu_decomposition_matrix_det: calculation "
"of determinant is not supported for this template parameter");
}
template<typename Numeric>
inline bool solve_linear_system(const Numeric* /* in */,
unsigned /* dim */,
const Numeric* /* rhs */,
Numeric* /* solution */)
{
throw std::invalid_argument(
"In npstat::solve_linear_system: solution of "
"linear systems is not supported for this template parameter");
}
template<>
inline void sym_matrix_eigenvalues<double>(const double* in,
const unsigned dim,
double* eigenvalues)
{
// The following code is not thread safe!
static std::vector<double> bufv;
assert(in);
assert(dim);
assert(eigenvalues);
char UPLO[] = "U", JOBZ[] = "N";
int N = dim, LDA = dim, INFO = 0, LWORK = 34*dim;
if (bufv.size() < dim*dim + LWORK)
bufv.resize(dim*dim + LWORK);
double* buf = &bufv[0];
double* work = buf + dim*dim;
for (unsigned i=0; i<dim; ++i)
for (unsigned j=0; j<=i; ++j)
{
const unsigned idx = i*dim + j;
buf[idx] = (in[idx] + in[j*dim + i])/2.0;
}
dsyev_(JOBZ, UPLO, &N, buf, &LDA, eigenvalues,
work, &LWORK, &INFO, 1, 1);
Private::check_lapack_status(INFO, "sym_matrix_eigenvalues", "DSYEV");
}
template<>
inline void sym_matrix_eigensystem<double>(const double* in,
const unsigned dim,
double* eigenvalues,
double* eigenvectors)
{
// The following code is not thread safe!
static std::vector<double> bufv;
assert(in);
assert(dim);
assert(eigenvalues);
assert(eigenvectors);
char UPLO[] = "U", JOBZ[] = "V";
int N = dim, LDA = dim, INFO = 0, LWORK = 34*dim;
for (unsigned i=0; i<dim; ++i)
for (unsigned j=0; j<=i; ++j)
{
const unsigned idx = i*dim + j;
eigenvectors[idx] = (in[idx] + in[j*dim + i])/2.0;
}
if (bufv.size() < static_cast<unsigned>(LWORK))
bufv.resize(LWORK);
dsyev_(JOBZ, UPLO, &N, eigenvectors, &LDA, eigenvalues,
&bufv[0], &LWORK, &INFO, 1, 1);
Private::check_lapack_status(INFO, "sym_matrix_eigensystem", "DSYEV");
}
template<>
inline void sym_matrix_eigenvalues<float>(const float* in,
const unsigned dim,
float* eigenvalues)
{
// The following code is not thread safe!
static std::vector<float> bufv;
assert(in);
assert(dim);
assert(eigenvalues);
char UPLO[] = "U", JOBZ[] = "N";
int N = dim, LDA = dim, INFO = 0, LWORK = 34*dim;
if (bufv.size() < dim*dim + LWORK)
bufv.resize(dim*dim + LWORK);
float* buf = &bufv[0];
float* work = buf + dim*dim;
for (unsigned i=0; i<dim; ++i)
for (unsigned j=0; j<=i; ++j)
{
const unsigned idx = i*dim + j;
buf[idx] = (in[idx] + in[j*dim + i])/2.0;
}
ssyev_(JOBZ, UPLO, &N, buf, &LDA, eigenvalues,
work, &LWORK, &INFO, 1, 1);
Private::check_lapack_status(INFO, "sym_matrix_eigenvalues", "SSYEV");
}
template<>
inline void sym_matrix_eigensystem<float>(const float* in,
const unsigned dim,
float* eigenvalues,
float* eigenvectors)
{
// The following code is not thread safe!
static std::vector<float> bufv;
assert(in);
assert(dim);
assert(eigenvalues);
assert(eigenvectors);
char UPLO[] = "U", JOBZ[] = "V";
int N = dim, LDA = dim, INFO = 0, LWORK = 34*dim;
for (unsigned i=0; i<dim; ++i)
for (unsigned j=0; j<=i; ++j)
{
const unsigned idx = i*dim + j;
eigenvectors[idx] = (in[idx] + in[j*dim + i])/2.0;
}
if (bufv.size() < static_cast<unsigned>(LWORK))
bufv.resize(LWORK);
ssyev_(JOBZ, UPLO, &N, eigenvectors, &LDA, eigenvalues,
&bufv[0], &LWORK, &INFO, 1, 1);
Private::check_lapack_status(INFO, "sym_matrix_eigensystem", "SSYEV");
}
template<>
inline void sym_matrix_eigenvalues_dc<double>(const double* in,
const unsigned dim,
double* eigenvalues)
{
// The following code is not thread safe!
static std::vector<double> bufv;
static std::vector<int> ibufv;
static unsigned maxDim = 0;
assert(in);
assert(dim);
assert(eigenvalues);
if (dim > maxDim)
{
// Calculate optimal buffer sizes and allocate the memory
char UPLO[] = "U", JOBZ[] = "N";
int N = dim, LDA = dim, INFO = 0, LWORK = -1, LIWORK = -1;
double work = 0.0, buf = 0.0;
int IWORK = 0;
dsyevd_(JOBZ, UPLO, &N, &buf, &LDA, eigenvalues,
&work, &LWORK, &IWORK, &LIWORK, &INFO, 1, 1);
assert(INFO == 0);
assert(IWORK > 0);
assert(work > 0.0);
bufv.resize(dim*dim + static_cast<unsigned>(work));
ibufv.resize(IWORK);
maxDim = dim;
}
char UPLO[] = "U", JOBZ[] = "N";
int N = dim, LDA = dim, INFO = 0, LWORK = bufv.size() - dim*dim;
int LIWORK = ibufv.size();
double* buf = &bufv[0];
double* work = buf + dim*dim;
int* IWORK = &ibufv[0];
for (unsigned i=0; i<dim; ++i)
for (unsigned j=0; j<=i; ++j)
{
const unsigned idx = i*dim + j;
buf[idx] = (in[idx] + in[j*dim + i])/2.0;
}
dsyevd_(JOBZ, UPLO, &N, buf, &LDA, eigenvalues,
work, &LWORK, IWORK, &LIWORK, &INFO, 1, 1);
Private::check_lapack_status(INFO, "sym_matrix_eigenvalues_dc",
"DSYEVD");
}
template<>
inline void sym_matrix_eigenvalues_dc<float>(const float* in,
const unsigned dim,
float* eigenvalues)
{
// The following code is not thread safe!
static std::vector<float> bufv;
static std::vector<int> ibufv;
static unsigned maxDim = 0;
assert(in);
assert(dim);
assert(eigenvalues);
if (dim > maxDim)
{
// Calculate optimal buffer sizes and allocate the memory
char UPLO[] = "U", JOBZ[] = "N";
int N = dim, LDA = dim, INFO = 0, LWORK = -1, LIWORK = -1;
float work = 0.f, buf = 0.f;
int IWORK = 0;
ssyevd_(JOBZ, UPLO, &N, &buf, &LDA, eigenvalues,
&work, &LWORK, &IWORK, &LIWORK, &INFO, 1, 1);
assert(INFO == 0);
assert(IWORK > 0);
assert(work > 0.f);
bufv.resize(dim*dim + static_cast<unsigned>(work));
ibufv.resize(IWORK);
maxDim = dim;
}
char UPLO[] = "U", JOBZ[] = "N";
int N = dim, LDA = dim, INFO = 0, LWORK = bufv.size() - dim*dim;
int LIWORK = ibufv.size();
float* buf = &bufv[0];
float* work = buf + dim*dim;
int* IWORK = &ibufv[0];
for (unsigned i=0; i<dim; ++i)
for (unsigned j=0; j<=i; ++j)
{
const unsigned idx = i*dim + j;
buf[idx] = (in[idx] + in[j*dim + i])/2.0;
}
ssyevd_(JOBZ, UPLO, &N, buf, &LDA, eigenvalues,
work, &LWORK, IWORK, &LIWORK, &INFO, 1, 1);
Private::check_lapack_status(INFO, "sym_matrix_eigenvalues_dc",
"SSYEVD");
}
template<>
inline void sym_matrix_eigenvalues_rrr<double>(const double* in,
const unsigned dim,
double* eigenvalues)
{
// The following code is not thread safe!
static std::vector<double> bufv;
static std::vector<int> ibufv;
static unsigned maxDim = 0;
assert(in);
assert(dim);
assert(eigenvalues);
if (dim > maxDim)
{
// Calculate optimal buffer sizes and allocate the memory
char UPLO[] = "U", RANGE[] = "A", JOBZ[] = "N";
int N = dim, LDA = dim, INFO = 0, LWORK = -1, LIWORK = -1;
double work = 0.0, buf = 0.0, VL = 0.0, VU = 0.0, Z = 0.0;
int IWORK = 0, IL = 0, IU = 0, M = 0, LDZ = dim, ISUPPZ = 0;
double ABSTOL = 0.0;
dsyevr_(JOBZ, RANGE, UPLO, &N, &buf, &LDA,
&VL, &VU, &IL, &IU, &ABSTOL, &M, eigenvalues,
&Z, &LDZ, &ISUPPZ,
&work, &LWORK, &IWORK, &LIWORK, &INFO, 1, 1, 1);
assert(INFO == 0);
assert(IWORK > 0);
assert(work > 0.0);
bufv.resize(dim*dim + static_cast<unsigned>(work));
ibufv.resize(IWORK + 2*dim);
maxDim = dim;
}
char UPLO[] = "U", RANGE[] = "A", JOBZ[] = "N";
int N = dim, LDA = dim, INFO = 0, LWORK = bufv.size() - dim*dim;
int LIWORK = ibufv.size() - 2*dim;
double VL = 0.0, VU = 0.0, Z = 0.0;
int IL = 0, IU = 0, M = 0, LDZ = dim;
double ABSTOL = 0.0;
double* buf = &bufv[0];
double* work = buf + dim*dim;
int* ISUPPZ = &ibufv[0];
int* IWORK = ISUPPZ + 2*dim;
for (unsigned i=0; i<dim; ++i)
for (unsigned j=0; j<=i; ++j)
{
const unsigned idx = i*dim + j;
buf[idx] = (in[idx] + in[j*dim + i])/2.0;
}
dsyevr_(JOBZ, RANGE, UPLO, &N, buf, &LDA,
&VL, &VU, &IL, &IU, &ABSTOL, &M, eigenvalues,
&Z, &LDZ, ISUPPZ,
work, &LWORK, IWORK, &LIWORK, &INFO, 1, 1, 1);
Private::check_lapack_status(INFO, "sym_matrix_eigenvalues_rrr",
"DSYEVR");
assert(M == N);
}
template<>
inline void sym_matrix_eigenvalues_rrr<float>(const float* in,
const unsigned dim,
float* eigenvalues)
{
// The following code is not thread safe!
static std::vector<float> bufv;
static std::vector<int> ibufv;
static unsigned maxDim = 0;
assert(in);
assert(dim);
assert(eigenvalues);
if (dim > maxDim)
{
// Calculate optimal buffer sizes and allocate the memory
char UPLO[] = "U", RANGE[] = "A", JOBZ[] = "N";
int N = dim, LDA = dim, INFO = 0, LWORK = -1, LIWORK = -1;
float work = 0.f, buf = 0.f, VL = 0.f, VU = 0.f, Z = 0.f;
int IWORK = 0, IL = 0, IU = 0, M = 0, LDZ = dim, ISUPPZ = 0;
float ABSTOL = 0.f;
ssyevr_(JOBZ, RANGE, UPLO, &N, &buf, &LDA,
&VL, &VU, &IL, &IU, &ABSTOL, &M, eigenvalues,
&Z, &LDZ, &ISUPPZ,
&work, &LWORK, &IWORK, &LIWORK, &INFO, 1, 1, 1);
assert(INFO == 0);
assert(IWORK > 0);
assert(work > 0.f);
bufv.resize(dim*dim + static_cast<unsigned>(work));
ibufv.resize(IWORK + 2*dim);
maxDim = dim;
}
char UPLO[] = "U", RANGE[] = "A", JOBZ[] = "N";
int N = dim, LDA = dim, INFO = 0, LWORK = bufv.size() - dim*dim;
int LIWORK = ibufv.size() - 2*dim;
float VL = 0.f, VU = 0.f, Z = 0.f;
int IL = 0, IU = 0, M = 0, LDZ = dim;
float ABSTOL = 0.f;
float* buf = &bufv[0];
float* work = buf + dim*dim;
int* ISUPPZ = &ibufv[0];
int* IWORK = ISUPPZ + 2*dim;
for (unsigned i=0; i<dim; ++i)
for (unsigned j=0; j<=i; ++j)
{
const unsigned idx = i*dim + j;
buf[idx] = (in[idx] + in[j*dim + i])/2.0;
}
ssyevr_(JOBZ, RANGE, UPLO, &N, buf, &LDA,
&VL, &VU, &IL, &IU, &ABSTOL, &M, eigenvalues,
&Z, &LDZ, ISUPPZ,
work, &LWORK, IWORK, &LIWORK, &INFO, 1, 1, 1);
Private::check_lapack_status(INFO, "sym_matrix_eigenvalues_rrr",
"SSYEVR");
assert(M == N);
}
template<>
inline void sym_matrix_eigensystem_dc<double>(const double* in,
const unsigned dim,
double* eigenvalues,
double* eigenvectors)
{
// The following code is not thread safe!
static std::vector<double> bufv;
static std::vector<int> ibufv;
static unsigned maxDim = 0;
assert(in);
assert(dim);
assert(eigenvalues);
assert(eigenvectors);
if (dim > maxDim)
{
// Calculate optimal buffer sizes and allocate the memory
char UPLO[] = "U", JOBZ[] = "V";
int N = dim, LDA = dim, INFO = 0, LWORK = -1, LIWORK = -1;
double work = 0.0, buf = 0.0;
int IWORK = 0;
dsyevd_(JOBZ, UPLO, &N, &buf, &LDA, eigenvalues,
&work, &LWORK, &IWORK, &LIWORK, &INFO, 1, 1);
assert(INFO == 0);
assert(IWORK > 0);
assert(work > 0.0);
bufv.resize(static_cast<unsigned>(work));
ibufv.resize(IWORK);
maxDim = dim;
}
char UPLO[] = "U", JOBZ[] = "V";
int N = dim, LDA = dim, INFO = 0, LWORK = bufv.size();
int LIWORK = ibufv.size();
double* work = &bufv[0];
int* IWORK = &ibufv[0];
for (unsigned i=0; i<dim; ++i)
for (unsigned j=0; j<=i; ++j)
{
const unsigned idx = i*dim + j;
eigenvectors[idx] = (in[idx] + in[j*dim + i])/2.0;
}
dsyevd_(JOBZ, UPLO, &N, eigenvectors, &LDA, eigenvalues,
work, &LWORK, IWORK, &LIWORK, &INFO, 1, 1);
Private::check_lapack_status(INFO, "sym_matrix_eigensystem_dc",
"DSYEVD");
}
template<>
inline void sym_matrix_eigensystem_dc<float>(const float* in,
const unsigned dim,
float* eigenvalues,
float* eigenvectors)
{
// The following code is not thread safe!
static std::vector<float> bufv;
static std::vector<int> ibufv;
static unsigned maxDim = 0;
assert(in);
assert(dim);
assert(eigenvalues);
assert(eigenvectors);
if (dim > maxDim)
{
// Calculate optimal buffer sizes and allocate the memory
char UPLO[] = "U", JOBZ[] = "V";
int N = dim, LDA = dim, INFO = 0, LWORK = -1, LIWORK = -1;
float work = 0.f, buf = 0.f;
int IWORK = 0;
ssyevd_(JOBZ, UPLO, &N, &buf, &LDA, eigenvalues,
&work, &LWORK, &IWORK, &LIWORK, &INFO, 1, 1);
assert(INFO == 0);
assert(IWORK > 0);
assert(work > 0.f);
bufv.resize(static_cast<unsigned>(work));
ibufv.resize(IWORK);
maxDim = dim;
}
char UPLO[] = "U", JOBZ[] = "V";
int N = dim, LDA = dim, INFO = 0, LWORK = bufv.size();
int LIWORK = ibufv.size();
float* work = &bufv[0];
int* IWORK = &ibufv[0];
for (unsigned i=0; i<dim; ++i)
for (unsigned j=0; j<=i; ++j)
{
const unsigned idx = i*dim + j;
eigenvectors[idx] = (in[idx] + in[j*dim + i])/2.0;
}
ssyevd_(JOBZ, UPLO, &N, eigenvectors, &LDA, eigenvalues,
work, &LWORK, IWORK, &LIWORK, &INFO, 1, 1);
Private::check_lapack_status(INFO, "sym_matrix_eigensystem_dc",
"SSYEVD");
}
template<>
inline void sym_matrix_eigensystem_rrr<double>(const double* in,
const unsigned dim,
double* eigenvalues,
double* eigenvectors)
{
// The following code is not thread safe!
static std::vector<double> bufv;
static std::vector<int> ibufv;
static unsigned maxDim = 0;
assert(in);
assert(dim);
assert(eigenvalues);
assert(eigenvectors);
if (dim > maxDim)
{
// Calculate optimal buffer sizes and allocate the memory
char UPLO[] = "U", RANGE[] = "A", JOBZ[] = "V";
int N = dim, LDA = dim, INFO = 0, LWORK = -1, LIWORK = -1;
double work = 0.0, buf = 0.0, VL = 0.0, VU = 0.0, Z = 0.0;
int IWORK = 0, IL = 0, IU = 0, M = 0, LDZ = dim, ISUPPZ = 0;
double ABSTOL = 0.0;
dsyevr_(JOBZ, RANGE, UPLO, &N, &buf, &LDA,
&VL, &VU, &IL, &IU, &ABSTOL, &M, eigenvalues,
&Z, &LDZ, &ISUPPZ,
&work, &LWORK, &IWORK, &LIWORK, &INFO, 1, 1, 1);
assert(INFO == 0);
assert(IWORK > 0);
assert(work > 0.0);
bufv.resize(dim*dim + static_cast<unsigned>(work));
ibufv.resize(IWORK + 2*dim);
maxDim = dim;
}
char UPLO[] = "U", RANGE[] = "A", JOBZ[] = "V";
int N = dim, LDA = dim, INFO = 0, LWORK = bufv.size() - dim*dim;
int LIWORK = ibufv.size() - 2*dim;
double VL = 0.0, VU = 0.0;
int IL = 0, IU = 0, M = 0, LDZ = dim;
double ABSTOL = 0.0;
double* buf = &bufv[0];
double* work = buf + dim*dim;
int* ISUPPZ = &ibufv[0];
int* IWORK = ISUPPZ + 2*dim;
for (unsigned i=0; i<dim; ++i)
for (unsigned j=0; j<=i; ++j)
{
const unsigned idx = i*dim + j;
buf[idx] = (in[idx] + in[j*dim + i])/2.0;
}
dsyevr_(JOBZ, RANGE, UPLO, &N, buf, &LDA,
&VL, &VU, &IL, &IU, &ABSTOL, &M, eigenvalues,
eigenvectors, &LDZ, ISUPPZ,
work, &LWORK, IWORK, &LIWORK, &INFO, 1, 1, 1);
Private::check_lapack_status(INFO, "sym_matrix_eigensystem_rrr",
"DSYEVR");
assert(M == N);
}
template<>
inline void sym_matrix_eigensystem_rrr<float>(const float* in,
const unsigned dim,
float* eigenvalues,
float* eigenvectors)
{
// The following code is not thread safe!
static std::vector<float> bufv;
static std::vector<int> ibufv;
static unsigned maxDim = 0;
assert(in);
assert(dim);
assert(eigenvalues);
assert(eigenvectors);
if (dim > maxDim)
{
// Calculate optimal buffer sizes and allocate the memory
char UPLO[] = "U", RANGE[] = "A", JOBZ[] = "V";
int N = dim, LDA = dim, INFO = 0, LWORK = -1, LIWORK = -1;
float work = 0.f, buf = 0.f, VL = 0.f, VU = 0.f, Z = 0.f;
int IWORK = 0, IL = 0, IU = 0, M = 0, LDZ = dim, ISUPPZ = 0;
float ABSTOL = 0.f;
ssyevr_(JOBZ, RANGE, UPLO, &N, &buf, &LDA,
&VL, &VU, &IL, &IU, &ABSTOL, &M, eigenvalues,
&Z, &LDZ, &ISUPPZ,
&work, &LWORK, &IWORK, &LIWORK, &INFO, 1, 1, 1);
assert(INFO == 0);
assert(IWORK > 0);
assert(work > 0.f);
bufv.resize(dim*dim + static_cast<unsigned>(work));
ibufv.resize(IWORK + 2*dim);
maxDim = dim;
}
char UPLO[] = "U", RANGE[] = "A", JOBZ[] = "V";
int N = dim, LDA = dim, INFO = 0, LWORK = bufv.size() - dim*dim;
int LIWORK = ibufv.size() - 2*dim;
float VL = 0.f, VU = 0.f;
int IL = 0, IU = 0, M = 0, LDZ = dim;
float ABSTOL = 0.f;
float* buf = &bufv[0];
float* work = buf + dim*dim;
int* ISUPPZ = &ibufv[0];
int* IWORK = ISUPPZ + 2*dim;
for (unsigned i=0; i<dim; ++i)
for (unsigned j=0; j<=i; ++j)
{
const unsigned idx = i*dim + j;
buf[idx] = (in[idx] + in[j*dim + i])/2.0;
}
ssyevr_(JOBZ, RANGE, UPLO, &N, buf, &LDA,
&VL, &VU, &IL, &IU, &ABSTOL, &M, eigenvalues,
eigenvectors, &LDZ, ISUPPZ,
work, &LWORK, IWORK, &LIWORK, &INFO, 1, 1, 1);
Private::check_lapack_status(INFO, "sym_matrix_eigensystem_rrr",
"SSYEVR");
assert(M == N);
}
template<>
inline void invert_posdef_sym_matrix<double>(const double* in,
const unsigned dim,
double* out)
{
// The following code is not thread safe!
static std::vector<double> bufv;
assert(in);
assert(out);
if (bufv.size() < dim*dim)
bufv.resize(dim*dim);
double* buf = &bufv[0];
for (unsigned i=0; i<dim; ++i)
for (unsigned j=0; j<=i; ++j)
{
const unsigned idx = i*dim + j;
buf[idx] = (in[idx] + in[j*dim + i])/2.0;
}
int N = dim, LDA = dim, INFO = 0;
char U[] = "U";
dpotrf_(U, &N, buf, &LDA, &INFO, 1);
Private::check_lapack_status(INFO,"invert_posdef_sym_matrix","DPOTRF");
dpotri_(U, &N, buf, &LDA, &INFO, 1);
Private::check_lapack_status(INFO,"invert_posdef_sym_matrix","DPOTRI");
for (unsigned i=0; i<dim; ++i)
for (unsigned j=0; j<=i; ++j)
{
const unsigned idx = i*dim + j;
out[idx] = buf[idx];
out[j*dim + i] = buf[idx];
}
}
template<>
inline void invert_posdef_sym_matrix<float>(const float* in,
const unsigned dim,
float* out)
{
// The following code is not thread safe!
static std::vector<float> bufv;
assert(in);
assert(out);
if (bufv.size() < dim*dim)
bufv.resize(dim*dim);
float* buf = &bufv[0];
for (unsigned i=0; i<dim; ++i)
for (unsigned j=0; j<=i; ++j)
{
const unsigned idx = i*dim + j;
buf[idx] = (in[idx] + in[j*dim + i])/2.0f;
}
int N = dim, LDA = dim, INFO = 0;
char U[] = "U";
spotrf_(U, &N, buf, &LDA, &INFO, 1);
Private::check_lapack_status(INFO,"invert_posdef_sym_matrix","SPOTRF");
spotri_(U, &N, buf, &LDA, &INFO, 1);
Private::check_lapack_status(INFO,"invert_posdef_sym_matrix","SPOTRI");
for (unsigned i=0; i<dim; ++i)
for (unsigned j=0; j<=i; ++j)
{
const unsigned idx = i*dim + j;
out[idx] = buf[idx];
out[j*dim + i] = buf[idx];
}
}
template<>
inline void invert_sym_matrix<double>(const double* in,
const unsigned dim,
double* out)
{
// The following code is not thread safe!
static std::vector<double> bufv;
static std::vector<int> bufPiv;
assert(in);
assert(out);
if (bufv.size() < dim*(dim + 32U))
bufv.resize(dim*(dim + 32U));
double* buf = &bufv[0];
double* WORK = buf + dim*dim;
for (unsigned i=0; i<dim; ++i)
for (unsigned j=0; j<=i; ++j)
{
const unsigned idx = i*dim + j;
buf[idx] = (in[idx] + in[j*dim + i])/2.0;
}
if (bufPiv.size() < dim)
bufPiv.resize(dim);
int* IPIV = &bufPiv[0];
int N = dim, LDA = dim, INFO = 0, LWORK = dim*32;
char U[] = "U";
dsytrf_(U, &N, buf, &LDA, IPIV, WORK, &LWORK, &INFO, 1);
Private::check_lapack_status(INFO, "invert_sym_matrix", "DSYTRF");
dsytri_(U, &N, buf, &LDA, IPIV, WORK, &INFO, 1);
Private::check_lapack_status(INFO, "invert_sym_matrix", "DSYTRI");
for (unsigned i=0; i<dim; ++i)
for (unsigned j=0; j<=i; ++j)
{
const unsigned idx = i*dim + j;
out[idx] = buf[idx];
out[j*dim + i] = buf[idx];
}
}
template<>
inline void invert_sym_matrix<float>(const float* in,
const unsigned dim,
float* out)
{
// The following code is not thread safe!
static std::vector<float> bufv;
static std::vector<int> bufPiv;
assert(in);
assert(out);
if (bufv.size() < dim*(dim + 32U))
bufv.resize(dim*(dim + 32U));
float* buf = &bufv[0];
float* WORK = buf + dim*dim;
for (unsigned i=0; i<dim; ++i)
for (unsigned j=0; j<=i; ++j)
{
const unsigned idx = i*dim + j;
buf[idx] = (in[idx] + in[j*dim + i])/2.0;
}
if (bufPiv.size() < dim)
bufPiv.resize(dim);
int* IPIV = &bufPiv[0];
int N = dim, LDA = dim, INFO = 0, LWORK = dim*32;
char U[] = "U";
ssytrf_(U, &N, buf, &LDA, IPIV, WORK, &LWORK, &INFO, 1);
Private::check_lapack_status(INFO, "invert_sym_matrix", "SSYTRF");
ssytri_(U, &N, buf, &LDA, IPIV, WORK, &INFO, 1);
Private::check_lapack_status(INFO, "invert_sym_matrix", "SSYTRI");
for (unsigned i=0; i<dim; ++i)
for (unsigned j=0; j<=i; ++j)
{
const unsigned idx = i*dim + j;
out[idx] = buf[idx];
out[j*dim + i] = buf[idx];
}
}
template<>
inline double lu_decomposition_matrix_det<double>(const double* in,
const unsigned dim)
{
// The following code is not thread safe!
static std::vector<double> bufv;
static std::vector<int> bufPiv;
assert(in);
const unsigned len = dim*dim;
if (bufv.size() < len)
bufv.resize(len);
double* buf = &bufv[0];
for (unsigned i=0; i<len; ++i)
buf[i] = in[i];
if (bufPiv.size() < dim)
bufPiv.resize(dim);
int* IPIV = &bufPiv[0];
int M = dim, N = dim, LDA = dim, INFO = 0;
dgetrf_(&M, &N, buf, &LDA, IPIV, &INFO);
if (INFO > 0)
return 0.0;
else
{
Private::check_lapack_status(INFO, "lu_decomposition_matrix_det",
"DGETRF");
long double prod = 1.0L;
for (unsigned i=0; i<dim; ++i)
prod *= buf[i*dim + i];
for (int i=0; i<static_cast<int>(dim); ++i)
if (IPIV[i] != i + 1)
prod *= -1.0L;
return prod;
}
}
template<>
inline float lu_decomposition_matrix_det<float>(const float* in,
const unsigned dim)
{
// The following code is not thread safe!
static std::vector<float> bufv;
static std::vector<int> bufPiv;
assert(in);
const unsigned len = dim*dim;
if (bufv.size() < len)
bufv.resize(len);
float* buf = &bufv[0];
for (unsigned i=0; i<len; ++i)
buf[i] = in[i];
if (bufPiv.size() < dim)
bufPiv.resize(dim);
int* IPIV = &bufPiv[0];
int M = dim, N = dim, LDA = dim, INFO = 0;
sgetrf_(&M, &N, buf, &LDA, IPIV, &INFO);
if (INFO > 0)
return 0.0;
else
{
Private::check_lapack_status(INFO, "lu_decomposition_matrix_det",
"SGETRF");
long double prod = 1.0L;
for (unsigned i=0; i<dim; ++i)
prod *= buf[i*dim + i];
for (int i=0; i<static_cast<int>(dim); ++i)
if (IPIV[i] != i + 1)
prod *= -1.0L;
return prod;
}
}
template<>
inline void invert_general_matrix<double>(const double* in,
const unsigned dim,
double* out)
{
// The following code is not thread safe!
static std::vector<double> bufv;
static std::vector<int> bufPiv;
assert(in);
assert(out);
const unsigned len = dim*dim;
if (bufv.size() < dim*(dim + 32U))
bufv.resize(dim*(dim + 32U));
double* buf = &bufv[0];
double* WORK = buf + dim*dim;
for (unsigned i=0; i<len; ++i)
buf[i] = in[i];
if (bufPiv.size() < dim)
bufPiv.resize(dim);
int* IPIV = &bufPiv[0];
int M = dim, N = dim, LDA = dim, INFO = 0, LWORK = dim*32;
dgetrf_(&M, &N, buf, &LDA, IPIV, &INFO);
Private::check_lapack_status(INFO, "invert_general_matrix", "DGETRF");
dgetri_(&N, buf, &LDA, IPIV, WORK, &LWORK, &INFO);
Private::check_lapack_status(INFO, "invert_general_matrix", "DGETRI");
for (unsigned i=0; i<len; ++i)
out[i] = buf[i];
}
template<>
inline void invert_general_matrix<float>(const float* in,
const unsigned dim,
float* out)
{
// The following code is not thread safe!
static std::vector<float> bufv;
static std::vector<int> bufPiv;
assert(in);
assert(out);
const unsigned len = dim*dim;
if (bufv.size() < dim*(dim + 32U))
bufv.resize(dim*(dim + 32U));
float* buf = &bufv[0];
float* WORK = buf + dim*dim;
for (unsigned i=0; i<len; ++i)
buf[i] = in[i];
if (bufPiv.size() < dim)
bufPiv.resize(dim);
int* IPIV = &bufPiv[0];
int M = dim, N = dim, LDA = dim, INFO = 0, LWORK = dim*32;
sgetrf_(&M, &N, buf, &LDA, IPIV, &INFO);
Private::check_lapack_status(INFO, "invert_general_matrix", "SGETRF");
sgetri_(&N, buf, &LDA, IPIV, WORK, &LWORK, &INFO);
Private::check_lapack_status(INFO, "invert_general_matrix", "SGETRI");
for (unsigned i=0; i<len; ++i)
out[i] = buf[i];
}
template<>
inline bool solve_linear_system(const double* in, const unsigned dim,
const double* rhs, double* solution)
{
// The following code is not thread safe!
static std::vector<double> bufv;
static std::vector<int> bufPiv;
assert(in);
assert(rhs);
assert(solution);
const unsigned len = dim*dim;
if (bufv.size() < len)
bufv.resize(len);
double* buf = &bufv[0];
for (unsigned i=0; i<len; ++i)
buf[i] = in[i];
if (bufPiv.size() < dim)
bufPiv.resize(dim);
int* IPIV = &bufPiv[0];
int M = dim, N = dim, LDA = dim, INFO = 0;
dgetrf_(&M, &N, buf, &LDA, IPIV, &INFO);
if (INFO > 0)
return false;
else
{
Private::check_lapack_status(INFO, "solve_linear_system", "DGETRF");
for (unsigned i=0; i<len; ++i)
solution[i] = rhs[i];
int NRHS = 1, LDB = dim;
char T[] = "T";
dgetrs_(T, &N, &NRHS, buf, &LDA, IPIV,
solution, &LDB, &INFO, 1);
Private::check_lapack_status(INFO, "solve_linear_system", "DGETRS");
return true;
}
}
template<>
inline bool solve_linear_system(const float* in, const unsigned dim,
const float* rhs, float* solution)
{
// The following code is not thread safe!
static std::vector<float> bufv;
static std::vector<int> bufPiv;
assert(in);
assert(rhs);
assert(solution);
const unsigned len = dim*dim;
if (bufv.size() < len)
bufv.resize(len);
float* buf = &bufv[0];
for (unsigned i=0; i<len; ++i)
buf[i] = in[i];
if (bufPiv.size() < dim)
bufPiv.resize(dim);
int* IPIV = &bufPiv[0];
int M = dim, N = dim, LDA = dim, INFO = 0;
sgetrf_(&M, &N, buf, &LDA, IPIV, &INFO);
if (INFO > 0)
return false;
else
{
Private::check_lapack_status(INFO, "solve_linear_system", "SGETRF");
for (unsigned i=0; i<len; ++i)
solution[i] = rhs[i];
int NRHS = 1, LDB = dim;
char T[] = "T";
sgetrs_(T, &N, &NRHS, buf, &LDA, IPIV,
solution, &LDB, &INFO, 1);
Private::check_lapack_status(INFO, "solve_linear_system", "SGETRS");
return true;
}
}
}
File Metadata
Details
Attached
Mime Type
text/x-c
Expires
Tue, Sep 30, 4:41 AM (1 d, 3 h)
Storage Engine
blob
Storage Format
Raw Data
Storage Handle
6561604
Default Alt Text
lapack_interface.icc (38 KB)
Attached To
Mode
rNPSTATSVN npstatsvn
Attached
Detach File
Event Timeline
Log In to Comment