topic LU factorization problem in C/C++ in Intel® oneAPI Math Kernel Library

LU factorization problem in C/C++

mpbl — Fri, 03 Feb 2012 11:28:23 GMT

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]

LU factorization problem in C/C++

Gennady_F_Intel — Fri, 03 Feb 2012 12:29:54 GMT

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

LU factorization problem in C/C++

mpbl — Fri, 03 Feb 2012 12:45:49 GMT

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.

LU factorization problem in C/C++

TimP — Fri, 03 Feb 2012 15:40:07 GMT

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?

LU factorization problem in C/C++

mpbl — Fri, 03 Feb 2012 16:10:30 GMT

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

LU factorization problem in C/C++

mecej4 — Fri, 03 Feb 2012 18:29:38 GMT

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.

LU factorization problem in C/C++

mpbl — Sun, 05 Feb 2012 10:33:39 GMT

Thanks mecej4,

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

Matt.