
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
  
