DGESV LAPACK HELP

[Farzad Tatar]

DEAR ALL,
I am trying to use LAPack libraries in Fortran to invert a square matrix. I am new in Fortran and also LAPack. To rely on the results after the inversion, I used a simple 2-2 matrix to check if I used the subroutine correctly. However, I am receiving the wrong results. Here is the code.

real(DP):: TEST(2,2) = TRANSPOSE(reshape((/1,2,3,4/), (/2,2/)))
real(DP):: TEST2(2) = (/1,2/)
PRINT*, 'TEST , TEST2 ARE:', TEST, TEST2
CALL DGESV(2, 1, TEST, 2, IPIV, TEST2, 2, INFO)
PRINT*, 'AFTER THE INVERSION, TEST and TEST2 ARE:', TEST, TEST2

In the output, I get this:

TEST , TEST2 ARE: 1.00000000000000 3.00000000000000
2.00000000000000 4.00000000000000 1.00000000000000
2.00000000000000
AFTER THE INVERSION, TEST , TEST2 ARE: 3.00000000000000
0.333333333333333 4.00000000000000 0.666666666666667
0.000000000000000E+000 0.500000000000000

while it should be this, Test_inv=(/-2, 1, 1.5, -0.5/).

For the following input, however, the result is ok, even though TEST after inversion is still the same.
real(DP):: TEST(2,2) = TRANSPOSE(reshape((/1,0,0,3/), (/2,2/)))
real(DP):: TEST2(2) = (/1,6/)

output:
TEST , TEST2 ARE: 1.00000000000000 0.000000000000000E+000
0.000000000000000E+000 3.00000000000000 1.00000000000000
6.00000000000000
AFTER THE INVERSION, TEST , TEST2 ARE: 1.00000000000000
0.000000000000000E+000 0.000000000000000E+000 3.00000000000000
1.00000000000000 2.00000000000000


Thank you in advance for your help.

Cheers,

Farzad Tatar

[Giorgio Pastore]

Il 31/08/23 16:31, Farzad Tatar ha scritto:

> real(DP):: TEST(2,2) = TRANSPOSE(reshape((/1,2,3,4/), (/2,2/)))
> real(DP):: TEST2(2) = (/1,2/)
> PRINT*, 'TEST , TEST2 ARE:', TEST, TEST2
> CALL DGESV(2, 1, TEST, 2, IPIV, TEST2, 2, INFO)
> PRINT*, 'AFTER THE INVERSION, TEST and TEST2 ARE:', TEST, TEST2

You should be aware that dgesv doesn't invert a matrix. It solves a
linear system where test is the matrix of coefficients, and test2 the
array of the known terms. In aoutput test2 is the array containing the
solution.

Giorgio

[gah4]

On Thursday, August 31, 2023 at 7:31:38 AM UTC-7, Farzad Tatar wrote:
> DEAR ALL,
> I am trying to use LAPack libraries in Fortran to invert a square matrix.
> I am new in Fortran and also LAPack. To rely on the results after the inversion,
> I used a simple 2-2 matrix to check if I used the subroutine correctly.
> However, I am receiving the wrong results. Here is the code.

Looks to me like the answer is right.

It does not return the inverse in the third argument, but the LU
decomposition, not including the diagonal of L, which is all 1.

I am not fast enough doing LU by hand, to tell if it is right.
It does seem the LU for a diagonal matrix is the same, as you see.

But multiplying the original matrix by the result does give the
original values for TEST2.

It is, I believe, easy to get the actual inverse from LU, though.

You don't show the declaration or values for IPIV.
It will sometimes exchange the rows, which you will only know if
you show IPIV.

[gah4]

On Thursday, August 31, 2023 at 7:31:38 AM UTC-7, Farzad Tatar wrote:

> I am trying to use LAPack libraries in Fortran to invert a square matrix.

(snip)

> real(DP):: TEST(2,2) = TRANSPOSE(reshape((/1,2,3,4/), (/2,2/)))
> real(DP):: TEST2(2) = (/1,2/)
> PRINT*, 'TEST , TEST2 ARE:', TEST, TEST2
> CALL DGESV(2, 1, TEST, 2, IPIV, TEST2, 2, INFO)
> PRINT*, 'AFTER THE INVERSION, TEST and TEST2 ARE:', TEST, TEST2

Seems to me that there is an easy way to get the inverse matrix.

Set TEST2 to an identity matrix, and so the solution will be the inverse.

[Farzad Tatar]

>It does not return the inverse in the third argument, but the LU
>decomposition, not including the diagonal of L, which is all 1.

Thanks for telling me this, I had totally ignored this. probably this is the reason. I will check it out.

> Set TEST2 to an identity matrix, and so the solution will be the inverse.

Since the dimensions are not the same, i.e., I don't have two columns of unknowns, this solution will not help.

>You don't show the declaration or values for IPIV.
>It will sometimes exchange the rows, which you will only know if
>you show IPIV.

Here is one example of what I believed was not correct.

A. CODE:

real:: TEST(2,2) = TRANSPOSE(reshape((/1.,2.,3.,4./), (/2,2/)))
real:: TEST2(2) = (1.,1.)
integer :: ipiv2(2)
integer, parameter :: sizet=2
goto 999
...

999 NRHS=1
print*, 'before solving', test, test2
CALL DGESV(sizet, nrhs, TEST, sizet, IPIV2, TEST2, sizet, INFO)
print*, ' after solving:', test, test2
print*, ' ipiv2 vector: ' ipiv2

B. RESULTS IN CMD:
before solving 1.000000 3.000000 2.000000 4.000000
1.000000 1.000000
after solving: 2.000000 4.000000 -1.7014126E+38 1.375000
0.0000000E+00 0.0000000E+00
ipiv2 vector: 2 2
OUTPUTS OF SOLVING LS: INFO 0

C. WHILE I EXPECT THESE VALUES, YOU CAN PUT THEM IN THE EQUATION AND CHECK.
X1=-1, X2=1 BUT THE OUTPUT IS 0 AND 0.

[gah4]

On Friday, September 1, 2023 at 12:37:16 AM UTC-7, Farzad Tatar wrote:

(snip, I wrote)

> > Set TEST2 to an identity matrix, and so the solution will be the inverse.

> Since the dimensions are not the same, i.e., I don't have two columns of unknowns, this solution will not help.

You can have any number of columns of unknowns (greater than zero).
So, with a test2(2,2), you could have two columns and set nrhs to 2.

[gah4]

On Friday, September 1, 2023 at 12:37:16 AM UTC-7, Farzad Tatar wrote:

(snip)

> real:: TEST(2,2) = TRANSPOSE(reshape((/1.,2.,3.,4./), (/2,2/)))
> real:: TEST2(2) = (1.,1.)
> integer :: ipiv2(2)
> integer, parameter :: sizet=2
> goto 999
> ...
>
> 999 NRHS=1
> print*, 'before solving', test, test2
> CALL DGESV(sizet, nrhs, TEST, sizet, IPIV2, TEST2, sizet, INFO)
> print*, ' after solving:', test, test2
> print*, ' ipiv2 vector: ' ipiv2

(snip)

DGESV needs double precision 3rd and 6th argument.

If you give it single precision, it will try to access, in this case,
four double precision variables in an array big enough for two.

The -1.7014126E+38 is the most negative single precision value
on many systems. That is a sign that something is wrong.

If you change to SGESV it should work.