Wrong Results in Matrix Inversion

Wrong Results in Matrix Inversion

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):

#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");
}


}

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

6 posts / 0 new
Last post
For more complete information about compiler optimizations, see our Optimization Notice.

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:

#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;
}

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

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):

#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");
}


}

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.

Hello,
There is a special button in editor for this:

Thanks,
Art

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

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

Leave a Comment

Please sign in to add a comment. Not a member? Join today