Chimie Théorique » OpenPath

SUBROUTINE GINVSE(A, MA, NA, AIN, NAIN, M, N, TOL, V, NV, FAIL, NOUT)

  Use Io_module

  ! INVERSE.OF.A = V * INVERSE.OF.S * TRANSPOSE.OF.U

  ! VIA SINGULAR VALUE DECOMPOSITION USING ALGORITHMS 1 AND 2 OF

  ! NASH J C, COMPACT NUMERICAL METHODS FOR COMPUTERS, 1979,

  !    ADAM HILGER, BRISTOL OR HALSTED PRESS, N.Y.

  !

  !

  ! CALLING SEQUENCE

  !

  !    CALL GINVSE (A,MA,NA,AIN,NAIN,M,N,TOL,V,NV,FAIL,NOUT)

  !

  ! A      = MATRIX OF SIZE M BY N TO BE 'INVERTED'

  !          A IS DESTROYED BY THIS ROUTINE AND BECOMES THE

  !          U MATRIX OF THE SINGULAR VALUE DECOMPOSITION

  ! MA     = ROW DIMENSION OF ARRAY A

  !          UNCHANGED BY THIS SUB-PROGRAM

  ! NA     = COLUMN DIMENSION OF ARRAY A

  !          UNCHANGED BY THIS SUB-PROGRAM

  ! AIN    = COMPUTED 'INVERSE' (N ROWS BY M COLS.)

  ! NAIN   = ROW DIMENSION OF ARRAY AIN

  !          ( COLUMN DIMENSION ASSUMED TO BE AT LEAST M )

  !          UNCHANGED BY THIS SUB-PROGRAM

  ! M      = NUMBER OF ROWS IN MATRIX A AND

  !          THE NUMBER OF COLUMNS IN MATRIX AIN

  ! N      = NUMBER OF COLUMNS IN MATRIX A AND

  !          THE NUMBER OF ROWS IN MATRIX AIN

  !          UNCHANGED BY THIS SUB-PROGRAM

  ! TOL    = MACHINE PRECISION FOR COMPUTING ENVIRONMENT

  !          UNCHANGED BY THIS SUB-PROGRAM

  ! V      = WORK MATRIX WHICH BECOMES THE V MATRIX (N BY N)

  !          OF THE SINGULAR VALUE DECOMPOSITION

  ! NV     = DIMENSIONS OF V  (BOTH ROWS AND COLUMNS)

  !          UNCHANGED BY THIS SUB-PROGRAM

  ! FAIL   = ERROR FLAG - .TRUE. IMPLIES FAILURE OF GINVSE

  ! NOUT   = OUTPUT CHANNEL NUMBER

  !          UNCHANGED BY THIS SUB-PROGRAM

  !

  ! NOTE.  THIS ROUTINE IS NOT DESIGNED TO PRODUCE A 'GOOD'

  !     GENERALIZED INVERSE. IT IS MEANT ONLY TO FURNISH TEST

  !     DATA FOR THE ROUTINES  GMATX, ZIELKE, AND PTST

  !

  INTEGER(KINT), INTENT(IN) :: MA, NA, NAIN,N, M,NV, NOUT

  INTEGER(KINT) :: ICOUNT, NM1, I, J, JP1, K,SWEEP, LIMIT

  LOGICAL, INTENT(INOUT) :: FAIL

  REAL(KREAL), INTENT(IN) :: TOL

  REAL(KREAL) :: P, Q, R, VV, C, S, TEMP

  REAL(KREAL), INTENT(INOUT) :: V(NV,NV)

  REAL(KREAL), INTENT(INOUT) :: AIN(NAIN,M)

  REAL(KREAL), INTENT(INOUT) :: A(MA,NA)

  LOGICAL :: debug

  INTERFACE

     function valid(string) result (isValid)

       CHARACTER(*), intent(in) :: string

       logical                  :: isValid

     END function VALID

  END INTERFACE

  debug=valid('Ginvse')

  if (debug) WRITE(*,*) '=======Entering Ginvse======='

  !Print *, 'A='

  DO J=1,MA

     !WRITE(*,'(50(1X,F10.2))') AIN(:,J)

     !Print *, A(:,J)

  END DO

  !Print *, 'End of A'

  !

  ! INITIALLY SET FAIL=.FALSE. TO IMPLY CORRECT EXECUTION

  !

  FAIL = .FALSE.

  !

  ! TESTS ON VALIDITY OF DIMENSIONS

  !

  IF (M.LT.2 .OR. MA.LT.2 .OR. N.LT.2 .OR. NA.LT.2) GO TO 230

  IF (N.GT.NA) GO TO 210

  IF (M.GT.MA) GO TO 220

  !

  ! N MUST NOT EXCEED M, OTHERWISE GINVSE WILL FAIL

  !

  IF (N.GT.M) GO TO 200

  !

  ! SWEEP COUNTER  INITIALIZED TO ZERO

  !

  SWEEP = 0

  !

  ! SET SWEEP LIMIT

  !

  LIMIT = MAX0(N,6)

  !

  ! MAX0(N,6) WAS  CHOSEN FROM EXPERIENCE

  !

  !

  ! V(I,J) INITIALLY N BY N IDENTITY

  !

  DO 20 I=1,N

     DO 10 J=1,N

        V(I,J) = 0.0

