matrix inverse - dgetri

[smtp12357]

Hi everybody,

I can not find mistake with dgetri. Please see my code:

[cpp]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, lworkspace = -1;
int *ipiv = new int [SIZE];
dgetrf(&SIZE, &SIZE, a, &SIZE, ipiv, &info);

double workspace;
dgetri(&SIZE, a, &SIZE, ipiv, &workspace, &lworkspace, &info);[/cpp]

true answer a=
| 3.75 -6.5 2.25 |
| -5 9 -3 |
| 1.5 -3 1 |

but I have a=
| 12 12 10 |
| 0 -2 -6 |
| 0.(3) 0.5 -0.(3) |

both after dgetrf and dgetri. What's wrong?

[Gennady_F_Intel]

Hi,
I 've slightly changed your code:

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

===
3.750000 -6.500000 2.250000
-5.000000 9.000000 -3.000000
1.500000 -3.000000 1.000000
Press any key to continue . . .

- Gennady

[smtp12357]

Gennady,

Thank you for the help!