Wrong Results in Matrix Inversion

[rsisnett]

Hi, I'm trying to invert a very large matrix (625 x 625) using MKL.

However I'm getting several problems:

First of all something as simple as (taken from a previous post):

[cpp]#include 
#include 


int main(){
int SIZE = 3;
double *a = new double [SIZE*SIZE];

for(int i = 0; i < SIZE; i++ ) {
for(int j = 0; j < SIZE; j++ ) {
a[j*SIZE+i] = (double) (i+j+1)*(2*i-j);
}
}

int info = 0;
int lworkspace = SIZE;//-1;    !!!!!!
int *ipiv = new int [SIZE];
dgetrf(&SIZE, &SIZE, a, &SIZE, ipiv, &info);

//double workspace;  !!!!! this is array
double* workspace = new double [lworkspace * sizeof(double)];
dgetri(&SIZE, a, &SIZE, ipiv, workspace, &lworkspace, &info);

for(int i = 0; i < SIZE; i++ ) {
   for(int j = 0; j < SIZE; j++ ) {
      printf(" %lf", a[j*SIZE+i] );
   }
   printf("n");
}


}
[/cpp]

keeps giving me segmentation fault.

The real application, on the other hand is giving me wrong results.

Original Matrix

1.51225 -0.488091 0
-0.0795172 1.44725 -0.367735
0 -0.192743 1.44582

PIVOTS

1 0 2 0 3 0

Factored after dgetrf
0.512247 -0.952843 0
-0.0795172 0.371484 -0.989907
0 -0.192743 0.255018

Inverse after dgetri
2.64824 4.484 3.69867
0.73051 4.70592 3.88172
0.315827 2.03455 3.9213

However this doesn't seems right.

Any pointers will be thanked.

Sisnett

[Artem_V_Intel]

Hello Sisnett,

Could you please specify platform and architecture you are working on? Also could you specify link line you used?

I tried to build and run example you submitted on Linux Intel64 and doesn't obtain any issues. Please see the results below:

$ make
icpc -o test main.cpp -I/opt/intel/Compiler/11.1/046/mkl/include -Wl,--start-group /opt/intel/Compiler/11.1/046/mkl/lib/em64t/libmkl_intel_lp64.a /opt/intel/Compiler/11.1/046/mkl/lib/em64t/libmkl_intel_thread.a /opt/intel/Compiler/11.1/046/mkl/lib/em64t/libmkl_core.a -Wl,--end-group -openmp -lpthread
$ ./test
3.750000 -6.500000 2.250000
-5.000000 9.000000 -3.000000
1.500000 -3.000000 1.000000

Also I wrote a simple program that computes the inverse of matrix you submitted:

[cpp]#include 
#include 

int main(int argc, char *argv[])
{
	int		n = 3, ipiv[3], info = 0, lwork = 3, i, j;
	double	a[9] = {1.51225, -0.0795172, 0.0, 
					-0.488091, 1.44725, -0.192743,
					0.0, -0.367735, 1.44582};
	double	work[3];

	for(i = 0; i < n; i++)
	{
		for(j = 0; j < n; j++)
		{
			printf("%15.13ft", a[j*n + i]);
		}
		printf("n");
	}
	printf("n");

	dgetrf(&n, &n, a, &n, ipiv, &info);
	if(info)
	{
		printf("info = %dn", info);
		return 1;
	}

	for(i = 0; i < n; i++)
	{
		for(j = 0; j < n; j++)
		{
			printf("%15.13ft", a[j*n + i]);
		}
		printf("n");
	}
	printf("n");

	dgetri(&n, a, &n, ipiv, work, &lwork, &info);
	if(info)
	{
		printf("info = %dn", info);
		return 1;
	}
	
	for(i = 0; i < n; i++)
	{
		for(j = 0; j < n; j++)
		{
			printf("%15.13ft", a[j*n + i]);
		}
		printf("n");
	}
	return 0;
}
[/cpp]

I built and ran this on Win32 in sequential mode.
The result of the program run:
1.5122500000000 -0.4880910000000 0.0000000000000
-0.0795172000000 1.4472500000000 -0.3677350000000
0.0000000000000 -0.1927430000000 1.4458200000000

