Real Symmetric Matrix Inversion : Fortran

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

"rjmagyar" <rjmagyar@gmail.com> wrote in message
news:d8fc7e65-caa3-4307-8ac4-8f45f05244fa@e60g2000hsh.googlegroups.com...
> 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

For n < 50
below is PROVEN fastest and most accurate..


! --------------------------------------------------------------------
SUBROUTINE Gauss (a,n) ! Invert matrix by Gauss method
! --------------------------------------------------------------------
IMPLICIT NONE
INTEGER :: n
REAL(8) :: a(n,n), b(n,n), c, d, temp(n)
INTEGER :: j, k, m, imax(1), ipvt(n)

b = a
ipvt = [ 1:n ]

DO k = 1,n
imax = MAXLOC(ABS(b(k:n,k)))
m = k-1+imax(1)
IF (m /= k) THEN
ipvt( [m,k] ) = ipvt( [k,m] )
b( [m,k], = b( [k,m],
END IF
d = 1/b(k,k)
temp = b(:,k)
DO j = 1, n
c = b(k,j)*d
b(:,j) = b(:,j)-temp*c
b(k,j) = c
END DO
b(:,k) = temp*(-d)
b(k,k) = d
END DO
a(:,ipvt) = b
END SUBROUTINE Gauss

On 2008-03-20 16:04:19 -0300, rjmagyar <rjmagyar@gmail.com> said:

> 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

The usual advice in order would be:

1. You are more likely to get advice on numerical methods from a
numerical methods newsgroup. Try sci.math.num-****ysis. c.l.f
has many folks who are knowlegable in NA but that is more good
luck on your part than good planning.

2. Do not invert unless you have extremely good reasons to so.
Solving equations by inversion is both slower and less accurrate
than it need be.

3. Try the netlib collection.

4. Consult your local numerical ****ysis expert. That rarely means
the guy you met standing in line for lunch who you overheard talking
about programming.

If your matrix is a graded as you suggest then you are in bad need of
an expert and may need much more than merely conversion of some found
code into double precision. Is your matrix only symmetric or is it
positive definite symmetric? It matters a lot and the fact that you did
not mention this issue reinforces the notion that you need real advice.

rjmagyar wrote:

> 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.

I don't know the specific routine, but there are iterative methods
which do the inversion, then correct it iteratively. My guess is
that with 10 orders of magnitude even that won't be enough with
only double precision. You didn't say how big it is, though.

I would suggest a multiple precision (software) package, though I
don't know that I have ever seen matrix inversion written for one.

As others have said, the numerical ****ysis newsgroup will
give better answers.

-- glen

"rjmagyar" <rjmagyar@gmail.com> schrieb im Newsbeitrag news:d8fc7e65-caa3-4307-8ac4-8f45f05244fa@e60g2000hsh.googlegroups.com...
> 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

Rudy,

Have a look at LAPACK.
http://www.netlib.org/lapack/

The have several different inversion methods. Some of those come with an automatic "rescaling" of the matrix to cope with
difference in magnitude (but only of the diagonal elements I believe).

However, please take the advice of the previous posters seriously! Only invert if you really have to. Solving a matrix is (much)
faster and more precise.

Regarding LAPACK (which also needs BLAS). Most compilers have their "own" version of this so there is normally no need to install
it yourself (e.g. on Sun it is contained in the "perflib", the intel compiler has it in its "mkl" library). However, installing it
yourself is very easy and almost completely automatic.

Good luck!
Tim
http://gnss.servolux.nl

David Frank wrote:
> "rjmagyar" <rjmagyar@gmail.com> wrote in message
> news:d8fc7e65-caa3-4307-8ac4-8f45f05244fa@e60g2000hsh.googlegroups.com...
>> 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
>
> For n < 50
> below is PROVEN fastest and most accurate..

For a 1x1 matrix, your claim above is demonstrably false. I imagine there are
other holes in your claim.

>
>
> ! --------------------------------------------------------------------
> SUBROUTINE Gauss (a,n) ! Invert matrix by Gauss method
> ! --------------------------------------------------------------------
> IMPLICIT NONE
> INTEGER :: n
> REAL(8) :: a(n,n), b(n,n), c, d, temp(n)
> INTEGER :: j, k, m, imax(1), ipvt(n)
>
> b = a
> ipvt = [ 1:n ]
>
> DO k = 1,n
> imax = MAXLOC(ABS(b(k:n,k)))
> m = k-1+imax(1)
> IF (m /= k) THEN
> ipvt( [m,k] ) = ipvt( [k,m] )
> b( [m,k], = b( [k,m],
> END IF
> d = 1/b(k,k)
> temp = b(:,k)
> DO j = 1, n
> c = b(k,j)*d
> b(:,j) = b(:,j)-temp*c
> b(k,j) = c
> END DO
> b(:,k) = temp*(-d)
> b(k,k) = d
> END DO
> a(:,ipvt) = b
> END SUBROUTINE Gauss
>
>
>

Rich Townsend <rhdt@barVOIDtol.udel.edu> wrote:

> David Frank wrote:

> > For n < 50
> > below is PROVEN fastest and most accurate..
>
> For a 1x1 matrix, your claim above is demonstrably false. I imagine there are
> other holes in your claim.

No doubt. I have had DF killfilled for a long time. I don't usually even
see his postings unless someone quotes them. He has a long history of
doing severely biased tests, just ignoring any data that he doesn't
like, severely misinterpreting what data he gets (which is why lots of
people stopped sending him data, which, of course, he then misinterprets
as "proof" that they have nothing to contradict him), and refusing to
study the work of actual experts who have spent their careers on the
issues in question.

I don't consider it worth even looking to see what DF's misstatements
are any more; he has used up more than his quota of my time.

--
Richard Maine | Good judgement comes from experience;
email: last name at domain . net | experience comes from bad judgement.
domain: summertriangle | -- Mark Twain

Actually I kind of have the same problem... My matrix has about 500 elements, all real dp, and span a range from 10e5 to 10e-3, with a lot of zero elements (its for a hamiltonian). Anything to recommend?