10      CONTINUE

        V(I,I) = 1.0

20      CONTINUE

        !

        ! INITIALIZE ROTATION COUNTER (COUNTS DOWN TO 0)

        !

30      ICOUNT = N*(N-1)/2

        !

        ! COUNT SWEEP

        !

        SWEEP = SWEEP + 1

        NM1 = N - 1

        DO 110 J=1,NM1

           JP1 = J + 1

           DO 100 K=JP1,N

              P = 0.

              Q = 0.

              R = 0.

              DO 40 I=1,M

                 !

           ! TEST FOR AND AVOID UNDERFLOW

           ! NOT NEEDED FOR MACHINES WHICH UNDERFLOW TO ZERO WITHOUT

           ! ERROR MESSAGE

           !

           IF (ABS(A(I,J)).GE.TOL) Q = Q + A(I,J)*A(I,J)

           !Print *, '1 Q=',Q

           !Print *, 'A(',I,J,')=',A(I,J)

           IF (ABS(A(I,K)).GE.TOL) R = R + A(I,K)*A(I,K)

           IF (ABS(A(I,J)/TOL)*ABS(A(I,K)/TOL).GE.1.0) P = P + A(I,J)*A(I,K)

40         CONTINUE

           !Print *, '2 Q=',Q

           IF (Q.GE.R) GO TO 50

           C = 0.

           S = 1.

           GO TO 60

50         IF (ABS(Q).LT.TOL .AND. ABS(R).LT.TOL) GO TO 90

           IF (R.EQ.0.0) GO TO 90

           IF ((P/Q)*(P/R).LT.TOL) GO TO 90

           !

           ! CALCULATE THE SINE AND COSINE OF THE ANGLE OF ROTATION

           !

           !Print *, '3 Q=',Q

           !Print *, 'R=',R

           Q = Q - R

           VV = SQRT(4.*P*P+Q*Q)

           !Print *, 'VV=',VV

           !Print *, 'P=',P

           !Print *, 'Q=',Q

           C = SQRT((VV+Q)/(2.*VV))

           !Print *, 'C=', C

           S = P/(VV*C)

           !Print *, 'S=', S

           !

           ! APPLY THE ROTATION TO A

           !

60         DO 70 I=1,M

              R = A(I,J)

              A(I,J) = R*C + A(I,K)*S

              A(I,K) = -R*S + A(I,K)*C

70            CONTINUE

              !if (debug) WRITE(*,*) '=======After rotation of A======='

              !

              ! APPLY THE ROTATION TO V

              !

              DO 80 I=1,N

                 R = V(I,J)

                 V(I,J) = R*C + V(I,K)*S

                 !Print *, 'V(I,J)=', V(I,J)

                 !Print *, 'R=', R

                 V(I,K) = -R*S + V(I,K)*C

                 !Print *, 'V(I,K)=', V(I,K)

80               CONTINUE

                 GO TO 100

90               ICOUNT = ICOUNT - 1

