Calculating Inverse of a Large Matrix in FORTRAN 90

[mirzay...@gmail.com]

Hello Everyone !

First of all, I am glad to see this group for help in FORTRAN. And I hope to get valuable suggestions and help here.
I am working in FORTRAN 90 where I need to calculate Inverse of a 7x7 Matrix accurately. I have used a subroutine in main program which augmenting the matrix with Identity matrix and then perform calculations for Inverse determination by pivoting. But the problem is that most of the time I do not get the correct answer. I verify my inverse using fact that : matrix A multiplied by its Inverse = Identity Matrix . But it do not gives Identity matrix when I use the Inverse calculated by the subroutine. The code is attached at end (named: MatInv.f90) .
I am working in Ubuntu 16.04 LTS. I have installed FORTRAN 90 and I compile my programs by command " gfortran porgram.f90 -o program" in terminal.

I searched on INTERNET that LAPACK is a library which can efficiently calculates inverse of large matrices. So following instruction in a video on YOUTUBE, I installed LAPACK by commands in terminal as follows:

1) tar zxvf lapack-3.6.0.gz
2) cd lapack-3.6.0

◇ gfortran
cp make.inc.example make.inc

4) make blaslib
5) make lapacklib
6) sudo ln -s $HOME/lapack-3.6.0/librefblas.a /usr/local/lib/libblas.a
7) sudo ln -s $HOME/lapack-3.6.0/liblapack.a /usr/local/lib/liblapack.a

