LU factorization problem in C/C++

[mpbl]

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.
icc main.cpp -lmkl_lapack -lmkl -lmkl_intel_thread -liomp5 -lpthread -lm -O0
But when compiling without optimization indication or with "-O1", "-O2", or "-O3"
icc main.cpp -lmkl_lapack -lmkl -lmkl_intel_thread -liomp5 -lpthread -lm
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]

[Gennady_F_Intel]

Hello Matt,

I cannot seelibmkl_core linked with these examples. Please check with Linker Adviser how to select recommended libraries for linking with mkl 10.1.

--Gennady

[mpbl]

Thanks Gennady,

As you suggest I modified the command line as :
icc main.cpp -L/opt/intel/Compiler/11.0/083/mkl/lib/em64t -lmkl_intel_lp64 -lmkl_intel_thread -lmkl_core -openmp -lpthread -lm -I/opt/intel/Compiler/11.0/083/mkl/include -O0

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

Matt.

[TimP]

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?

[mpbl]

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

[mecej4 · Accepted solution]

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.

[mpbl]

Thanks mecej4,

It works and I am going to buy a 1st grade reading book.

Matt.