100              CONTINUE

110              CONTINUE

                 !if (debug) WRITE(*,*) '=======After rotation of V======='

                 !

                 ! PRINT THE NUMBER OF SWEEPS AND ROTATIONS

                 !

                 IF (NOUT.GT.0) WRITE (NOUT,9997) SWEEP, ICOUNT

                 !

                 ! CHECK NUMBER OF SWEEPS AND ROTATIONS (TERMINATION TEST)

                 !

                 IF (ICOUNT.GT.0 .AND. SWEEP.LT.LIMIT) GO TO 30

                 IF (SWEEP.GE.LIMIT .AND. NOUT.GT.0) WRITE (NOUT,9996)

                 DO 160 J=1,N

                    Q = 0.

                    DO 120 I=1,M

                       Q = Q + A(I,J)*A(I,J)

120                    CONTINUE

                       !

                       ! ARBITRARY RANK DECISION

                       !

                       IF (J.EQ.1) VV = 1.0E-3*SQRT(Q)

                       IF (SQRT(Q).LE.VV) GO TO 140

                       DO 130 I=1,M

                          A(I,J) = A(I,J)/Q

130                       CONTINUE

                          GO TO 160

140                       DO 150 I=1,M

                             A(I,J) = 0.0

150                          CONTINUE

160                          CONTINUE

                                   !if (debug) WRITE(*,*) '=======After termination test======='

                                   DO 190 I=1,N

                                      DO 180 J=1,M

                                         TEMP = 0.0

                                         DO 170 K=1,N

                                            TEMP = TEMP + V(I,K)*A(J,K)

                                            !Print *, 'V(I,K)*A(J,K)=', V(I,K)*A(J,K)

                                            !Print *, 'V(I,K)=', V(I,K)

                                            !Print *, 'A(J,K)=', A(J,K)

170    CONTINUE

       AIN(I,J) = TEMP

       !Print *, 'AIN(I,J)=', AIN(I,J)

180    CONTINUE

190    CONTINUE

       RETURN

       !Print *, 'AIN='

       !DO J=1,NAIN

       !WRITE(*,'(50(1X,F10.2))') AIN(:,J)

       !Print *, AIN(:,J)

       !END DO

       !Print *, 'End of AIN'

       !

       ! IF AN ERROR OCCURED THE APPROPRIATE MESSAGE IS PRINTED

       !    AND FAIL IS SET TO  .TRUE.

       !

200    IF (NOUT.GT.0) WRITE (NOUT,9995)

       FAIL = .TRUE.

       RETURN

210    IF (NOUT.GT.0) WRITE (NOUT,9999) N, NA

       FAIL = .TRUE.

       RETURN

220    IF (NOUT.GT.0) WRITE (NOUT,9998) M, MA

       FAIL = .TRUE.

       RETURN

230    IF (NOUT.GT.0) WRITE (NOUT,9994) M, MA, N, NA

       FAIL = .TRUE.

       RETURN

       !

       ! * * * * * * * * * * * * * * * * * * *

                                            !

9999    FORMAT (11H0**ERROR** , I5, 24H EXCEEDS ALLOWABLE VALUE,23H FOR N IN GINVSE   NA= , I5, 5(/))

9998        FORMAT (11H0**ERROR** , I5, 24H EXCEEDS ALLOWABLE VALUE,23H FOR M IN GINVSE   MA= , I5, 5(/))

9997        FORMAT (8H SWEEP =, I4, 2H  , I4, 19H JACOBI ROTATIONS P,8HERFORMED)

9996        FORMAT (21H SWEEP LIMIT REACHED )

9995        FORMAT (38H0**ERROR** N GREATER THAN M IS INVALID,21H AND GINVSE WILL FAIL, 5(/))

9994        FORMAT (46H0**ERROR** SOME DIMENSIONS ARE LESS THAN TWO I,8HN GINVSE &

                 /10X, 3HM =, I5, 5H MA =,I5, 4H N =, I5,5H NA =, I5, 5(/))

            !

            ! * * * END OF SUBROUTINE GINVSE * * *

            !

          END SUBROUTINE GINVSE