How to calculate the inverse of a matrix using lapack

[JReinoso]

Good afternoon,

I am wondering how to calculate the inverse of a matrix of size n using the lapack library and implement it into my user element subroutine.
I want to use the library lapack as external which has an utility to calculate the inverse in a very robust way.
This motivation came because I am passing a user element subroutine from ABAQUS to FEAP. In ABAQUS it works without problems (inter compiler) but when passing to FEAP I have problems when calculating the inverse of a matrix of size 7x7, getting: "Program received signal SIGILL: Illegal instruction".
Attached I send the subroutine that I have created, in which I want to introduce the EAS method to my element.
Any routes of how to achieve this would be so much appreciated.


Thanks in advance

[Prof. S. Govindjee]

Do you really want to invert a matrix or solve linear equations?  If you want to invert and do not care too much about the method you can try calling program/invert.f
See the comments in the file on how to call.  If you are solving equations you can call the LAPACK routine DGESV.  The most basic call is

Code: [Select]

Call DGSEV(n,1,A,n,ipiv,rhs,n,info)where n is the size of A (n x n); rhs is the right-hand size and is replaced by the solution on exit. See https://netlib.org/lapack/explore-html-3.6.1/d7/d3b/group__double_g_esolve_ga5ee879032a8365897c3ba91e3dc8d512.html for details.

[JReinoso]

Many thanks for your reply. I saw the routine in /program and it says that only works for small matrices.
I need to invert square matrix for the application of the enhanced assumed strain method
Actually the code for that in the numerical recipes of Fortran was working for feap8.4
Then,,I wonder whether the code from the following link can be used
https://www.netlib.org/lapack/explore-html-3.6.1/dd/d9a/group__double_g_ecomputational_ga56d9c860ce4ce42ded7f914fdb0683ff.html

I mean, shall I copy the corresponding folder into my user elemt and call dgetri function?
Thanks

[Prof. S. Govindjee]

There is no need to copy lapack routines into your working folder.  You should just link to the LAPACK libraries within makefile.in by adding -lblas and llapack to LDOPTIONS (I am pretty sure that BLAS is needed for LAPACK to work).

[JReinoso]

Thanks again for your prompt and helpful reply.
I just call the corresponding routines and there is no compilation error

I am running a Ubuntu machine with gfortran. I enclosed you the makefile.in
The question is that I am trying to replicate in FEAP the benchmark example of calculating the inverse found in
https://support.nag.com/numeric/nl/nagdoc_26/nagdoc_fl26/pdf/f07/f07ajf.pdf

for doing so, I am creating a dummy user element where I create the entrey matrix aaux(4,4)
when I try to compute the inverse via


Call dgetrf(n,n,a,lda,ipiv,info)
Call dgetri(n,a,lda,ipiv,work,lwork,info)


I GET a complete different result in comparison with the becnhmark.
I am not able to see where I made the mistake
Enclosed you have the makefile.in, the I1tension (input) and the user element (elm01.f)

[Prof. S. Govindjee]

This is very curious.

If you add a test to your code to multiply the computed inverse against the original AAUX, then you get back the identity! (see the attached modified elmt01.f)

Code: [Select]

    Matrix: AAUX
    row/col       1           2           3           4
       1  1.8000E+00  2.8800E+00  2.0500E+00 -8.9000E-01
       2  5.2500E+00 -2.9500E+00 -9.5000E-01 -3.8000E+00
       3  1.5800E+00 -2.6900E+00 -2.9000E+00 -1.0400E+00
       4 -1.1100E+00 -6.6000E-01 -5.9000E-01 -8.0000E-01
 info after LU           0

    Matrix: AAUX-Inv
    row/col       1           2           3           4
       1  5.2823E-02  4.8454E-02  1.0911E-01 -4.3076E-01
       2  4.9394E-01 -2.5811E-01  4.0255E-01  1.5320E-01
       3 -3.4937E-01  2.9036E-01 -6.5994E-01 -1.3261E-01
       4 -2.2313E-01 -6.8430E-02  3.2162E-03 -6.8091E-01
 info after inv           0

    Matrix: Idtest-1
    row/col       1           2           3           4
       1  1.0000E+00  1.5959E-16 -2.8449E-16  0.0000E+00
       2 -1.1102E-16  1.0000E+00 -6.2450E-17  4.4409E-16
       3  1.3878E-16 -1.3878E-16  1.0000E+00  0.0000E+00
       4  1.1102E-16 -6.9389E-17  3.3393E-17  1.0000E+00

    Matrix: Idtest-2
    row/col       1           2           3           4
       1  1.0000E+00 -5.5511E-17  0.0000E+00  5.5511E-17
       2 -5.5511E-17  1.0000E+00  2.7756E-17  4.1633E-17
       3 -1.1102E-16 -1.3878E-17  1.0000E+00 -2.2204E-16
       4  1.1102E-16 -1.1102E-16  0.0000E+00  1.0000E+00


The two tests are AAUX*AAUXINV and AAUXINV*AAUX.  Thus the computed inverse is correct. 

However the inverse as shown in the NAG support page also seems to be correct; I tested it in python with numpy and scipy.

Your matrix AAUX is not singular; the determinant is roughly 4.06.  I am not sure what is going on.

I also used FEAP's invert and it matches what your code is computing and checking that it is the inverse gives the identity.  I also punched in the NAG support page inverse and computed the identity check but did not receive the identity; however I have only used the given digits.

 