1.5122500000000 -0.4880910000000 0.0000000000000
-0.0525820466193 1.4215851762836 -0.3677350000000
0.0000000000000 -0.1355831526774 1.3959613293502

0.6736309983309 0.2351500961549 0.0598089116277
0.0383094079300 0.7285644130096 0.1853056635114
0.0051070466674 0.0971252926759 0.7163522219240

That seems correct.

Also may be this article will be useful.

Thanks,
Art

[yjyincj]

Quoting - rsisnett

Hi, I'm trying to invert a very large matrix (625 x 625) using MKL.

However I'm getting several problems:

First of all something as simple as (taken from a previous post):

[cpp]#include 
#include 


int main(){
int SIZE = 3;
double *a = new double [SIZE*SIZE];

for(int i = 0; i < SIZE; i++ ) {
for(int j = 0; j < SIZE; j++ ) {
a[j*SIZE+i] = (double) (i+j+1)*(2*i-j);
}
}

int info = 0;
int lworkspace = SIZE;//-1;    !!!!!!
int *ipiv = new int [SIZE];
dgetrf(&SIZE, &SIZE, a, &SIZE, ipiv, &info);

//double workspace;  !!!!! this is array
double* workspace = new double [lworkspace * sizeof(double)];
dgetri(&SIZE, a, &SIZE, ipiv, workspace, &lworkspace, &info);

for(int i = 0; i < SIZE; i++ ) {
   for(int j = 0; j < SIZE; j++ ) {
      printf(" %lf", a[j*SIZE+i] );
   }
   printf("n");
}


}
[/cpp]

keeps giving me segmentation fault.

The real application, on the other hand is giving me wrong results.

Original Matrix

1.51225 -0.488091 0
-0.0795172 1.44725 -0.367735
0 -0.192743 1.44582

PIVOTS

1 0 2 0 3 0

Factored after dgetrf
0.512247 -0.952843 0
-0.0795172 0.371484 -0.989907
0 -0.192743 0.255018

Inverse after dgetri
2.64824 4.484 3.69867
0.73051 4.70592 3.88172
0.315827 2.03455 3.9213

However this doesn't seems right.

Any pointers will be thanked.

Sisnett

Hi there. I am just curious how do you get your codes highlighted? I found the website of syntaxhighlighter but have no idea how to use it. Thanks.

[Artem_V_Intel]

Hello,

There is a special button in editor for this:

Thanks,

Art

[rsisnett]

My compile Script is:

icc $1 -o $2 -g -w -L$MKLPATH $MKLPATH/libmkl_solver_ilp64_sequential.a -Wl,--start-group -lmkl_intel_ilp64 -lmkl_sequential -lmkl_core -Wl,--end-group -lpthread -I/share/apps/intel/ict/impi/3.2.1.009/include64 -L/share/apps/intel/ict/impi/3.2.1.009/lib64 -lmpigc4 -Xlinker --enable-new-dtags -Xlinker -rpath -Xlinker $libdir -Xlinker -rpath -Xlinker -lmpi -lmpigf -lmpigi -lrt -lpthread -ldl

(I'm also using MPI so half of the script links to mpi libraries)

The environment I'm running in:

Red Hat Linux Enterprise Edition 64-bits Architecture

Xeon Cuad Core Processors

MKL 10.2.0.013

[Artem_V_Intel]

Hello Sisnett,

First of all, if you plan to use cluster version of MKL LAPACK (ScaLAPACK) your link line should be like this:

$MKLPATH/libmkl_scalapack_ilp64.a $MKLPATH/libmkl_solver_ilp64_sequential.a -Wl,--start-group $MKLPATH/libmkl_intel_ilp64.a $MKLPATH/libmkl_sequential.a $MKLPATH/libmkl_core.a $MKLPATH/libmkl_blacs_intelmpi_ilp64.a -Wl,--end-group -lpthread

Here MKLPATH is a full path to MKL lib/em64t directory.

Please try to use this tool to adjust correct linking parameters:
http://software.intel.com/en-us/articles/intel-mkl-link-line-advisor/

As I can see from your post you tried to use ILP64 model. In this case you should add -DMKL_ILP64 compiler option for ILP64 interface support in MKL. Also you should use MKL_INT type instead of int.

Thanks,
Art