Compute LU factorization; then solve system

[diegosp]

Hi,

I have a loop and, at each step, I have to solve a linear system, Ax=b. However, at each step of the loop, only b changes, and A always remains the same.
I am using LAPACK, and I've found some subroutines like "sgesv" which solve linear systems by performing first the LU factorization of the coefficients Matrix. However, I don't want to use this subroutine at each step of the loop. Instead, since A is always the same, I would like to compute just the LU factorization in the first step, store L and U, and then at the next steps solve the systems with the factorization already computed. I am guessing this would be far more efficient than solving the system at each step (I'm dealing with very large matrices and a long loop)

I have find searching for this info for some time, and I have not been able to find the answer. For example, in MatLab, performing just the LU factorization is as simple as writing [L,U]=lu(A). Is there not something like this in Fortran / LAPACK? I've written my own LU and solver subroutines, which is not hard, but I guess LAPACK's are very optimized, and I guess what I want do has to be something doable in the context of LAPACK.

Thank you

 

[mikrom]

Hi diegosp,

We can look into the source of SGESV - for example here

http://www.netlib.org/lapack/explore-3.1.1-html:

modified code snippet from sgesv.f

...
[COLOR=#0000ff]*     Step 1: Compute the LU factorization of A.[/color]
      [COLOR=#008b8b]CALL[/color] SGETRF(N, N, A, LDA, IPIV, INFO)
      [COLOR=#a52a2a][b]IF[/b][/color]( INFO[COLOR=#a52a2a][b].EQ.[/b][/color][COLOR=#ff00ff]0[/color] ) [COLOR=#a52a2a][b]THEN[/b][/color]
[COLOR=#0000ff]*       Step 2: Solve the system A*X = B, overwriting B with X.[/color]
        [COLOR=#008b8b]CALL[/color] SGETRS([COLOR=#ff00ff]'No transpose'[/color], N, NRHS, A, LDA, IPIV, B, LDB, INFO)
      [COLOR=#a52a2a][b]END IF[/b][/color]
...

As we can see from the code snippet above, the routine SGESV consists of 2 main steps:
In 1.step it calls SGETRF, which computes LU factorization and then
in 2.step it calls SGETRS which computes the solution using the LU factorization computed by SGETRF

IMO, you could do the same in your program:
Instead of calling SGESV in the loop, you could call only once before the loop SGETRF to obtain the LU decomposition of your matrix and then in the loop call SGETRS with your LU decomposited matrix.