Stiffness matrix inversion: discrepancy to Matlab

1

Hi community
I am inverting a (positive semi-definite) stiffness matrix with Fortran by using the Doolittle LU method, whic is the the second example that appears on a browsing engine.

The results are very similar to these from the inv command in Matlab, however with increasing dimension, the results start to deviate towards 1% or even more:

Fortran:

Column 1 of ivnverse of stiffness matrix:
3.9852017322009894E-009
-5.4468248058316286E-010
3.5966467990581952E-013
3.8489376774482402E-009
-5.4462950677681609E-010
4.2316496536804665E-013
3.7126770659939833E-009
-5.4447790039172174E-010
5.3906969589759813E-013
3.5764461728524985E-009
-5.4424163688506264E-010
5.7557562377480925E-013
3.4402786639993834E-009
-5.4393418810505790E-010
…
…
…
…

Matlab:

invKs(:,1)

ans =

  4.03129045022659e-09
 -5.37505393356795e-10
 -2.79340670005005e-22
  3.89691410188743e-09
 -5.37505393356825e-10
 -9.73811300825123e-23
  3.76253775354828e-09
 -5.37505393356911e-10
 -3.90319389199924e-22
  3.62816140520909e-09
 -5.37505393357039e-10

…
…
…
…

Is there an explanation to this? Is it advisable to switch to the Lapack routines instead?
Thank you

2

It is unfortunate this code is so high in the ranks of the search engine. The reason you see differences is because the Fortran code does not use pivoting.

I would recommend the LAPACK routine dpstrf instead. Do you need the inverse explicitly, or only the solution to a system of equations?

3

Thanks @ivanpribec for clarifying this issue. I need the inverse explicitly, indeed.
Already tried to link the code to the linear algebra library by:
C:\sandbox>gfortran FEM.f90 -o FEM.exe -L/C:/sandbox -lapack

where this is the directory I have placed the lapack.f and blas.f files, but this returns:
c:/mingw/bin/…/lib/gcc/mingw32/6.3.0/…/…/…/…/mingw32/bin/ld.exe: cannot find -lapack
collect2.exe: error: ld returned 1 exit status

so, something must be missing from the command, right?

4

Not exactly. Unfortunately my Windows PC is down, so I cannot guide you through the process. But you can try to follow the instructions at http://icl.cs.utk.edu/lapack-for-windows/lapack/ Seeing that you use mingw you should look under the section “Prebuilt dynamic libraries using Mingw”.

In case you have experience using CMake, you can try and follow the “Easy Windows Build” path (this will compile the LAPACK library on your computer). Depending on where your LAPACK library gets installed, you might need to modify your PATH environment variable.

Then you should be able to compile your example with:

gfortran FEM.f90 -o FEM.exe -lapack -lblas

Personally, I like to use the Windows Subsystem for Linux which I find is much easier to set up for Fortran development. But I would be happy to learn from other Windows users on our Discourse, how do they get their Fortran development environment set up with Lapack, BLAS, fftw, etc. ? (assuming they are not using Simply Fortran, or Visual Studio + Intel Compiler)

1 Like
5

I use MSYS2 on Windows, which has a native Windows (mingw-w64) OpenBLAS package for LAPACK/BLAS functions. Like WSL, it’s very easy to get set up for development using the package manager:

$ pacman -S mingw-w64-x86_64-openblas