Hello,
I have to invert a large matrix ( N > 10000 x 10000) this matrix is
square and complex. I need to use many times its inverse.
Do you really need the inverse or do you just want to solve several times
a linear system of equations but with the same coefficient matrix?
[1] If you really want the inverse: First call LAPACK ZGETRF, this will
perform an in-place LU factorization of your coefficient matrix. Then call
ZGETRI, this will use the inplace LU-factorization of your initial matrix
and returns you the inverse of the matrix A. Same thing in place as well.
So basically the code is two LAPACK calls. (ZGETRF+ZGETRI). When you want
to use the inverse, for example to solve a linear system of equations, you
just need to multiply by the matrix inverse. This can be done with ZGEMV
(one right-hand sides) or ZGEMM (multiple right-hand sides). This approach
costs: 2/3*N^3 for ZGETRF, 4/3*N^3 for ZGETRI and 2*N^2 per right-hand
sides (ZGEMV or ZGEMM).
[2] However if all you want to do is to solve one linear system of
equations with several right-hand sides, this is not the recommended way
to go. A more stable and less expensive approach is the following. First
call ZGETRF. This performs the LU factorization of your coefficient matrix
in-place and cost 2/3*N^3 FLOPS. Keep the N-by-N array computed by LAPACK
as long as you want to solve a linear system with the initial matrix. Each
time you want to solve for a given right-hand side, you need to call
ZGETRS, ZGETRF uses the LU factorization output form SGETRF. (ZGETRS can
handle one right-hand side or multiple, if multiple, the routine is then
more efficient than several call to this routine.) The cost is: 2/3*N^3
for ZGETRF and 2*N^2 per right-hand sides. Same thing two lines calls to
LAPACK. As long as you do not touch the output from ZGETRF, you are free
to reuse ZGETRS, i.e. you are free to solve another set of right-hand
sides.
So [2] is less expensive and moreover is more stable. (Skip the stability
discussion.)
If you want to see a call to ZGETRF followed to call to ZGETRS look at
the LAPACK routine ZGESV. That's basically all that routine does.
With respect to CLAPACK. Sure you can use CLAPACK. Say that if you have
LAPACK already installed on your system then you better off using it.
CLAPACK is done for arch/OS that do not have a (free/cheap/at-all) Fortran
compiler.
For this kind of sizes (N=10,000). Do not forget to use a good BLAS behind
LAPACK, if possible multithreaded.
-j
On Sun, 10 Jun 2007, Michel Yvert wrote:
Dear CLAPACK helpers,
This mail is to ask for an advice: I have to invert a large matrix ( N >
10000 x 10000) this matrix is square and complex. I need to use many
times its inverse. I asked Roldan Pozo who suggested to use CLAPACK (I
work in C++ and would like to incorporate this inversion in a C++ job).
Could you indicate me which is the best approach (the routine to use)
for a best efficiency regarding the precision and the calculation time ?
I would appreciate very much a simple example working on a similar task.
Many thanks in advance
Michel Yvert
--
=======================================================
Michel Yvert
Laboratoire d'Annecy-le-Vieux de Physique des Particules
B.P. 110
74941 Annecy-le-Vieux Cedex
France
Tel : +33 4 50 09 16 75 / 04 50 09 16 75
Fax : +33 4 50 27 94 95 / 04 50 27 94 95
http://wwwlapp.in2p3.fr
_______________________________________________
Lapack mailing list
Lapack@Domain.Removed
http://lists.cs.utk.edu/listinfo/lapack
|