Line data    Source code

       1             : MODULE m_rinvgj
       2             : 
       3             :    CONTAINS
       4             : 
       5         148 :    SUBROUTINE rinvgj(ainv,aMat,arraydim,n)
       6             : !   ********************************************************************
       7             : !   *                                                                  *
       8             : !   *                      AINV = A**(-1)                              *
       9             : !   *                                                                  *
      10             : !   *  invert aMat using the GAUSS-JORDAN - algorithm                     *
      11             : !   *  the 1- matrix is not set up and use is made of its structure    *
      12             : !   *                                                                  *
      13             : !   *                    REAL*8 VERSION                                *
      14             : !   *                                                                  *
      15             : !   ********************************************************************
      16             : 
      17             :       IMPLICIT NONE
      18             : 
      19             :       INTEGER, INTENT (IN) :: arraydim,n
      20             : 
      21             :       REAL, INTENT (INOUT) :: aMat(arraydim,arraydim)
      22             :       REAL, INTENT (INOUT) :: ainv(arraydim,arraydim)
      23             : 
      24             :       REAL t, t1
      25             :       INTEGER icol, l, ll
      26             : 
      27         598 :       DO icol = 1,n
      28             : !                                               make A(ICOL,ICOL) = 1
      29         450 :          t1 = 1.0 / aMat(icol,icol)
      30         983 :          DO l = (icol+1),n
      31         983 :             aMat(icol,l) = aMat(icol,l)*t1
      32             :          END DO
      33         983 :          DO l = 1, (icol-1)
      34         983 :             ainv(icol,l) = ainv(icol,l)*t1
      35             :          END DO
      36         450 :          ainv(icol,icol) = t1
      37             : !                                    make A(LL,ICOL) = 0 for LL<>ICOL
      38        2114 :          DO ll = 1,n
      39        1966 :             IF (ll.NE.icol) THEN
      40        1066 :                t = aMat(ll,icol)
      41        2523 :                DO l = (icol+1),n
      42        2523 :                   aMat(ll,l) = aMat(ll,l) - aMat(icol,l)*t
      43             :                END DO
      44        2523 :                DO l = 1, (icol-1)
      45        2523 :                   ainv(ll,l) = ainv(ll,l) - ainv(icol,l)*t
      46             :                END DO
      47        1066 :                ainv(ll,icol) = -t1*t
      48             :             END IF
      49             :          END DO
      50             :       END DO
      51             : 
      52         148 :    END SUBROUTINE rinvgj
      53             : 
      54             : END MODULE m_rinvgj