INTELLIGENT WORK FORUMS
FOR COMPUTER PROFESSIONALS

Are you a
Computer / IT professional?
Join Tek-Tips Forums!
• Talk With Other Members
• Be Notified Of Responses
• Keyword Search
Favorite Forums
• Automated Signatures
• Best Of All, It's Free!

*Tek-Tips's functionality depends on members receiving e-mail. By joining you are opting in to receive e-mail.

#### Posting Guidelines

Promoting, selling, recruiting, coursework and thesis posting is forbidden.

# Compute LU factorization; then solve system

 Forum Search FAQs Links MVPs

## Compute LU factorization; then solve system

(OP)
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

### RE: Compute LU factorization; then solve system

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

#### CODE

...
*     Step 1: Compute the LU factorization of A.
CALL SGETRF(N, N, A, LDA, IPIV, INFO)
IF( INFO.EQ.0 ) THEN
*       Step 2: Solve the system A*X = B, overwriting B with X.
CALL SGETRS('No transpose', N, NRHS, A, LDA, IPIV, B, LDB, INFO)
END IF
... 

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.

#### Red Flag This Post

Please let us know here why this post is inappropriate. Reasons such as off-topic, duplicates, flames, illegal, vulgar, or students posting their homework.

#### Red Flag Submitted

Thank you for helping keep Tek-Tips Forums free from inappropriate posts.
The Tek-Tips staff will check this out and take appropriate action.

#### Posting in the Tek-Tips forums is a member-only feature.

Click Here to join Tek-Tips and talk with other members!

Close Box

# Join Tek-Tips® Today!

Join your peers on the Internet's largest technical computer professional community.
It's easy to join and it's free.

Here's Why Members Love Tek-Tips Forums:

• Talk To Other Members
• Notification Of Responses To Questions
• Favorite Forums One Click Access
• Keyword Search Of All Posts, And More...

Register now while it's still free!