invmtx.f

Go to the documentation of this file.

!                   *****************
                    SUBROUTINE invmtx
!                   *****************
!
     &(am,bm,np)
!
!***********************************************************************
! BIEF   V6P1                                   21/08/2010
!***********************************************************************
!
!brief    INVERTS A MATRIX OF NP BY NP.
!+
!+            BM IS THE INVERSION OF AM.
!
!note     THIS SUBROUTINE CALLS LUDCMP AND LUBKSB, WHICH WERE
!+         COPIED FROM "NUMERIC RECIPES" -- A WELL-KNOWN BOOK.
!
!history  CHUN WANG
!+        28/07/2006
!+        V5P7
!+
!
!history  N.DURAND (HRW), S.E.BOURBAN (HRW)
!+        13/07/2010
!+        V6P0
!+   Translation of French comments within the FORTRAN sources into
!+   English comments
!
!history  N.DURAND (HRW), S.E.BOURBAN (HRW)
!+        21/08/2010
!+        V6P0
!+   Creation of DOXYGEN tags for automated documentation and
!+   cross-referencing of the FORTRAN sources
!
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!| AM             |-->| MATRIX TO BE INVERTED
!| BM             |<--| INVERTED MATRIX
!| NP             |-->| RANK OF MATRICES AM AND BM.
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!
      USE declarations_special
      IMPLICIT NONE
!
!+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
!
      INTEGER, INTENT(IN)             :: NP
      DOUBLE PRECISION, INTENT(INOUT) :: AM(np,np),BM(np,np)
!
!+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
!
      INTEGER INDX(500),N,I,J
!
!-----------------------------------------------------------------------
!
      IF(np.GT.500) THEN
        WRITE(lu,*) 'NP MUST BE LESS THAN 500 IN INVMTX'
        CALL plante(1)
        stop
      ENDIF
!
      n = np
      DO i=1,n !SET UP IDENTITY MATRIX
        DO j=1,n
          bm(i,j)=0.d0
        ENDDO
        bm(i,i)=1.d0
      ENDDO
!
!     DECOMPOSES THE MATRIX JUST ONCE
!
      CALL ludcmp(am,n,np,indx)
!
!     FINDS INVERSE BY COLUMNS
!
      DO j=1,n
        CALL lubksb(am,n,np,indx,bm(1,j))
      ENDDO
!
!-----------------------------------------------------------------------
!
      RETURN
      END