Matrix Inversion Suggestions for Fortran

rjmagyar

unread,

Mar 21, 2008, 12:10:32 AM3/21/08

to

Hi, I was wondering if any one can recommend a robust and fast matrix
inversion routine. I have a matrix with vastly different (sometimes
over 10 orders of magnitude) elements that needs to be inverted. I
have modified the matinv routine by C. Reeve to allow double precision
calculations, but this doesn't seems to be enough. The inverse I get
is only approximately correct. Does any one know of either a general
routine that can handle this, or does any one know a scheme around
this?

Best, Rudy

user923005

unread,

Mar 21, 2008, 12:22:16 AM3/21/08

to

The MIRACL library can invert a large Hilbert matrix, because it uses
multiple precision rational arithmetic.
http://www.shamus.ie/

Are you sure that the matrices are not actually singular?

Evgenii Rudnyi

unread,

Mar 21, 2008, 4:55:23 AM3/21/08

to

For what purpose do you need the matrix inverse? Could it be that you
need to solve just a system of linear equations?

As for accuracy, you have to check the matrix condition number. For
example, compute its singular values and then the ration sigma_max/
sigma_min. If this ratio is high, then you have to rethink your
original problem.

none

unread,

Mar 21, 2008, 11:09:26 AM3/21/08

to

'Numercail Recipes' will make a good start, offering a variety of
alternatives if you know the particular problem that you have. If you
really don't, then worth looking up are iterative refinement - which
DEMANDS that if you solve the inverse using DOUBLE (DOUBLE PRECISION)
then you MUST have acess to LONG FLOAT (REAL*10 or better) for the
iterative refinement. A useful alternate, which will help give you a 'good
enough for your working precision' alternative is SVD.

aruzinsky

unread,

Mar 21, 2008, 1:44:18 PM3/21/08

to

In the rare cases I needed an inverse, I would use LUP decomposition
and solve

LUP A^-1 = I

rmau...@googlemail.com

unread,

Mar 21, 2008, 4:27:08 PM3/21/08

to

Numerical Recipes also has a chapter on iterative improvement.

http://www.nrbook.com/a/PDF4.gif

rmau...@googlemail.com

unread,

Mar 21, 2008, 4:30:20 PM3/21/08

to

The Fortran version of this chapter (from an obsolete edition) is
here:

http://www.nrbook.com/a/bookfpdf/f2-5.pdf

aruzinsky

unread,

Mar 22, 2008, 12:37:58 PM3/22/08

to

I can't help but think that if iterative improvement worked with LUP
decomposition, iterative improvement would be a lot more popular.

In my experience, Numerical Recipies, is not relatively good. Many
better fundamental algorithms can be found in Collected Algorithms
From ACM Volumes I - III. The algorithms date back to the 1960s and
many were written in algol.

Ken

unread,

Mar 22, 2008, 12:49:04 PM3/22/08

to

On Mar 20, 9:10 pm, rjmagyar <rjmag...@gmail.com> wrote:

Since you mentioned Fortran, LAPACK has robust routines that includes
matrix inversion. It also has routines for scaling to try to reduce
the condition number.

http://www.netlib.org/lapack/double/

Also, netlib has multi-precision packages.

Ken

unread,

Mar 22, 2008, 1:04:05 PM3/22/08

to

NAG has some excellent documentation on the use of LAPACK routines:

http://www.hep.phy.cam.ac.uk/NAGdoc/fl/html/F07_fl19.html

Victor Eijkhout

unread,

Mar 23, 2008, 4:21:43 PM3/23/08

to

rjmagyar <rjma...@gmail.com> wrote:

> Hi, I was wondering if any one can recommend a robust and fast matrix
> inversion routine.

They don't exist. Inversion is intrinsically less stable than LU
factorizaitons. If you only need to solve a linear system, do not invert
it.

Victor.
--
Victor Eijkhout -- eijkhout at tacc utexas edu

aruzinsky

unread,

Mar 24, 2008, 1:13:27 PM3/24/08

to

On Mar 23, 2:21 pm, s...@sig.for.address (Victor Eijkhout) wrote:

How can that be when inversion is accomplished by LU solution on the
columns of I?

Evgenii Rudnyi

unread,

Mar 24, 2008, 2:42:34 PM3/24/08

to

This is a typical way to find an inverse. From the definition

A A^-1 = I

you have a system of linear equations with A and with many right hand
sides formed by I. As an answer to this, one obtains A^-1. Clearly the
way to do it, is first to make the LU decomposition of A first and
then use back substitutions.

If you like something unusual, look at this paper

Tibor Mazúch, Jan Kozánek
New recurrent algorithm for a matrix inversion
Journal of Computational and Applied Mathematics
Volume 136 , Issue 1-2 (November 2001)
Pages: 219 - 226
Year of Publication: 2001

aruzinsky

unread,

Mar 25, 2008, 11:49:01 AM3/25/08

to

On Mar 24, 12:42 pm, Evgenii Rudnyi <use...@rudnyi.ru> wrote:
> On 24 Mrz., 18:13, aruzinsky <aruzin...@general-cathexis.com> wrote:
>
>
>
>
>
> > On Mar 23, 2:21 pm, s...@sig.for.address (Victor Eijkhout) wrote:
>
> > > rjmagyar <rjmag...@gmail.com> wrote:
> > > > Hi, I was wondering if any one can recommend a robust and fast matrix
> > > > inversion routine.
>
> > > They don't exist. Inversion is intrinsically less stable than LU
> > > factorizaitons. If you only need to solve a linear system, do not invert
> > > it.
>
> > > Victor.
> > > --
> > > Victor Eijkhout -- eijkhout at tacc utexas edu
>
> > How can that be when inversion is accomplished by LU solution on the
> > columns of I?
>
> This is a typical way to find an inverse. From the definition
>
> A A^-1 = I
>
> you have a system of linear equations with A and with many right hand
> sides formed by I. As an answer to this, one obtains A^-1. Clearly the
> way to do it, is first to make the LU decomposition of A first and
> then use back substitutions.
>

I already said that so what is your point?

Evgenii Rudnyi

unread,

Mar 26, 2008, 2:57:03 AM3/26/08

to

> I already said that so what is your point?-

I am very sorry, I have misinterpreted your question.

Evgenii