Code: [Select]

    Matrix: FEAP-invert
    row/col       1           2           3           4
       1  5.2823E-02  4.8454E-02  1.0911E-01 -4.3076E-01
       2  4.9394E-01 -2.5811E-01  4.0255E-01  1.5320E-01
       3 -3.4937E-01  2.9036E-01 -6.5994E-01 -1.3261E-01
       4 -2.2313E-01 -6.8430E-02  3.2162E-03 -6.8091E-01

    Matrix: Idtest-3(feap)
    row/col       1           2           3           4
       1  1.0000E+00  5.5511E-17  5.5511E-17  5.5511E-17
       2  2.2204E-16  1.0000E+00  5.5511E-17 -1.3878E-16
       3 -3.0531E-16  3.7470E-16  1.0000E+00  5.5511E-17
       4  1.1102E-16  5.5511E-17  1.1102E-16  1.0000E+00

    Matrix: NAGinv
    row/col       1           2           3           4
       1  1.7720E+00  5.7569E-01  8.4325E-02  4.8155E+00
       2 -1.1747E-01 -4.4562E-01  4.1136E-01 -1.7126E+00
       3  1.7986E-01  4.5266E-01 -6.6757E-01  1.4824E+00
       4  2.4944E+00  7.6498E-01 -3.5954E-02  7.6119E+00

    Matrix: Idtest-4(nag)
    row/col       1           2           3           4
       1  1.0000E+00  7.7000E-09  4.1000E-09 -7.7048E+00
       2  2.7600E-08  1.0000E+00 -4.3000E-09  2.7401E+00
       3 -1.7500E-08  9.2000E-09  1.0000E+00 -2.3718E+00
       4  1.8600E-08 -1.9700E-08 -1.4600E-08 -1.1179E+01

[Prof. S. Govindjee]

Here is the modified elmt01.f

[Prof. S. Govindjee]

I also checked the condition number of your AAUX and it is 78.62.  So there should not be any sensitivities of the form you are seeing.

[Prof. S. Govindjee]

Ignore everything.  All is fine.  When you punched in AAUX you had a sign error in the (4,4) entry.  Everything works now:

Code: [Select]

    Matrix: AAUX
    row/col       1           2           3           4
       1  1.8000E+00  2.8800E+00  2.0500E+00 -8.9000E-01
       2  5.2500E+00 -2.9500E+00 -9.5000E-01 -3.8000E+00
       3  1.5800E+00 -2.6900E+00 -2.9000E+00 -1.0400E+00
       4 -1.1100E+00 -6.6000E-01 -5.9000E-01  8.0000E-01
 info after LU           0

    Matrix: AAUX-Inv
    row/col       1           2           3           4
       1  1.7720E+00  5.7569E-01  8.4325E-02  4.8155E+00
       2 -1.1747E-01 -4.4562E-01  4.1136E-01 -1.7126E+00
       3  1.7986E-01  4.5266E-01 -6.6757E-01  1.4824E+00
       4  2.4944E+00  7.6498E-01 -3.5954E-02  7.6119E+00
 info after inv           0

    Matrix: Idtest-1
    row/col       1           2           3           4
       1  1.0000E+00  1.1102E-16 -3.8164E-16 -8.8818E-16
       2  0.0000E+00  1.0000E+00 -1.1102E-16 -3.5527E-15
       3  0.0000E+00 -3.3307E-16  1.0000E+00 -1.7764E-15
       4  0.0000E+00 -1.1102E-16  1.1449E-16  1.0000E+00

    Matrix: Idtest-2
    row/col       1           2           3           4
       1  1.0000E+00 -8.8818E-16  0.0000E+00  0.0000E+00
       2  0.0000E+00  1.0000E+00  0.0000E+00  0.0000E+00
       3  0.0000E+00  0.0000E+00  1.0000E+00  0.0000E+00
       4  0.0000E+00  0.0000E+00  8.8818E-16  1.0000E+00

    Matrix: FEAP-invert
    row/col       1           2           3           4
       1  1.7720E+00  5.7569E-01  8.4325E-02  4.8155E+00
       2 -1.1747E-01 -4.4562E-01  4.1136E-01 -1.7126E+00
       3  1.7986E-01  4.5266E-01 -6.6757E-01  1.4824E+00
       4  2.4944E+00  7.6498E-01 -3.5954E-02  7.6119E+00

    Matrix: Idtest-3(feap)
    row/col       1           2           3           4
       1  1.0000E+00 -4.4409E-16  0.0000E+00  0.0000E+00
       2  6.6613E-16  1.0000E+00  0.0000E+00 -4.4409E-16
       3 -4.4409E-16  2.2204E-16  1.0000E+00  0.0000E+00
       4 -1.7764E-15 -8.8818E-16 -8.8818E-16  1.0000E+00

    Matrix: NAGinv
    row/col       1           2           3           4
       1  1.7720E+00  5.7569E-01  8.4325E-02  4.8155E+00
       2 -1.1747E-01 -4.4562E-01  4.1136E-01 -1.7126E+00
       3  1.7986E-01  4.5266E-01 -6.6757E-01  1.4824E+00
       4  2.4944E+00  7.6498E-01 -3.5954E-02  7.6119E+00

    Matrix: Idtest-4(nag)
    row/col       1           2           3           4
       1  1.0000E+00  7.7000E-09  4.1000E-09  1.5900E-08
       2  2.7600E-08  1.0000E+00 -4.3000E-09 -1.8100E-08
       3 -1.7500E-08  9.2000E-09  1.0000E+00  1.2900E-08
       4  1.8600E-08 -1.9700E-08 -1.4600E-08  1.0000E+00


[JReinoso]

Thanks a lot and very sorry for bothering with this issue!!
Now I believe we can finalize the enhanced assumed strain implementation with fracture and plasticity