ludcmp.f

Go to the documentation of this file.
00001 C:\opentelemac\v7p0\sources\utils\bief\ludcmp.f
00002 !
00067                      SUBROUTINE LUDCMP
00068 !                    *****************
00069 !
00070      &(A,N,NP,INDX)
00071 !
00072 !***********************************************************************
00073 ! BIEF   V6P2                                   21/08/2010
00074 !***********************************************************************
00075 !
00076 !
00077 !
00078 !
00079 !
00080 !
00081 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
00082 !| A              |-->| MATRIX OF THE SYSTEM
00083 !| INDX           |-->| ADDRESS IN RIGHT-HAND SIDE
00084 !| N              |-->| SIZE OF B
00085 !| NP             |-->| RANK OF MATRIX A
00086 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
00087 !
00088       IMPLICIT NONE
00089       INTEGER LNG,LU
00090       COMMON/INFO/LNG,LU
00091 !
00092 !+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
00093 !
00094       INTEGER, INTENT(IN)             :: N,NP
00095       INTEGER, INTENT(INOUT)          :: INDX(N)
00096       DOUBLE PRECISION, INTENT(INOUT) :: A(NP,NP)
00097 !
00098 !+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
00099 !
00100       DOUBLE PRECISION D
00101       INTEGER I,IMAX,J,K
00102       DOUBLE PRECISION AAMAX,DUM,XSOM,VV(500)
00103       DOUBLE PRECISION, PARAMETER:: TINY=1.D-20
00104 !
00105 !------------------------------------------------------------------------
00106 !
00107       D=1.D0 ! NO ROW INTERCHANGES YET
00108 !
00109 !     LOOP OVER ROWS TO GET THE IMPLICIT SCALING INFORMATION
00110 !
00111       DO I=1,N
00112         AAMAX=0.D0
00113         DO J=1,N
00114           IF(ABS(A(I,J)).GT.AAMAX) AAMAX=ABS(A(I,J))
00115         ENDDO
00116         IF(AAMAX.LT.TINY) THEN
00117           WRITE(LU,*) 'SINGULAR MATRIX IN LUDCMP'
00118           CALL PLANTE(1)
00119           STOP
00120         ENDIF
00121         VV(I)=1.D0/AAMAX ! SAVES THE SCALING
00122       ENDDO
00123 !
00124 !     LOOP OVER COLUMNS OF CROUT'S METHOD
00125 !
00126       DO J=1,N
00127         DO I=1,J-1 ! EQUATION (2.3.12) EXCEPT FOR I = J
00128           XSOM=A(I,J)
00129           DO K=1,I-1
00130              XSOM=XSOM-A(I,K)*A(K,J)
00131           ENDDO
00132           A(I,J)=XSOM
00133         ENDDO
00134         AAMAX=0.D0 ! INITIALISES FOR THE SEARCH OF LARGEST PIVOT ELEMENT
00135         DO I=J,N   ! THIS IS I = J OF EQUATION (2.3.12) AND
00136                   ! I = J +1 : : : N OF EQUATION (2.3.13)
00137           XSOM=A(I,J)
00138           DO K=1,J-1
00139             XSOM=XSOM-A(I,K)*A(K,J)
00140           ENDDO
00141           A(I,J)=XSOM
00142           DUM=VV(I)*ABS(XSOM) ! FIGURE OF MERIT FOR THE PIVOT
00143           IF (DUM.GE.AAMAX) THEN ! IS IT BETTER THAN THE BEST SO FAR?
00144             IMAX=I
00145             AAMAX=DUM
00146           ENDIF
00147         ENDDO
00148         IF (J.NE.IMAX) THEN ! NEEDS TO INTERCHANGE ROWS?
00149           DO K=1,N
00150             DUM=A(IMAX,K)
00151             A(IMAX,K)=A(J,K)
00152             A(J,K)=DUM
00153           ENDDO
00154           D=-D !...AND CHANGES THE PARITY OF D
00155           VV(IMAX)=VV(J) ! ALSO INTERCHANGES THE SCALE FACTOR
00156         ENDIF
00157         INDX(J)=IMAX
00158         IF(ABS(A(J,J)).LT.TINY) A(J,J)=TINY
00159 !
00160 !  IF THE PIVOT ELEMENT IS 0 THE MATRIX IS SINGULAR (AT LEAST TO THE
00161 !  PRECISION OF THE ALGORITHM)
00162 !
00163         IF(J.NE.N) THEN ! DIVIDES BY THE PIVOT ELEMENT
00164           DUM=1.D0/A(J,J)
00165           DO I=J+1,N
00166             A(I,J)=A(I,J)*DUM
00167           ENDDO
00168         ENDIF
00169 !
00170       ENDDO
00171 !
00172 !------------------------------------------------------------------------
00173 !
00174       RETURN
00175       END

Generated on Fri Aug 31 2013 18:12:58 by S.E.Bourban (HRW) using doxygen 1.7.0