Intel® oneAPI Math Kernel Library
Ask questions and share information with other developers who use Intel® Math Kernel Library.

## LU factorization problem in C/C++

Beginner
1,442 Views
Hello,
I messed around for a while now and I can't find the solution to a strange bug that I have.
I use the 10.1.1 version of the mkl and the 11.0 version of icc. I wrote a mock-up for LU factorization (?gttrf / ? gttrs).
When compiling with the "-O0" option, all works fine.
But when compiling without optimization indication or with "-O1", "-O2", or "-O3"
I've got a segfault. Did I missed something ?
I've attached the code below.

Thanks,

Matt.

[cpp]#include
#include
#include
#include

#include "mkl_types.h"
#include "mkl_lapack.h"

int main(int argc, char **argv) {
MKL_INT n = 5;
double *dl  = new double[4];
double *d   = new double[5];
double *du  = new double[4];
double *du2 = new double[3];
MKL_INT *ipiv  = new MKL_INT[5];
double *t   = new double[5];
double tmin, tmax;

double alpha = 1.e-6f;
double dt    = 0.001f;
double dx    = 1.0e-2f;
double lambda = dt/(dx*dx);

for (int i = 0; i < 4; ++i) {
dl = -alpha*lambda;
du = -alpha*lambda;
}
for (int i = 0; i < 5; ++i) {
d = 1.f + 2*alpha*lambda;
}

for (int i = 0; i < 5; ++i) {
t = 1.f;
}

MKL_INT info;
MKL_INT ldb = 5;
char trans='N';
MKL_INT nrhs = 5;
double error = 1.f;

// LU factorization of a tridiagonal matrix
dgttrf(&n, dl, d, du, du2, ipiv, &info);
if (info != 0) {
std::cerr << "Error in the LU factorization." << std::endl;
return -1;
}

for (int i = 0; i < 5; ++i) {
std::cout << t << " ";
}
std::cout << std::endl;
// Solve the tridiag LU system
dgttrs(&trans, &n, &nrhs, dl, d, du, du2, ipiv, t, &ldb, &info);
if (info != 0) {
std::cerr << "Error in the LU solver." << std::endl;
return -1;
}

for (int i = 0; i < 5; ++i) {
std::cout << t << " ";
}
std::cout << std::endl;
std::cout << std::endl;

delete [] dl;
delete [] d;
delete [] du;
delete [] du2;
delete [] ipiv;
delete [] t;

return 0;
}
[/cpp]
1 Solution
Honored Contributor III
1,442 Views
The seg-fault is caused by this bug in your program: you are solving with a single RHS vector, yet you give

MKL_INTnrhs=5;

When you make the call to dgttrs, the library routine will try to use 5 X 5 matrices for X and B in

A X = B

instead of, as you probably intended, 5 X 1 vectors x and b such that

A x = b

Change the value of nrhs from 5 to 1.

6 Replies
Moderator
1,442 Views
Hello Matt,
Beginner
1,442 Views

As you suggest I modified the command line as :

but the problem remain the same. It segfaults when not using "-O0".

Matt.
Honored Contributor III
1,442 Views
Did you check your stack settings and other usual suspects, such as data overruns or uninitialized data? Why would you use an older MKL installation than the one provided with your compiler?
Beginner
1,442 Views
Thanks TimP,

The stack size looks ok.
stack size (kbytes, -s) unlimited

\$ versionquery/_results/intel_lp64_parallel_em64t_so/getversionstring_c.out
Intel Math Kernel Library Version 10.1.1 Product Build 082212.12 for Intel 64 architecture applications

That's the version I found in: /opt/intel/Compiler/11.0/083/mkl

I changed it for the version in /opt/intel/mkl/10.1.2.024/lib/em64t, but it gives the same error.

I tried to use valgrind memcheck on the output. For this, i changed the mkl to sequential:
icc main.cpp -L/opt/intel/mkl/10.1.2.024/lib/em64t -lmkl_intel_lp64 -lmkl_sequential -lmkl_core -lpthread -lm -I/opt/intel/mkl/10.1.2.024/include -O0

when the following line is commented (only dgttrf is executed) (line 54)
dgttrs(&trans, &n, &nrhs, dl, d, du, du2, ipiv, t, &ldb, &info);
valgrind says there is 0 errors.

when using dgttrs in addition with the "-O0" option, memcheck gives 52 errors but no conditionnal jump depends on ...().
Address 0x6704280 is not stack'd, malloc'd or (recently) free'd
where the dgttrs occures.

Is this observation relevant ?

Matt

Honored Contributor III
1,443 Views
The seg-fault is caused by this bug in your program: you are solving with a single RHS vector, yet you give

MKL_INTnrhs=5;

When you make the call to dgttrs, the library routine will try to use 5 X 5 matrices for X and B in

A X = B

instead of, as you probably intended, 5 X 1 vectors x and b such that

A x = b

Change the value of nrhs from 5 to 1.

Beginner
1,442 Views
Thanks mecej4,