ludcmp.f

Go to the documentation of this file.

!                   *****************
                    SUBROUTINE ludcmp
!                   *****************
!
     &(a,n,np,indx)
!
!***********************************************************************
! BIEF   V6P2                                   21/08/2010
!***********************************************************************
!
!brief    GIVEN A MATRIX A(1:N,1:N), WITH PHYSICAL DIMENSION NP
!+                BY NP, THIS ROUTINE REPLACES IT BY THE LU FACTORISATION
!+                OF A ROWWISE PERMUTATION OF ITSELF. A AND N ARE INPUT.
!+                A IS OUTPUT, ARRANGED AS IN EQUATION (2.3.14) ABOVE;
!+                INDX(1:N) IS AN OUTPUT VECTOR THAT RECORDS THE ROW
!+                PERMUTATION EFFECTED BY THE PARTIAL PIVOTING; D IS
!+                OUTPUT AS SIGMA1 DEPENDING ON WHETHER THE NUMBER OF
!+                ROW INTERCHANGES WAS EVEN OR ODD, RESPECTIVELY. THIS
!+                SUBROUTINE IS USED IN COMBINATION WITH LUBKSB TO SOLVE
!+                LINEAR EQUATIONS OR INVERT A MATRIX.
!
!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
!
!history  J-M HERVOUET (LNHE)
!+        25/07/2012
!+        V6P2
!+   Correction of one test of division by 0.
!
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!| A              |-->| MATRIX OF THE SYSTEM
!| INDX           |-->| ADDRESS IN RIGHT-HAND SIDE
!| N              |-->| SIZE OF B
!| NP             |-->| RANK OF MATRIX A
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!
      USE declarations_special
      IMPLICIT NONE
!
!+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
!
      INTEGER, INTENT(IN)             :: N,NP
      INTEGER, INTENT(INOUT)          :: INDX(n)
      DOUBLE PRECISION, INTENT(INOUT) :: A(np,np)
!
!+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
!
      DOUBLE PRECISION D
      INTEGER I,IMAX,J,K
      DOUBLE PRECISION AAMAX,DUM,XSOM,VV(500)
      DOUBLE PRECISION, PARAMETER:: CHOUIA=1.d-20
!
!------------------------------------------------------------------------
!
      d=1.d0 ! NO ROW INTERCHANGES YET
!
!     LOOP OVER ROWS TO GET THE IMPLICIT SCALING INFORMATION
!
      DO i=1,n
        aamax=0.d0
        DO j=1,n
          IF(abs(a(i,j)).GT.aamax) aamax=abs(a(i,j))
        ENDDO
        IF(aamax.LT.chouia) THEN
          WRITE(lu,*) 'SINGULAR MATRIX IN LUDCMP'
          CALL plante(1)
          stop
        ENDIF
        vv(i)=1.d0/aamax ! SAVES THE SCALING
      ENDDO
!
!     LOOP OVER COLUMNS OF CROUT'S METHOD
!
      DO j=1,n
        DO i=1,j-1 ! EQUATION (2.3.12) EXCEPT FOR I = J
          xsom=a(i,j)
          DO k=1,i-1
            xsom=xsom-a(i,k)*a(k,j)
          ENDDO
          a(i,j)=xsom
        ENDDO
        aamax=0.d0 ! INITIALISES FOR THE SEARCH OF LARGEST PIVOT ELEMENT
        DO i=j,n   ! THIS IS I = J OF EQUATION (2.3.12) AND
                  ! I = J +1 : : : N OF EQUATION (2.3.13)
          xsom=a(i,j)
          DO k=1,j-1
            xsom=xsom-a(i,k)*a(k,j)
          ENDDO
          a(i,j)=xsom
          dum=vv(i)*abs(xsom) ! FIGURE OF MERIT FOR THE PIVOT
          IF (dum.GE.aamax) THEN ! IS IT BETTER THAN THE BEST SO FAR?
            imax=i
            aamax=dum
          ENDIF
        ENDDO
        IF (j.NE.imax) THEN ! NEEDS TO INTERCHANGE ROWS?
          DO k=1,n
            dum=a(imax,k)
            a(imax,k)=a(j,k)
            a(j,k)=dum
          ENDDO
          d=-d !...AND CHANGES THE PARITY OF D
          vv(imax)=vv(j) ! ALSO INTERCHANGES THE SCALE FACTOR
        ENDIF
        indx(j)=imax
        IF(abs(a(j,j)).LT.chouia) a(j,j)=chouia
!
!  IF THE PIVOT ELEMENT IS 0 THE MATRIX IS SINGULAR (AT LEAST TO THE
!  PRECISION OF THE ALGORITHM)
!
        IF(j.NE.n) THEN ! DIVIDES BY THE PIVOT ELEMENT
          dum=1.d0/a(j,j)
          DO i=j+1,n
            a(i,j)=a(i,j)*dum
          ENDDO
        ENDIF
!
      ENDDO
!
!------------------------------------------------------------------------
!
      RETURN
      END