(on checking in /usr/local/bin these libraries are present there)
Now I want to use the program which uses LAPACK to find inverse of a large matrix but I do not know how to compile the code using these libraries. After searching on INTERNET and using file matrix_inverse.f90 (attached at end named: inverse_mat.f90) when I use command :
gfortran my_program.f90 -llapack -lblas
It says:
(.text+0x20): undefined reference to `main'collect2: error: ld returned 1 exit status

I am completely stuck here. I do not know certain things as follows:

(I). Whether I need to use LAPACK or not to find Inverse of a 7x7 matrix
(II). Why my program (without LAPACK) do not give correct result for inverse. If the inverse is accurate it should satisfy condition:
matrix A multiplied by its Inverse = Identity Matrix
but it do not satisfy this condition.
(III). If I use LAPACK, then how to accuratly use it to write and compile a program for inverse
(IV). What is the reason of the error mentioned above
(V). How to check whether I have installed LAPACK correctly or not
(VI). How I can install and use LAPACK, step by step instruction.
(VII). Finally, without knowing answers of all above questions, how can I find ACCURATE Inverse of a 7x7 martrix in FORTRAN 90?

I shall be highly grateful for any suggestions and help.
Best Regards,
Mirza Younis Baig

####################################################################

File:
{code}
SUBROUTINE MATINV(MATSUM,INVERSE)
IMPLICIT NONE
Double Precision, INTENT(IN),DIMENSION(100,100) :: MATSUM
Double Precision, INTENT(OUT),DIMENSION(100,100) :: INVERSE
INTEGER :: I,J,NROWS,NCOLS,M,N,K,p,q
REAL :: summ,xmult,xmult2,xmult3
Double Precision, DIMENSION(100,100) :: A, Ainv, Ident

NROWS=7
NCOLS=7
Do I=1,NROWS
DO J=1,NCOLS
A(I,J)=MATSUM(I,J)
END Do
END DO
PRINT*,''
PRINT*,''
PRINT*,' ---------------------------------------------------'
!Echoing Input Matrix
PRINT*,'>>>>>>>>Input Matrix<<<<<<<<'
DO I=1,NROWS
Print '(7(1x,F9.4))', (A(I,J),J=1,NCOLS)
END DO
PRINT*,' ---------------------------------------------------'
!Defining Identity Matrix
Do I=1,NROWS
Do J=1,NCOLS
IF(I.EQ.J) THEN
Ident(I,J)=1.0
ELSE
Ident(I,J)=0.0
END IF
END DO
END DO
!Row Operations
! PRINT*,'>>>>>>>>Row Operations<<<<<<<<'
DO I=1,NROWS
DO J=1,NCOLS
IF(I.EQ.J) THEN
xmult=A(I,J)
DO K=1,NCOLS
A(I,K)=A(I,K)/xmult
Ident(I,K)=Ident(I,K)/xmult
END DO
DO M=I+1,NROWS
xmult2=A(M,J)
DO N=1,NCOLS
A(M,N)=A(M,N)-xmult2*A(I,N)
Ident(M,N)=Ident(M,N)-xmult2*Ident(I,N)
END DO
END DO
IF(I.GT.1) then
DO p=I-1,1,-1
xmult3=A(p,J)
Do q=1,NCOLS
A(p,q)=A(p,q)-xmult3*A(I,q)
Ident(p,q)=Ident(p,q)-xmult3*Ident(I,q)
END DO
END DO
END IF
END IF

END DO
! PRINT*,I,'>>>>'
! PRINT*,''
! DO K=1,NROWS
! Print '(12(1x,F10.4))',(A(K,J),J=1,NCOLS)
! END DO
! PRINT*,'======================iii=========================='
! DO K=1,NROWS
! Print '(12(1x,F10.4))',(Ident(K,J),J=1,NCOLS)
! END DO

! PRINT*,'==================================================='
END DO



! PRINT*,'>>>>>>>>Matrix Inverse<<<<<<<<'
! DO I=1,NROWS
! Print '(10(1x,F10.4))',(Ident(I,J),J=1,NCOLS)
! END DO

! Back Substitution
DO J=1,NCOLS
DO I=NROWS,1,-1
IF(I.EQ.NROWS) THEN
Ainv(I,J)=Ident(I,J)/A(I,I)
ELSE
summ=0.0
DO K=I+1,NROWS
summ=summ+(A(I,K)*Ainv(K,J))
END DO
Ainv(I,J)=(Ident(I,J)-summ)/A(I,I)
END IF
END DO
END DO

Do I=1,NROWS
DO J=1,NCOLS
INVERSE(I,J)=Ainv(I,J)
END Do
END DO


RETURN
END SUBROUTINE MATINV


{/code}

####################################################################

File:inverse_mat.f90
{code}

! Returns the inverse of a matrix calculated by finding the LU
! decomposition. Depends on LAPACK.
function inv(A) result(Ainv)
real(dp), dimension(:,:), intent(in) :: A
real(dp), dimension(size(A,1),size(A,2)) :: Ainv

real(dp), dimension(size(A,1)) :: work ! work array for LAPACK
integer, dimension(size(A,1)) :: ipiv ! pivot indices
integer :: n, info

! External procedures defined in LAPACK
external DGETRF
external DGETRI

! Store A in Ainv to prevent it from being overwritten by LAPACK
Ainv = A
n = size(A,1)

! DGETRF computes an LU factorization of a general M-by-N matrix A
! using partial pivoting with row interchanges.
call DGETRF(n, n, Ainv, n, ipiv, info)

if (info /= 0) then
stop 'Matrix is numerically singular!'
end if

! DGETRI computes the inverse of a matrix using the LU factorization
! computed by DGETRF.
call DGETRI(n, Ainv, n, ipiv, work, n, info)

if (info /= 0) then
stop 'Matrix inversion failed!'
end if
end function inv

{/code}

[Gordon Sande]

The flippant answer is that you need to take a course on numerial analysis.
For numerical linear algebra that would probably be a third course. A good text
will not be very thin and expects you to know a good deal of
mathematics typical of
a third year math student.

A common recomendation to to read "Numerical Recipes in XXX" where XXX
could be Fortran
or a couple other languages. It is a good read and often even humourous
if you know the material. But be warned that MANY experts view it as
seriously flawed to the point of
sugesting it be avoided.

For immediate advice you should understand that 7x7 is usually
considered to be very
small. The other is you did not say be how much the product failed you
expectations.
So 1. find the largest absolute value in your matrix 2. find the
largest absolute value in your calculated inverse and 3. multiply those
two numbers togeher with your roundoff
level. That should be a decent measure of the expected error level.
When you have taken
your numerical analysis course you will see why I suggest this test.
You just calculated two norms and a condition number.

If you were to post the 7x7 matrix I epect quite a few eager folks
would be happy to
calculate both the inverse and how accurate it is reasonalbe to expect
it to be.
Some may even use a system that would provide an exact answer in
rational arithmetic.
(The numerators and denomiantors are likely to be rather long so it is
not actually a
very useful thing for most purposes.)

With a bit more work you might even find the source of the LAPACK
routines you need to
avoid the fuss of installing libraries for just two routines. LAPACK is
a good choice.

> END DO! PRINT*,I,'>>>>'

[herrman...@gmail.com]

On Wednesday, October 19, 2016 at 4:35:49 PM UTC-7, Gordon Sande wrote:

(snip on matrix inversion)

> A common recomendation to to read "Numerical Recipes in XXX"
> where XXX could be Fortran or a couple other languages.
> It is a good read and often even humourous if you know the material.
> But be warned that MANY experts view it as seriously flawed to
> the point of sugesting it be avoided.

Many experts have noted flaws in the programs, not so many
in the explanations. I believe that those experts have not
understood the way the book is to be used.

Use it like a textbook, to learn about the subject.
Don't pretend it is the final answer.

What you are specifically not supposed to do is to
use the programs like black boxes. Use them like
the answers in the back of the book in your college
textbooks.

> For immediate advice you should understand that 7x7
> is usually considered to be very small.

Yes it is. But matrix inverse can fail on even smaller
matrices.

The common problem with homemade matrix inversion
is not doing pivoting. I am sure that NR explains the
reason for, and how pivoting works.

Try your program on:

0 1
1 0

a nice small matrix.

But the first question is, can you use LU decomposition
instead of inversion. (But you still need pivoting.)

[mirzay...@gmail.com]

Dear Gordon,

Thank you for your generous advice. I have read "Numerical Analysis" by Burden & Faires 9th edition and know a little about matrix norms, conditioning and stability. But in my case it is bit difficult to implement the readings of book in programming. Therefore, I would greatly appreciate for help in programming.
As you said 7x7 is a small matrix then why the inverse is always not accurate as in my case? Is it my code fault or any other?
And secondly why it is not always good to find direct inverse of a matrix?
At last, can you refer me some good subroutine which calculate inverse correctly or mention error ,if any , in my program code.
Best,
Mirza

[Gordon Sande]

On 2016-10-20 00:19:34 +0000, mirzay...@gmail.com said:

> Dear Gordon,
>
> Thank you for your generous advice. I have read "Numerical Analysis" by
> Burden & Faires 9th edition and know a little about matrix norms,
> conditioning and stability. But in my case it is bit difficult to
> implement the readings of book in programming. Therefore, I would
> greatly appreciate for help in programming.
> As you said 7x7 is a small matrix then why the inverse is always not
> accurate as in my case? Is it my code fault or any other?

Did you calculate the numbers I suggested?

Otherwise you have just found a new way to say "It doesn't work" and completely
failed to provide any indication of what you expect or how it failed.
Yes you showed all the Fortran but you asked about numbers and there is
not a single number to be found.

> And secondly why it is not always good to find direct inverse of a matrix?

Take that bunercal analysis course. Finding the inverse is equivaent to
using the
worst possible numerical example while LU only is as bad as you actual example.
Technically it is a question of which condition number applies.

> At last, can you refer me some good subroutine which calculate inverse
> correctly or mention error ,if any , in my program code.

LAPACK is OK but building libraries etc is a long hard way to just use
a couple routines.
Just find the source of the desired routines, copy it into your prgram
and compi;e and go.

> Best,
> Mirza

[JB]

On 2016-10-19, mirzay...@gmail.com <mirzay...@gmail.com> wrote:
> Hello Everyone !
>
> First of all, I am glad to see this group for help in FORTRAN. And I hope to get valuable suggestions and help here.
> I am working in FORTRAN 90 where I need to calculate Inverse of a 7x7 Matrix accurately. I have used a subroutine in main program which augmenting the matrix with Identity matrix and then perform calculations for Inverse determination by pivoting. But the problem is that most of the time I do not get the correct answer. I verify my inverse using fact that : matrix A multiplied by its Inverse = Identity Matrix . But it do not gives Identity matrix when I use the Inverse calculated by the subroutine. The code is attached at end (named: MatInv.f90) .
> I am working in Ubuntu 16.04 LTS. I have installed FORTRAN 90 and I compile my programs by command " gfortran porgram.f90 -o program" in terminal.
>
> I searched on INTERNET that LAPACK is a library which can efficiently calculates inverse of large matrices. So following instruction in a video on YOUTUBE, I installed LAPACK by commands in terminal as follows:

Others have already explained why calculating the matrix inverse tends
to be a bad idea, or at least given you pointers so you can learn more
about the problem if you're interested.

Since you're using Ubuntu, instead of downloading and installing
lapack from sources, you can just install it from the package
repositories, like:

sudo apt-get install liblapack3 liblapack-dev libopenblas-dev libopenblas-base

(this also install OpenBLAS, a fairly fast open source BLAS library)
and then compile/link with

gfortran -O2 -g -Wall -fcheck=all my_program.f90 -llapack -lopenblas

(I also enabled optimizations, and added a few debugging options which
can be useful for finding potential bugs in your code)


--
JB

[Louis Krupp]

On Wed, 19 Oct 2016 15:43:13 -0700 (PDT), mirzay...@gmail.com
wrote:

<snip>

>6) sudo ln -s $HOME/lapack-3.6.0/librefblas.a /usr/local/lib/libblas.a
>7) sudo ln -s $HOME/lapack-3.6.0/liblapack.a /usr/local/lib/liblapack.a

I'm not sure this is a good practice. If you want to keep the
libraries under your home directory, you can link to them using the -L
option in your compile or link command. If you want to keep the
libraries in /usr/local/bin, then install them there.

Louis

[mecej4]

mirzay...@gmail.com wrote:

> Hello Everyone !
>
> First of all, I am glad to see this group for help in FORTRAN.
> And I hope to get valuable suggestions and help here. I am
> working in FORTRAN 90 where I need to calculate Inverse of a
> 7x7 Matrix accurately. I have used a subroutine in main
> program which augmenting the matrix with Identity matrix and
> then perform calculations for Inverse determination by
> pivoting. But the problem is that most of the time I do not
> get the correct answer. I verify my inverse using fact that :
> matrix A multiplied by its Inverse = Identity Matrix . But it
> do not gives Identity matrix when I use the Inverse calculated
> by the subroutine. The code is attached at end (named:
> MatInv.f90) . I am working in Ubuntu 16.04 LTS. I have
> installed FORTRAN 90 and I compile my programs by command "
> gfortran porgram.f90 -o program" in terminal.
>
> I searched on INTERNET that LAPACK is a library which can
> efficiently calculates inverse of large matrices. So following
> instruction in a video on YOUTUBE, I installed LAPACK by
> commands in terminal as follows:
>

<-- CUT -->

I took your code, changed your subroutine to a program,
dimensioned the arrays with size 7 X 7, populated the matrix
with a call to RANDOM_NUMBER, fixed up one bug, and ran the
program. The difference A.inv(A) - I contained no element larger
in magnitude than 2D-15.

As rural wisdom tells us, "If it ain't broke much, don't fix it
much".

If all you need is to invert a 7 X 7 matrix, you have done it.
No need for Lapack, etc.

-- mecej4

[campbel...@gmail.com]

On Thursday, October 20, 2016 at 9:26:38 PM UTC+11, mecej4 wrote:
>
>
> As rural wisdom tells us, "If it ain't broke much, don't fix it
> much".
>
> If all you need is to invert a 7 X 7 matrix, you have done it.
> No need for Lapack, etc.
>
> -- mecej4

From what was presented, I could not identify any pivoting to cope with poorly ordered equations. This could be where the problem is for the initial test equations being solved.

[herrman...@gmail.com]

On Thursday, October 20, 2016 at 4:27:18 AM UTC-7, campbel...@gmail.com wrote:

(snip on matrix inversion not working)

> From what was presented, I could not identify any pivoting to cope
> with poorly ordered equations. This could be where the problem is
> for the initial test equations being solved.

I remember in high school days, doing matrix inversion where you
write down the matrix, (in pencil), write an identity matrix to the
right, and then do row operations on the augmented matrix.

Written in pencil, you can erase each element and replace it with
the new value after each operation. The row operations are the
you can multiply a whole row by a constant, or add a multiple
of one row to another row.

The OP mentions pivoting, but it sounds like this is what he
calls the above row operations.

Still in high school, I have the IBM SSP MINV subroutine,
in both Fortran and PL/I. I remember also finding it interesting
that those routines do it in place.

Note that with the pencil and paper method, each time
you change a value on the right side, you set the value
in the matching place on the left to zero or one. The OP
stores both halves. I do remember tracing the code to
figure this out years ago.

Also, it is usual for matrix inversion routines to also
compute the determinant at the same time. It isn't
much extra work, and is useful to know how good the
result is.

I haven't read it recently, but I believe that the
Numerical Recipes books give a good explanation
on how to do matrix inversion, including pivoting and
when it is needed.

As before, try computing the inverse of:

0 1
1 0

with and without pivoting.

[mirzay...@gmail.com]

Dear Gordon,

I am very serious about my problem and I want to know the answer. I do not want to say that it simply doesn't work. I want to find the solution. Here is my
input matrix:
744.5047 909.5175 983.0251 1024.6988 1051.5530 1070.3036 1084.1393
909.5175 1112.7881 1203.5504 1255.0658 1288.2851 1311.4910 1328.6199
983.0251 1203.5504 1302.1218 1358.0994 1394.2078 1419.4376 1438.0635
1024.6988 1255.0658 1358.0994 1416.6292 1454.3913 1480.7802 1500.2639
1051.5530 1288.2851 1394.2078 1454.3913 1493.2255 1520.3664 1540.4069
1070.3036 1311.4910 1419.4376 1480.7802 1520.3664 1548.0349 1568.4664
1084.1393 1328.6199 1438.0635 1500.2639 1540.4069 1568.4664 1589.1877

Inverse Matrix:
194892.17 -2676397.98 12794703.89 -29030326.53 34287822.19 -20440033.21 4870719.46
-2678042.90 37238660.99 -179679469.67 410772414.76 -488277403.17 292695143.71 -70089779.98
12809584.58 -179776156.15 873816705.80 -2009970250.99 2401836380.15 -1446411088.00 347781821.34
-29078476.75 411190368.26 -2010921035.68 4649638868.64 -5581093528.71 3374226665.58 -814158094.21
34360162.28 -488984828.19 2403984004.76 -5583420620.73 6727982437.59 -4081509744.22 987817190.06
-20491621.21 293233933.11 -1448252941.57 3376890289.51 -4083018289.32 2484401012.38 -602897686.84
4884880.78 -70244030.64 348344857.91 -815075119.32 988510499.29 -603096794.27 146707751.59

Now according to your advice, the largest absolute value in Input matrix is 1589.1877 and in Inverse Matrix it is 6727982437.59. The product of these two values is 1.069202694×10¹³.
I apologies if my previous mail was bit out of way. I would appreciate If you guide me further.

Best,
Mirza

[mirzay...@gmail.com]

Hi mecej,

Thank you for your work but there is one problem. That is, the subroutine works for row operations by augmenting Identity Matrix and at last the matrix obtained is not exact inverse of the input matrix. Because if we multiply the matrix with its inverse then we should get Identity matrix. And we do not get Identity matrix after multiplication. Did you understood where is the error or problem?

Best,
Mirza

[Gordon Sande]

> product of these two values is 1.069202694��

And the product with the roundoff is ??
And the observed error when you multiplied them out was ??

The number you did give is a plausible approximation to a condition number.
For a 7x7 matrix it is high so you are going to see numerical problems.

The first 10 second look at the matrix shows it to be highly colinear
so a large
condition number is no surprise. IF your problem is a statistical problem (you
have been silent on the subject matter) then a first sensible step would be to
take the means out of all the observations.

Your problem is going to be more a subject matter problem than a
Fortram problem.
Although getting the Fortran wrong would be a major stumbling block. Using any
standard statistical package would give much more information about
your problem
(and someone else would have fixed the bugs long ago).

[mirzay...@gmail.com]

The product of these two numbers is of the order of E13. And I did not understand what is meant by "product with the round off" and the "error when multiplying". I mean this digit can be rounded off to two decimal places and how come there error?
I can send you the problem statement in private if you don't mind.

[Gordon Sande]

What is the roundoff level of the real arithmetic you are using.
Typically either 8 byte
reals or 16 byte reals (double precision). 10^-6 or 10^-16 are usual values.

> and the "error when multiplying".

Multiply the matrix by its inverse. Should be an identify. What are the
observed errors?

> I mean this digit can be rounded off to two decimal places and how come
> there error?
> I can send you the problem statement in private if you don't mind.

I think you need to find the local statistical advice service assuming
you are at some university. Otherwise you need some serious analytical
advice for whatever institution
you are in.

Using a statistical package will help with whatever advice you choose to get.

[mirzay...@gmail.com]

It is Double precision. How error can occur in Inverse beyond a certain round off level?

[Jos Bergervoet]

On 10/21/2016 6:10 PM, Gordon Sande wrote:
> On 2016-10-21 14:54:36 +0000, mirzay...@gmail.com said:

...
..

>> product of these two values is 1.069202694×10¹³.

>
> And the product with the roundoff is ??
> And the observed error when you multiplied them out was ??
>
> The number you did give is a plausible approximation to a condition number.

Looking at the eigenvalues, I would say cond(Mat) = 238124263.5868,
so that is much more friendly. And if it is correct, then the inverse
can be computed with a correctly pivoting algorithm and using double
precision.

But since the input in Mirza's previous post is having *less accuracy*
than this condition number, we don't know if it's correct!
1) From the given input accuracy, this matrix's condition number is
indistinguishable from infinity.
2) From the given input accuracy, no inverse matrix can be computed.
3) If we could have the input with more accuracy, the (recomputed)
condition number might very well become even larger, so then we
still could not compute an inverse in a meaningful way.
4) It is, however, also possible that with more accurate input the
condition number remains low enough to compute the inverse. You
need to have the more accurate input first.
5) Or you could look at the way the input matrix is defined/computed,
which might indicate whether it is fundamentally singular to begin
with.

--
Jos

[Jos Bergervoet]

Non-pivoting algorithms is in fact not "broken much", in the
sense that a random matrix will hardly ever need pivoting.
Still, in practice, you often get matrices that *do* need it.
For several reasons they are not like random matrices.

Murphy's law surely defies rural wisdom in this case! (And
I'm not even convinced that for raising cattle you can get
around it so easily..)

--
Jos

[dpb]

On 10/21/2016 12:03 PM, Gordon Sande wrote:
> On 2016-10-21 16:35:31 +0000, mirzay...@gmail.com said:

...

...

> Multiply the matrix by its inverse. Should be an identify. What are the
> observed errors?
>
>> I mean this digit can be rounded off to two decimal places and how
>> come there error?
>> I can send you the problem statement in private if you don't mind.
>
> I think you need to find the local statistical advice service assuming
> you are at some university. Otherwise you need some serious analytical
> advice for whatever institution
> you are in.
>
> Using a statistical package will help with whatever advice you choose to
> get.

Just for grins, MATLAB with the above (approximate?) input gives--

>> m= ...
[744.5047 909.5175 983.0251 1024.6988 1051.5530 1070.3036 1084.1393

909.5175 1112.7881 1203.5504 1255.0658 1288.2851 1311.4910 1328.6199
983.0251 1203.5504 1302.1218 1358.0994 1394.2078 1419.4376 1438.0635
1024.6988 1255.0658 1358.0994 1416.6292 1454.3913 1480.7802 1500.2639
1051.5530 1288.2851 1394.2078 1454.3913 1493.2255 1520.3664 1540.4069
1070.3036 1311.4910 1419.4376 1480.7802 1520.3664 1548.0349 1568.4664

1084.1393 1328.6199 1438.0635 1500.2639 1540.4069 1568.4664 1589.1877 ];
>> cond(m)
ans =
5.4404e+08
>> det(m)
ans =
-3.0969e-15
>> inv(m)
ans =
1.0e+04 *
-0.0077 0.1348 -0.2930 0.0472 0.1597 0.0915 -0.1319
0.1348 -0.9685 1.6217 -0.0444 -0.8634 -0.6214 0.7424
-0.2930 1.6217 -2.5395 0.5472 0.8224 0.4896 -0.6549
0.0472 -0.0444 0.5472 -2.0793 1.4658 1.1719 -1.1048
0.1597 -0.8634 0.8224 1.4658 -2.1308 -0.0810 0.6303
0.0915 -0.6214 0.4896 1.1719 -0.0810 -2.8225 1.7719
-0.1319 0.7424 -0.6549 -1.1048 0.6303 1.7719 -1.2548
>> inv(m)*m'
ans =
1.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000
0 1.0000 0 0.0000 0.0000 0.0000 -0.0000
0.0000 0.0000 1.0000 0 0.0000 0.0000 0.0000
0.0000 0.0000 0.0000 1.0000 -0.0000 0 0
-0.0000 0.0000 -0.0000 0.0000 1.0000 0.0000 0.0000
0 0 -0.0000 -0.0000 -0.0000 1.0000 0.0000
0.0000 0 0.0000 0.0000 0.0000 0.0000 1.0000

Matlab uses "specially tweaked" versions of LAPACK but the documentation
for INV() specifically notes--

"... In practice, it is seldom necessary to form the explicit inverse
of a matrix. A frequent misuse of inv arises when solving the system of
linear equations Ax = b. One way to solve this is with x = inv(A)*b.
A better way, from both an execution time and numerical accuracy
standpoint, is to use the matrix division operator x = A\b. This
produces the solution using Gaussian elimination, without forming
the inverse. ..."

As recommended, using inverse is generally _a_bad_idea_(tm)_ and should
be avoided. It does behoove the beginning numerical analyst to become
intimately familiar with the "why's" of the recommendation.

[septc...@gmail.com]

Hi,

I have tried your code (in the first post), and I think you need to declare the work variables
as double precision (in consistent with input/output arrays). Specifically, this line

REAL :: summ,xmult,xmult2,xmult3

should be modified as

REAL*8 :: summ,xmult,xmult2,xmult3
or
double precision :: summ,xmult,xmult2,xmult3

etc (there are also other portable ways to write this, but the point is to use double precision consistently).

-----------
Test:

If I change your routine so that it receives ncols and nrows as arguments

SUBROUTINE MATINV(MATSUM,INVERSE, ncols, nrows )

IMPLICIT NONE
Double Precision, INTENT(IN),DIMENSION(100,100) :: MATSUM
Double Precision, INTENT(OUT),DIMENSION(100,100) :: INVERSE
INTEGER :: I,J,NROWS,NCOLS,M,N,K,p,q

REAL*8 :: summ,xmult,xmult2,xmult3

Double Precision, DIMENSION(100,100) :: A, Ainv, Ident

!-- NROWS=7
!-- NCOLS=7

...

and write the main program as follows,

program main
implicit none
integer, parameter :: n = 7
real*8 A( 100, 100 ), Ainv( 100, 100 ), t( n, n ) !! t = tmp array
integer i, j, itest

!>>>
itest = 2
!<<<

selectcase ( itest )

case ( 1 ) !! simple regular matrix

do i = 1, n
do j = 1, n
t( i, j ) = 1.0d0 / ( abs(i-j) + 1.0d0 )
enddo
enddo

case ( 2 ) !! OP's matrix in double precision

t(1,:) = [ 744.5047d0, 909.5175d0, 983.0251d0, 1024.6988d0, &
1051.5530d0, 1070.3036d0, 1084.1393d0 ]
t(2,:) = [ 909.5175d0, 1112.7881d0, 1203.5504d0, 1255.0658d0, &
1288.2851d0, 1311.4910d0, 1328.6199d0 ]
t(3,:) = [ 983.0251d0, 1203.5504d0, 1302.1218d0, 1358.0994d0, &
1394.2078d0, 1419.4376d0, 1438.0635d0 ]
t(4,:) = [ 1024.6988d0, 1255.0658d0, 1358.0994d0, 1416.6292d0, &
1454.3913d0, 1480.7802d0, 1500.2639d0 ]
t(5,:) = [ 1051.5530d0, 1288.2851d0, 1394.2078d0, 1454.3913d0, &
1493.2255d0, 1520.3664d0, 1540.4069d0 ]
t(6,:) = [ 1070.3036d0, 1311.4910d0, 1419.4376d0, 1480.7802d0, &
1520.3664d0, 1548.0349d0, 1568.4664d0 ]
t(7,:) = [ 1084.1393d0, 1328.6199d0, 1438.0635d0, 1500.2639d0, &
1540.4069d0, 1568.4664d0, 1589.1877d0 ]

case ( 3 ) !! OP's matrix in single precision

t(1,:) = [ 744.5047, 909.5175, 983.0251, 1024.6988, &
1051.5530, 1070.3036, 1084.1393 ]
t(2,:) = [ 909.5175, 1112.7881, 1203.5504, 1255.0658, &
1288.2851, 1311.4910, 1328.6199 ]
t(3,:) = [ 983.0251, 1203.5504, 1302.1218, 1358.0994, &
1394.2078, 1419.4376, 1438.0635 ]
t(4,:) = [ 1024.6988, 1255.0658, 1358.0994, 1416.6292, &
1454.3913, 1480.7802, 1500.2639 ]
t(5,:) = [ 1051.5530, 1288.2851, 1394.2078, 1454.3913, &
1493.2255, 1520.3664, 1540.4069 ]
t(6,:) = [ 1070.3036, 1311.4910, 1419.4376, 1480.7802, &
1520.3664, 1548.0349, 1568.4664 ]
t(7,:) = [ 1084.1393, 1328.6199, 1438.0635, 1500.2639, &
1540.4069, 1568.4664, 1589.1877 ]
endselect

A( 1:n, 1:n ) = t(:,:)

call matinv( A, Ainv, n, n )

print *, "Ainv = "
do i = 1, n
print "(*(f12.4))", Ainv( i, 1:n )
enddo

t = matmul( A( 1:n, 1:n ), Ainv( 1:n, 1:n ) )
print *, "A * Ainv = "
do i = 1, n
print "(*(f12.4))", t( i, 1:n )
enddo

end program

we can get the following results:

-----
For itest = 1 :

---------------------------------------------------
>>>>>>>>Input Matrix<<<<<<<<
1.0000 0.5000 0.3333 0.2500 0.2000 0.1667 0.1429
0.5000 1.0000 0.5000 0.3333 0.2500 0.2000 0.1667
0.3333 0.5000 1.0000 0.5000 0.3333 0.2500 0.2000
0.2500 0.3333 0.5000 1.0000 0.5000 0.3333 0.2500
0.2000 0.2500 0.3333 0.5000 1.0000 0.5000 0.3333
0.1667 0.2000 0.2500 0.3333 0.5000 1.0000 0.5000
0.1429 0.1667 0.2000 0.2500 0.3333 0.5000 1.0000
---------------------------------------------------
Ainv =
1.3601 -0.5884 -0.1056 -0.0546 -0.0364 -0.0287 -0.0350
-0.5884 1.6137 -0.5435 -0.0829 -0.0402 -0.0267 -0.0287
-0.1056 -0.5435 1.6213 -0.5400 -0.0812 -0.0402 -0.0364
-0.0546 -0.0829 -0.5400 1.6225 -0.5400 -0.0829 -0.0546
-0.0364 -0.0402 -0.0812 -0.5400 1.6213 -0.5435 -0.1056
-0.0287 -0.0267 -0.0402 -0.0829 -0.5435 1.6137 -0.5884
-0.0350 -0.0287 -0.0364 -0.0546 -0.1056 -0.5884 1.3601
A * Ainv =
1.0000 0.0000 -0.0000 0.0000 -0.0000 0.0000 0.0000
0.0000 1.0000 -0.0000 0.0000 -0.0000 -0.0000 0.0000
0.0000 0.0000 1.0000 0.0000 -0.0000 0.0000 -0.0000
0.0000 -0.0000 -0.0000 1.0000 0.0000 0.0000 -0.0000
-0.0000 0.0000 -0.0000 -0.0000 1.0000 -0.0000 0.0000
-0.0000 0.0000 -0.0000 -0.0000 -0.0000 1.0000 0.0000
-0.0000 -0.0000 -0.0000 -0.0000 -0.0000 0.0000 1.0000

----
For itest = 2:

---------------------------------------------------
>>>>>>>>Input Matrix<<<<<<<<

744.5047 909.5175 983.0251 1024.6988 1051.5530 1070.3036 1084.1393
909.5175 1112.7881 1203.5504 1255.0658 1288.2851 1311.4910 1328.6199
983.0251 1203.5504 1302.1218 1358.0994 1394.2078 1419.4376 1438.0635
1024.6988 1255.0658 1358.0994 1416.6292 1454.3913 1480.7802 1500.2639
1051.5530 1288.2851 1394.2078 1454.3913 1493.2255 1520.3664 1540.4069
1070.3036 1311.4910 1419.4376 1480.7802 1520.3664 1548.0349 1568.4664
1084.1393 1328.6199 1438.0635 1500.2639 1540.4069 1568.4664 1589.1877

---------------------------------------------------
Ainv =
-76.5341 1348.1077 -2930.3976 471.9861 1596.9489 914.7758 -1319.4799
1348.1077 -9685.4287 16217.1709 -444.2072 -8633.6875 -6213.9708 7423.6791
-2930.3976 16217.1709 -25394.9467 5472.2177 8223.5511 4896.2298 -6548.5709
471.9861 -444.2073 5472.2177 -20792.6764 14657.9507 11719.0732 -11047.5255
1596.9489 -8633.6874 8223.5510 14657.9506 -21308.0561 -810.0949 6302.8895
914.7758 -6213.9709 4896.2299 11719.0733 -810.0951 -28224.5399 17718.8556
-1319.4799 7423.6791 -6548.5710 -11047.5256 6302.8896 17718.8556 -12548.3770

A * Ainv =
1.0000 0.0000 0.0000 -0.0000 0.0000 0.0000 -0.0000
0.0000 1.0000 -0.0000 -0.0000 0.0000 0.0000 -0.0000
0.0000 0.0000 1.0000 -0.0000 0.0000 0.0000 -0.0000
0.0000 0.0000 0.0000 1.0000 0.0000 0.0000 -0.0000
0.0000 0.0000 0.0000 -0.0000 1.0000 -0.0000 -0.0000
-0.0000 0.0000 0.0000 -0.0000 0.0000 1.0000 -0.0000
-0.0000 0.0000 0.0000 -0.0000 0.0000 0.0000 1.0000

-----
For itest = 3:

---------------------------------------------------
>>>>>>>>Input Matrix<<<<<<<<
744.5047 909.5175 983.0251 1024.6989 1051.5530 1070.3036 1084.1393
909.5175 1112.7881 1203.5504 1255.0658 1288.2852 1311.4910 1328.6199

983.0251 1203.5504 1302.1218 1358.0994 1394.2078 1419.4376 1438.0635

1024.6989 1255.0658 1358.0994 1416.6292 1454.3914 1480.7802 1500.2639
1051.5530 1288.2852 1394.2078 1454.3914 1493.2255 1520.3665 1540.4069
1070.3036 1311.4910 1419.4376 1480.7802 1520.3665 1548.0349 1568.4664

1084.1393 1328.6199 1438.0635 1500.2639 1540.4069 1568.4664 1589.1877

---------------------------------------------------
Ainv =
1226.2858 -5381.5101 5434.5522 3357.7945 -4613.9139 -2745.1360 2756.5439
-5381.5101 25499.1781 -27231.1496 -17280.8722 24226.6250 15920.4866 -15887.3775
5434.5522 -27231.1496 32083.1816 14710.2002 -27275.0827 -13230.9528 15635.8056
3357.7945 -17280.8722 14710.2001 24727.0036 -18651.0769 -29212.4668 22412.1919
-4613.9139 24226.6250 -27275.0827 -18651.0769 23751.1226 23633.8067 -21165.6786
-2745.1360 15920.4865 -13230.9528 -29212.4668 23633.8067 30434.5518 -24832.8470
2756.5439 -15887.3775 15635.8056 22412.1919 -21165.6785 -24832.8470 21119.9518

A * Ainv =
1.0000 0.0000 0.0000 0.0000 -0.0000 -0.0000 0.0000
0.0000 1.0000 0.0000 0.0000 -0.0000 -0.0000 0.0000
0.0000 -0.0000 1.0000 0.0000 -0.0000 -0.0000 0.0000
0.0000 -0.0000 0.0000 1.0000 -0.0000 -0.0000 0.0000
0.0000 -0.0000 0.0000 0.0000 1.0000 -0.0000 0.0000
0.0000 -0.0000 0.0000 0.0000 -0.0000 1.0000 0.0000
0.0000 -0.0000 0.0000 0.0000 -0.0000 -0.0000 1.0000

-----

Please note that the input matrix elements also need to be given in double precision.
Indeed, in the above example, "Ainv" changes drastically if you enter the input matrix in
single precision.

-----

To see the reason, I also tried the same calculation using Julia.

# test.jl (note: slow because of the use of global variables)

put = println

A = [

744.5047 909.5175 983.0251 1024.6988 1051.5530 1070.3036 1084.1393 ;
909.5175 1112.7881 1203.5504 1255.0658 1288.2851 1311.4910 1328.6199 ;
983.0251 1203.5504 1302.1218 1358.0994 1394.2078 1419.4376 1438.0635 ;
1024.6988 1255.0658 1358.0994 1416.6292 1454.3913 1480.7802 1500.2639 ;
1051.5530 1288.2851 1394.2078 1454.3913 1493.2255 1520.3664 1540.4069 ;
1070.3036 1311.4910 1419.4376 1480.7802 1520.3664 1548.0349 1568.4664 ;

1084.1393 1328.6199 1438.0635 1500.2639 1540.4069 1568.4664 1589.1877 ]

# A = [ 1.0 / (abs(i-j)+1.0) for i=1:7, j=1:7 ] # for itest = 1

put( "A = " )
display( A )
put()

put( "inv(A) = " )
Ainv = inv( A )
display( Ainv )
put()

@show maxabs( A * Ainv - eye( size(A,1) ) )

put( "A * Ainv =" )
display( A * Ainv )
put()

D, V = eig( A )

put( "eigvals = " )
display( D )
put()

@show cond( A )
@show rank( A )

-----
Run:
$ julia test.jl

We get

A =
7x7 Array{Float64,2}:
744.505 909.518 983.025 1024.7 1051.55 1070.3 1084.14
909.518 1112.79 1203.55 1255.07 1288.29 1311.49 1328.62
983.025 1203.55 1302.12 1358.1 1394.21 1419.44 1438.06
1024.7 1255.07 1358.1 1416.63 1454.39 1480.78 1500.26
1051.55 1288.29 1394.21 1454.39 1493.23 1520.37 1540.41
1070.3 1311.49 1419.44 1480.78 1520.37 1548.03 1568.47
1084.14 1328.62 1438.06 1500.26 1540.41 1568.47 1589.19

inv(A) =
7x7 Array{Float64,2}:
-76.5341 1348.11 -2930.4 471.986 1596.95 914.776 -1319.48
1348.11 -9685.43 16217.2 -444.207 -8633.69 -6213.97 7423.68
-2930.4 16217.2 -25394.9 5472.22 8223.55 4896.23 -6548.57
471.986 -444.207 5472.22 -20792.7 14658.0 11719.1 -11047.5
1596.95 -8633.69 8223.55 14658.0 -21308.1 -810.095 6302.89
914.776 -6213.97 4896.23 11719.1 -810.095 -28224.5 17718.9
-1319.48 7423.68 -6548.57 -11047.5 6302.89 17718.9 -12548.4

maxabs(A * Ainv - eye(size(A,1))) = 1.4901161193847656e-8

A * Ainv =
7x7 Array{Float64,2}:
1.0 2.79397e-9 -5.58794e-9 … 7.45058e-9 0.0 -1.86265e-9
0.0 1.0 -1.86265e-9 1.02445e-8 -3.72529e-9 -1.86265e-9
-1.16415e-9 9.31323e-9 1.0 1.30385e-8 3.72529e-9 -7.45058e-9
-4.65661e-10 0.0 0.0 1.30385e-8 3.72529e-9 3.72529e-9
-2.32831e-10 3.72529e-9 -1.86265e-9 1.0 0.0 0.0
-4.65661e-10 9.31323e-9 -5.58794e-9 … 1.49012e-8 1.0 -7.45058e-9
-6.98492e-10 5.58794e-9 -1.30385e-8 1.11759e-8 0.0 1.0

eigvals =
7-element Array{Float64,1}:
-3.86428e-5
-2.78764e-5
-1.69139e-5
0.000363563
0.0108282
4.69211
9201.79

cond(A) = 5.440367717196925e8
rank(A) = 7

So we see that:
(1) Ainv agrees with the result of your code for itest=2, but not for
itest = 3; and,
(2) the condition number cond(A) is very large (on the order of 10^8),
so any contamination by single-precision literals or variables may change the result...

Best,
S

[septc...@gmail.com]

But if we add another test (itest = 0) with a trivial
block-diagonal matrix (as suggested above),

selectcase ( itest )

case ( 0 )
t(:,:) = 0.0d0
t( 1, 1:2 ) = [ 0.0d0, 1.0d0 ]
t( 2, 1:2 ) = [ 1.0d0, 0.0d0 ]
do i = 3, 7
t( i, i ) = 1.0d0
enddo

we get a lot of NaN!


---------------------------------------------------
>>>>>>>>Input Matrix<<<<<<<<
0.0000 1.0000 0.0000 0.0000 0.0000 0.0000 0.0000
1.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000
0.0000 0.0000 1.0000 0.0000 0.0000 0.0000 0.0000
0.0000 0.0000 0.0000 1.0000 0.0000 0.0000 0.0000

0.0000 0.0000 0.0000 0.0000 1.0000 0.0000 0.0000

0.0000 0.0000 0.0000 0.0000 0.0000 1.0000 0.0000

0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 1.0000

---------------------------------------------------
Ainv =
NaN NaN NaN NaN NaN NaN NaN
NaN NaN NaN NaN NaN NaN NaN
NaN NaN NaN NaN NaN NaN NaN
NaN NaN NaN NaN NaN NaN NaN
NaN NaN NaN NaN NaN NaN NaN
NaN NaN NaN NaN NaN NaN NaN
NaN NaN NaN NaN NaN NaN NaN

[mirzay...@gmail.com]

Hi S,

I followed your instructions and what I got is "Exact" Inverse. It Also yield Identity on multiplication with original matrix. You solved my problem. I was struck here my all day. I appreciate you effort and work.

Best,
Mirza

[mirzay...@gmail.com]

Hi S,

Most probably, I will not come across such a condition in which I have to deal with such a matrix as you described above and get NaN in output. But as a piece of knowledge, then what should I do when facing such matrix?

Best,
Mirza

[mirzay...@gmail.com]

On Saturday, October 22, 2016 at 12:45:43 AM UTC+5, septc...@gmail.com wrote:

>>In this case, should I use a check whether first element is =0.0 then stop?

[septc...@gmail.com]

Hi,

I have no experience to write a matrix inverse routine by myself, so have little knowledge
about the underlying algorithm (sorry!). But I guess the "pivot" thing that many people
above mentioned is probably important? (such as to deal with the case of 0 in the
diagonal elements etc). I guess texts in Numerical Recipes or other numerical analysis
books will help on this...

Best,
Sept

[campbel...@gmail.com]

The following is a simple test of the matrix you provided.
It appears to work to the accuracy of the values provided, as I inverted twice and reported the difference. (hopefully no errors !!)
{Program}
! program to test matrix invert for small dense matrix
!
integer*4 :: n ! dimension of matrix
real*8, allocatable :: A_orig(:,:) ! original matrix
real*8, allocatable :: A_inv(:,:) ! inverse matrix
real*8, allocatable :: A_new(:,:) ! inverse of inverse
real*8, allocatable :: A_dif(:,:) ! difference matrix
!
integer*4 i
!
! Read in the matrix
open (unit=11, file='matinv.txt', status='OLD', action='READ')
read (11,*) n
ALLOCATE ( A_orig(n,n) )
do i = 1,n
read (11,*) A_orig(i,1:n)
end do
call print_matrix ( n, A_orig, 'Original Matrix')
!
! Invert matrix
ALLOCATE ( A_inv(n,n) )
call invert_matrix ( n, A_orig, A_inv )
call print_matrix ( n, A_inv, 'Inverse Matrix')
!
! Multiply matrix
ALLOCATE ( A_new(n,n) )
A_new = matmul ( A_orig, A_inv )
call print_matrix ( n, A_new, 'Test Matrix')
!
call test_matrix ( n, a_new, 1.0d0 )
!
! Invert again
call invert_matrix ( n, A_inv, A_new )
call print_matrix ( n, A_new, 'Initial Matrix')
!
! Subtract A_orig - A_new
A_dif = A_orig - A_new
call test_matrix ( n, a_dif, 0.0d0 )
!
end

subroutine print_matrix ( n, A, name)
integer*4 n
real*8 A(n,n)
character name*(*)
!
integer*4 i
write (*,10) name,n,n
10 format (/'Matrix ',a,':',i0,'x',i0)
11 format (i2,99(2x,f0.4))
do i = 1,n
write (*,11) i,A(i,:)
end do
end subroutine print_matrix

subroutine invert_matrix ( n, A_orig, A_inv )
integer*4 n
real*8 A_orig(n,n), A_inv(n,n)
!
real*8, allocatable :: A(:,:)
integer*4 nn, i,j
!
! Set up work array [ A_orig | Identity ]
nn = n*2
allocate ( A(n,nn) )
!
do j = 1,n
do i = 1,n
A(i,j) = A_orig(i,j)
A(i,j+n) = 0
if ( i==j ) A(i,j+n) = 1
end do
end do
!
do i = 1,n
! reduce active row
a(i,i+1:nn) = a(i,i+1:nn)/a(i,i)
a(i,i) = 1
! remove this row from all other rows
do j = 1,n
if ( j==i ) cycle
a(j,i+1:nn) = a(j,i+1:nn) - a(i,i+1:nn)*a(j,i)
a(j,i) = 0
end do
end do
!
! copy inverse to A_inv
do j = 1,n
A_inv(:,j) = a(:,j+n)
end do
!
end subroutine invert_matrix

subroutine test_matrix ( n, a, d )
integer*4 n
real*8 A(n,n), d
!
integer*4 i,j
real*8 eii, aii
!
eii = 0
do i = 1,n
do j = 1,n
aii = A(i,j)
if ( i==j ) aii = aii-d
if ( abs(aii) > abs(eii) ) eii = aii
end do
end do
write (*,12) 'max error =',eii
12 format (a,es14.3)
END subroutine test_matrix

[matinv.txt]
7

[campbel...@gmail.com]

On Sunday, October 23, 2016 at 2:31:02 PM UTC+11, campbel...@gmail.com wrote:
> The following is a simple test of the matrix you provided.

I should not have edited the post without testing, as I omitted the ALLOCATE line.

! Subtract A_orig - A_new

ALLOCATE ( A_dif(n,n) )

A_dif = A_orig - A_new

These are not large matrices, so extra copies of the A matrix is not an issue.
As the final round-off error of my run is 2.2e-5, pivoting is not a requirement.
I could not reproduce the real*4 error case.

[campbel...@gmail.com]

You can introduce a form of partial pivoting, which solves for the identified (itest = 0) style of matrix with:
integer p
real*8 t
...

do i = 1,n
!

! find best active row
p = i
do j = i+1,n
if ( abs(a(j,i)) > abs(a(p,i)) ) p = j
end do
! swap rows
if ( p /= i ) then
write (*,*) 'Pivot : Swap rows',i,p
do j = 1,nn
t = a(i,j)
a(i,j) = a(p,j)
a(p,j) = t
end do
end if
!
! reduce active row

This works for:
3
0.,1.,0.
1.,0.,0.
.1,.1,1.

producing:
Matrix Original Matrix:3x3
1 0.0000 1.0000 0.0000
2 1.0000 0.0000 0.0000
3 0.1000 0.1000 1.0000
Pivot : Swap rows 1 2

Matrix Inverse Matrix:3x3
1 0.0000 1.0000 0.0000
2 1.0000 0.0000 0.0000
3 -0.1000 -0.1000 1.0000

Matrix Test Matrix:3x3
1 1.0000 0.0000 0.0000
2 0.0000 1.0000 0.0000
3 0.0000 0.0000 1.0000
max error = 0.000E+00
Pivot : Swap rows 1 2

Matrix Initial Matrix:3x3
1 0.0000 1.0000 0.0000
2 1.0000 0.0000 0.0000
3 0.1000 0.1000 1.0000
max error = 0.000E+00

Also, using pivoting the original 7x7 matrix produces:
Matrix Original Matrix:7x7
1 744.5047 909.5175 983.0251 1024.6988 1051.5530 1070.3036 1084.1393
2 909.5175 1112.7881 1203.5504 1255.0658 1288.2851 1311.4910 1328.6199
3 983.0251 1203.5504 1302.1218 1358.0994 1394.2078 1419.4376 1438.0635
4 1024.6988 1255.0658 1358.0994 1416.6292 1454.3913 1480.7802 1500.2639
5 1051.5530 1288.2851 1394.2078 1454.3913 1493.2255 1520.3664 1540.4069
6 1070.3036 1311.4910 1419.4376 1480.7802 1520.3664 1548.0349 1568.4664
7 1084.1393 1328.6199 1438.0635 1500.2639 1540.4069 1568.4664 1589.1877
Pivot : Swap rows 1 7
Pivot : Swap rows 2 7
Pivot : Swap rows 5 7
Pivot : Swap rows 6 7

Matrix Inverse Matrix:7x7
1 -76.5341 1348.1077 -2930.3976 471.9861 1596.9489 914.7758 -1319.4799
2 1348.1077 -9685.4287 16217.1709 -444.2072 -8633.6875 -6213.9709 7423.6791
3 -2930.3976 16217.1709 -25394.9467 5472.2177 8223.5511 4896.2298 -6548.5709
4 471.9861 -444.2072 5472.2177 -20792.6765 14657.9506 11719.0734 -11047.5256
5 1596.9489 -8633.6875 8223.5511 14657.9506 -21308.0560 -810.0951 6302.8896
6 914.7758 -6213.9709 4896.2298 11719.0734 -810.0951 -28224.5400 17718.8557
7 -1319.4799 7423.6791 -6548.5709 -11047.5256 6302.8896 17718.8557 -12548.3771

Matrix Test Matrix:7x7
1 1.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000
2 0.0000 1.0000 0.0000 0.0000 0.0000 0.0000 0.0000
3 0.0000 0.0000 1.0000 0.0000 0.0000 0.0000 0.0000
4 0.0000 0.0000 0.0000 1.0000 0.0000 0.0000 0.0000
5 0.0000 0.0000 0.0000 0.0000 1.0000 0.0000 0.0000
6 0.0000 0.0000 0.0000 0.0000 0.0000 1.0000 0.0000
7 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 1.0000
max error = -1.620E-08
Pivot : Swap rows 1 3
Pivot : Swap rows 3 4
Pivot : Swap rows 4 6

Matrix Initial Matrix:7x7
1 744.5047 909.5175 983.0251 1024.6988 1051.5530 1070.3036 1084.1393
2 909.5175 1112.7881 1203.5504 1255.0658 1288.2851 1311.4910 1328.6199
3 983.0251 1203.5504 1302.1218 1358.0994 1394.2078 1419.4376 1438.0635
4 1024.6988 1255.0658 1358.0994 1416.6292 1454.3913 1480.7802 1500.2639
5 1051.5530 1288.2851 1394.2078 1454.3913 1493.2255 1520.3664 1540.4069
6 1070.3036 1311.4910 1419.4376 1480.7802 1520.3664 1548.0349 1568.4664
7 1084.1393 1328.6199 1438.0635 1500.2639 1540.4069 1568.4664 1589.1877
max error = 1.910E-06

This looks reasonably stable, although full pivoting is still a possibility, which I have never required.
Calculating the inverse for small problems can be convenient and should not be an issue.