OilfieldToolsNet/opm-core
4
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

Honour exact interfaces of LAPACK routines.

Browse Source 
The size type, MAT_SIZE_T, is not necessarily equivalent to 'int'.
This commit is contained in:
Bård Skaflestad 2012-08-24 16:49:50 +02:00
parent 8a3593c372
commit 203728109f
1 changed files with 28 additions and 26 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

54
tests/test_lapack.cpp
View File

@ -24,34 +24,36 @@
namespace {
struct BandMatrixCoeff
{
BandMatrixCoeff(int N, int ku, int kl) : ku_(ku), kl_(kl), nrow_(2*kl + ku + 1), N_(N) {
BandMatrixCoeff(MAT_SIZE_T N, MAT_SIZE_T ku, MAT_SIZE_T kl)
: ku_(ku), kl_(kl), nrow_(2 * kl + ku + 1), N_(N)
{
}
// compute the position where to store the coefficient of a matrix A_{i,j} (i,j=0,...,N-1)
// in a array which is sent to the band matrix solver of LAPACK.
int operator ()(int i, int j) {
int operator ()(MAT_SIZE_T i, MAT_SIZE_T j) {
return kl_ + ku_ + i - j + j*nrow_;
}
const int ku_;
const int kl_;
const int nrow_;
const int N_;
const MAT_SIZE_T ku_;
const MAT_SIZE_T kl_;
const MAT_SIZE_T nrow_;
const MAT_SIZE_T N_;
};
} //end anonymous namespace
int main()
{
const int N = 5;
const int nrhs = 1;
const MAT_SIZE_T N = 5;
const MAT_SIZE_T nrhs = 1;
double DU[N-1] = { 2.1, -1.0, 1.9, 8.0 };
double D[N] = { 3.0, 2.3, -5.0, -0.9, 7.1 };
double DL[N-1] = { 3.4, 3.6, 7.0, -6.0 };
double B[N] = { 2.7, -0.5, 2.6, 0.6, 2.7 };
// double B[N] = { 2.7, -0.5, 2.6, 0.6, 2.7 };
int info = 0;
MAT_SIZE_T info = 0;
dgtsv_(&N, &nrhs, DL, D, DU, B, &N, &info);
if (info == 0) {
for (int i = 0; i < N; ++i) {
@ -64,12 +66,12 @@ int main()
std::cout << std::endl;
//test of dgbsv_
int ldb = N;
int lda = N;
std::vector<int> ipiv(N, 0);
MAT_SIZE_T ldb = N;
MAT_SIZE_T lda = N;
std::vector<MAT_SIZE_T> ipiv(N, 0);
std::vector<double> AA(N*N, 0.);
std::vector<double> BB(N, 0.);
for (int i = 0; i < N; ++i) {
for (MAT_SIZE_T i = 0; i < N; ++i) {
AA[i + i*N] = 10.;
if (i > 0) {
AA[i + (i - 1)*N] = i;
@ -78,31 +80,31 @@ int main()
BB[i] = i;
}
for (int i = 0; i < N; ++i) {
for (int j = 0; j < N; ++j) {
for (MAT_SIZE_T i = 0; i < N; ++i) {
for (MAT_SIZE_T j = 0; j < N; ++j) {
std::cout << " " << AA[i + j*N];
}
std::cout << " " << std::endl;
}
std::cout << std::endl;
int kl = 1;
int ku = 1;
int nrowAB = 2*kl + ku + 1;
int ldab = nrowAB;
MAT_SIZE_T kl = 1;
MAT_SIZE_T ku = 1;
MAT_SIZE_T nrowAB = 2*kl + ku + 1;
MAT_SIZE_T ldab = nrowAB;
std::vector<double> AB(nrowAB*N, 0.);
BandMatrixCoeff bmc(N, ku, kl);
// Rewrite AA matrix in band format AB
for (int i = 0; i < N; ++i) {
for (int j = -1; j < 2; ++j) {
for (MAT_SIZE_T i = 0; i < N; ++i) {
for (MAT_SIZE_T j = -1; j < 2; ++j) {
if (i + j > -1 && i + j < N)
AB[bmc(i, i + j)] = AA[i + N*(i + j)];
}
}
for (int i = 0; i < nrowAB; ++i) {
for (int j = 0; j < N; ++j) {
for (MAT_SIZE_T i = 0; i < nrowAB; ++i) {
for (MAT_SIZE_T j = 0; j < N; ++j) {
std::cout << " " << AB[i + j*nrowAB];
}
std::cout << " " << std::endl;
@ -113,7 +115,7 @@ int main()
dgesv_(&N ,&nrhs, &AA[0], &lda, &ipiv[0], &BB[0], &ldb, &info);
if (info == 0) {
for (int i = 0; i < N; ++i) {
for (MAT_SIZE_T i = 0; i < N; ++i) {
std::cout << BB[i] << ' ';
}
std::cout << std::endl;
@ -121,14 +123,14 @@ int main()
std::cerr << "Something went wrong in debsv_()\n";
}
for (int i = 0; i < N; ++i) {
for (MAT_SIZE_T i = 0; i < N; ++i) {
BB[i] = i;
}
dgbsv_(&N, &kl, &ku, &nrhs, &AB[0], &ldab, &ipiv[0], &BB[0], &ldb, &info);
if (info == 0) {
for (int i = 0; i < N; ++i) {
for (MAT_SIZE_T i = 0; i < N; ++i) {
std::cout << BB[i] << ' ';
}
std::cout << std::endl;