!! File  : matinv.f90
!!
!! Creation Date: 31.10.2015; Time: 10:18
!! Time-stamp: <2015-10-31 11:14:46 (andrzejbaran)> 
!!
!! Andrzej Baran
!! IF UMCS, ul. Radziszewskiego 10, 20-031 Lublin, PL
!! baran@kft.umcs.lublin.pl
!!
!! See: Num. Recipes., p40.
!
subroutine matinv(A,N,NP,flag)
  USE LU
  IMPLICIT NONE
  INTEGER, PARAMETER :: NMAX = 500
  INTEGER :: NP, N, D, FLAG, i, j
  REAL, DIMENSION(NP,NP) :: A, Y
  INTEGER, DIMENSION(NMAX) :: INDX
  !
  ! This subroutine calculates inverse of a matrix (see Num. Rec. p40)     
  ! A is the matrix to be inverted and is destroyed in the routine.
  ! on output matrix A contains its inverse A^(-1).
  ! Program stops if A is singular matrix.
  !
  ! Method: lu decomposition
  ! Routines used: ludcmp, lubksb (From module lu.f90)
  !
  IF (N.GT.NP.OR.N.GT.NMAX)  STOP 'matinv: n is to large. Increase np.'
  !
  Y = 0.
  DO I=1,N
     Y(I,I)=1.
  ENDDO
  !
  CALL LUDCMP(A,N,INDX,D,FLAG)
  if (flag.ne.0) stop "Matinv: Matrix is singular!"
  !
  DO J=1,N
     CALL LUBKSB(A,N,INDX,Y(1,J))
  ENDDO
  A = Y
  !
  RETURN
END subroutine matinv
