Computing inverse of a matrix using LAPACK,BLAS in Fortran90

Question

I am new to using LAPACK/BLAS. I want to compute the solution of an equation:

AU=F

I want to know what is the logical error in this part of the code. Where I use the solvers to input the matrix A which is of size ((xdiv-1)(ydiv-1),(xdiv-1)(ydiv-1)). Then, subsequently solve for the Equation. U=inverse(A) *f.

Where U and f are of the same size.(u((xdiv-1)(ydiv-1),1),f((xdiv-1)(ydiv-1),1)). I am getting a segmentation fault error while performing the matrix inverse.

Here is my code:

program main
double precision, allocatable :: A(:,:)
double precision, allocatable :: u(:,:), f(:,:)
double precision mesh(2), dx, dy
integer  xdiv, ydiv
 xdiv=55
 ydiv=55
 mesh(1)=.001
 mesh(2)=.001
 dx=mesh(1)
 dy=mesh(2)

allocate (A((xdiv-1)*(ydiv-1),(xdiv-1)*(ydiv-1))) 
allocate (Ainv((xdiv-1)*(ydiv-1),(xdiv-1)*(ydiv-1)))
allocate (u((xdiv-1)*(ydiv-1),1),f((xdiv-1)*(ydiv-1),1))

         do i =1,(xdiv-1)*(ydiv-1)
         A(i,i)=-2.d0*(1.d0/(dx**2)+1.d0/(dy**2))
         enddo
         do i=1,(xdiv-2)
         do j=1,(ydiv-1)
         A(i+(j-1)*(xdiv-1),i+(j-1)*(xdiv-1)+1)=1.d0/(dx**2)
         A(i+(j-1)*(xdiv-1)+1,i+(j-1)*(xdiv-1))=1.d0/(dx**2)
         enddo
         enddo
         do i=1,(xdiv-1)
         do j=1,(ydiv-2)
         A(i+(j-1)*(xdiv-1),i+(j)*(xdiv-1))=1.d0/(dy**2)
         A(i+(j)*(xdiv-1),i+(j-1)*(xdiv))=1.d0/(dy**2)
         enddo
         enddo

        do i=1,(xdiv-1)
         do j=1,(ydiv-1)
       xcoord = (i-1)*mesh(1)
       ycoord = (j-1)*mesh(2)
       xd=sin(2.0*3.14*xcoord)
       yd=sin(2.0*3.14*ycoord)
       f((i-1)*(xdiv-1)+j,1)= 8.0*3.1415*3.1415*xd*yd
         enddo
        enddo
    do i=1,(xdiv-1)*(ydiv-1)
     do j=1,(xdiv-1)*(ydiv-1)
      Ainv(i,j)=A(i,j)
     end do
    end do

 call DGETRF((xdiv-1)*(ydiv-1), (xdiv-1)*(ydiv-1), Ainv,  & 
  (xdiv-1)*(ydiv-1), ipiv, info)

 call DGETRI((xdiv-1)*(ydiv-1), Ainv, (xdiv-1)*(ydiv-1), ipiv, &
  work,(xdiv-1)*(ydiv-1), info)

 call DGEMV(N, (xdiv-1)*(ydiv-1), (xdiv-1)*(ydiv-1), 1.d0, Ainv, &
  (xdiv-1)*(ydiv-1), f, 1.d0, 0.d0, u, 1.d0)

  do i=1,(xdiv-1)*(ydiv-1)
   write(*,*) "u", u(i,1)
  enddo
 end program main

Basically, I compute the LU Decomposition, then invert it and then multiply. Please help me find the error and do suggest a better way to do this calculation if any.

NB: Some parts of variable defining/assigning values may seem redundant or inefficient that is because this is a part of a bigger code. I just extracted portions of it, to focus on the matrix inversion issue.

Please use tag fortran for all Fortran questions. Your code is very incomplete We must see how are all those variables declared and we should be able to compile your code and test it. Please read minimal reproducible example and How to Ask. — Vladimir F Героям слава
– Vladimir F Героям слава, Commented May 2, 2017 at 6:20
@VladimirF Thanks for your comments and suggestions. I have added the fortran tag and put up an entire code for anybody to compile and test. — Schneider
– Schneider, Commented May 2, 2017 at 7:20
Is there a good reason you are doig the factorisation and invert yourself rather than using dgesv for the whole process in one go? — Ian Bush
– Ian Bush, Commented May 2, 2017 at 7:22
Also if you want help I strongly suggest you use Implicit None in the code above - I can immediately see one mistake that including that might help you solve. — Ian Bush
– Ian Bush, Commented May 2, 2017 at 7:29
@IanBush No there isn't any specific reason for doing the decomposition by myself. Thanks for suggesting implicit none, could you mention what is the error? — Schneider
– Schneider, Commented May 2, 2017 at 23:53

b-fg · Accepted Answer · 2017-05-03 17:40:07Z

6

For the inversion of a general squared matrix have a look at this function (which uses LAPACK)

  function inv(A) result(Ainv)
    implicit none
    real,intent(in) :: A(:,:)
    real            :: Ainv(size(A,1),size(A,2))
    real            :: work(size(A,1))            ! work array for LAPACK
    integer         :: n,info,ipiv(size(A,1))     ! pivot indices

    ! Store A in Ainv to prevent it from being overwritten by LAPACK
    Ainv = A
    n = size(A,1)
    ! SGETRF computes an LU factorization of a general M-by-N matrix A
    ! using partial pivoting with row interchanges.
    call SGETRF(n,n,Ainv,n,ipiv,info)
    if (info.ne.0) stop 'Matrix is numerically singular!'
    ! SGETRI computes the inverse of a matrix using the LU factorization
    ! computed by SGETRF.
    call SGETRI(n,Ainv,n,ipiv,work,n,info)
    if (info.ne.0) stop 'Matrix inversion failed!'
end function inv

Also, you can find other useful Fortran subroutines and functions of matrix operations here.

Basically, this follows the same procedure as you did except that, as mentioned in some comments, the arguments you are passing to the LAPACK subroutines are not correct. Have a look at those and fix them or just use this function.

EDIT: Note that this is for single precision data types. You can easily adapt it for double precision.

edited May 3, 2017 at 17:40

answered May 3, 2017 at 17:34

b-fg

4,1972 gold badges31 silver badges45 bronze badges

Sign up to request clarification or add additional context in comments.

2 Comments

TheWhitestOfFangs Over a year ago

Given that a matrix A is diagonally-dominant, is there a way to improve this?

b-fg Over a year ago

I am not sure about that, sorry

Collectives™ on Stack Overflow

Computing inverse of a matrix using LAPACK,BLAS in Fortran90

1 Answer 1

2 Comments

Your Answer

Hot Network Questions

Collectives™ on Stack Overflow

1 Answer 1

2 Comments

Your Answer

Sign up or log in

Post as a guest

Related