I am using LAPACK to solve a linear sytem of equations A(n,n) * X(n,1) = B(n,1). A, X and B are real*8 (double precision).
First of all, I create A and B, initialize X to 0, and I then make a copy of A and B to A0 and B0, because A and B will be overwritten.
First of all, I use dgetrf to calculate LU factorization:
call dgetrf ( n, n, A, n, iPiv, info )
with integer*4 array iPiv(n). Matrix A gets overwritten by L and U, with unit diagonal elements not stored. This is clear.
I then continue with dgetrs to solve my system: