Community
cancel
Showing results for 
Search instead for 
Did you mean: 
mpbl
Beginner
141 Views

LU factorization problem in C/C++

Jump to solution
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]
0 Kudos
1 Solution
mecej4
Black Belt
141 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.

View solution in original post

6 Replies
Gennady_F_Intel
Moderator
141 Views
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
Beginner
141 Views
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
Black Belt
141 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?
mpbl
Beginner
141 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



mecej4
Black Belt
142 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.

View solution in original post

mpbl
Beginner
141 Views
Thanks mecej4,

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

Matt.
Reply