wanted fast routines

Can anyone provide me with with a superfast routine to find inverses of
symmetric and positive definite matrices. I am currently using the
subroutine provided by numerical recipes website.

Also, is the intrinsic Matrix Multiplication function (MATMUL) optimal?

In article <1118013746.166644.228610@g47g2000cwa.googlegroups.com>,
"student" <adarsh@stat.tamu.edu> writes:
>
> Can anyone provide me with with a superfast routine to find inverses of
> symmetric and positive definite matrices. I am currently using the
> subroutine provided by numerical recipes website.
>
> Also, is the intrinsic Matrix Multiplication function (MATMUL) optimal?

http://www.netlib.org/

Look for ATLAS and LAPACK.

--
Steve
http://troutmask.apl.washington.edu/~kargl/

student wrote:
> Can anyone provide me with with a superfast routine to find inverses of
> symmetric and positive definite matrices. I am currently using the
> subroutine provided by numerical recipes website.

A rule of thumb that was passed to me with regards to numerical matrix
programming is that if your algorithm needs to actually find the inverse
of a matrix, you're programming it wrong -- and avoiding the matrix
inversion will make it much faster.  This isn't always true, of course,
but it was claimed to be true nearly always.

- Brooks


--
The "bmoses-nospam" address is valid; no unmunging needed.

Brooks Moses wrote:
> student wrote:
> 
>
>
> A rule of thumb that was passed to me with regards to numerical matrix
> programming is that if your algorithm needs to actually find the inverse
> of a matrix, you're programming it wrong -- and avoiding the matrix
> inversion will make it much faster.  This isn't always true, of course,
> but it was claimed to be true nearly always.
>
> - Brooks
>
>

To be a little more specific for the OP, if you have the linear equation

A x = b,

where A is a matrix, and x and b are vectors or matrices, then it is
very inefficient to solve for x by (i) inverting A and then (ii)
multiplying b by A^-1. It is much better to use an algorithm that is
specifically targeted at solving linear equations, such as a
decomposition followed by a backsubstitution. This and related issues
are disussed in Numerical Recipes -- although when it comes to actual
code, you may be well served by avoiding NR and going for an appropriate
LAPACK routine.

cheers,

Rich

"student" <adarsh@stat.tamu.edu> wrote in message
news:1118013746.166644.228610@g47g2000cwa.googlegroups.com...
>
> Can anyone provide me with with a superfast routine to find inverses of
> symmetric and positive definite matrices. I am currently using the
> subroutine provided by numerical recipes website.
Evidently, most applications which require a conceptual inverse will be
better implemented by some solver which doesn't actually form the inverse.
Numerical Recipes has an introductory discussion for the case of
eigensystems, and I don't see them proposing any one method as fitting all
needs.
>
> Also, is the intrinsic Matrix Multiplication function (MATMUL) optimal?
Possibly, for a definition of optimality chosen by the individual who
implemented it for your compiler.  If you use gfortran, you have the
opportunity to examine the source, form your own opinion, and try
variations.  Remember that MATMUL is likely to need to perform well for
matrices with small sizes, while the programmer has the option to use
popular BLAS and other open source calls when those are known to be
suitable.  A few compilers will choose a call to the vendor's optimized BLAS
when the size of the matrix is known to be large enough, but expand in line
when that is likely to be preferable.

student wrote:

> Also, is the intrinsic Matrix Multiplication function (MATMUL) optimal?

Matrix multiplication is sensitive to the cache of the processor it
runs on, and so is often not optimal.  I believe there are some that
try to dynamically optimize for the cache, though.  You might look
for those.

-- glen

"student" <adarsh@stat.tamu.edu> wrote in message
news:1118013746.166644.228610@g47g2000cwa.googlegroups.com...
>
> Can anyone provide me with with a superfast routine to find inverses of
> symmetric and positive definite matrices. I am currently using the
> subroutine provided by numerical recipes website.
>
> Also, is the intrinsic Matrix Multiplication function (MATMUL) optimal?
>

For n < 100 below Fortran source for matrix inverse is BRIEFER/FASTER with
less ERROR
compared to any other Fortran source anyone will post.

! ---------------------------
SUBROUTINE Gauss (a,n)       ! Invert matrix by Gauss method
IMPLICIT NONE
INTEGER :: n
REAL(8) :: a(n,n), b(n,n), temp(n), c, d
INTEGER :: ipvt(n), j, k, m
! - - - - - - - - - - - - - -
b = a
ipvt = [ 1:n ]
DO k = 1,n
m = k-1+MAXLOC(ABS(b(k:n,k)),dim=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

David Frank wrote:
...
> For n < 100 below Fortran source ...

...
> ipvt = [ 1:n ]
...
>       ipvt( [m,k] ) = ipvt( [k,m] )
>       b( [m,k],:) = b( [k,m],:)
...

Uhm. Fortran source ?... with [..]  ?? What does it mean ?
If I am not wrong, standard Fortran 95 doesn't allow [] or {}.


Giorgio

When I had to multiply some triangular matrices with a vector, I made
some test, and the best result I got, was mixing a do loop for one
index, and sum for the other.

"Giorgio Pastore" <pastgio@univ.trieste.it> wrote in message
news:d81mdf$lqb$1@newsfeed.cineca.it...
> David Frank wrote:
> ... 
>
> ... 
> ... 
> ...
>
> Uhm. Fortran source ?... with [..]  ?? What does it mean ?
> If I am not wrong, standard Fortran 95 doesn't allow [] or {}.
>
>
> Giorgio

The is acceptable to CVF (a F95+ compiler)

I think replacing (/.../)  with [...]  is supported by a few other F95+
compilers, it definitely must
be accepted by a CURRENT STANDARD compiler..

ipvt = [ 1:n]     creates a 1....n  array and assigns it to ipvt
you can replace it with
ipvt = (/ (i,i=1,n) /)

again this is current standard fortran supported by a few F95+ compilers..