C This file: http://ftp.aset.psu.edu/pub/ger/fortran/hdk/ginv.for
C
      DOUBLE PRECISION A(3,3),U(3,3),AFLAG(3),ATEMP(3)
      INTEGER MR,NR,NC, I,J
      NR=3
      MR=3
      NC=3
      DO I=1,NR
        DO J=1,NC
          A(I,J)=0
        END DO
        A(I,I)=I
      END DO
      DO I=1,NR
       WRITE(*,*) (A(I,J),J=1,NC)
      END DO
      WRITE(*,*) ' '
      CALL GINV2(A,U,AFLAG,ATEMP,3,3,3)
      DO I=1,NC
       WRITE(*,*) (A(J,I),J=1,NR)
      END DO

      STOP
      END
      SUBROUTINE GINV2 (A,U,AFLAG,ATEMP,MR,NR,NC)
C A Simple Algorithm for Computing the Generalized Inverse of a Matrix
C  by B. RUST, W. R. BURRUS AND C. SCHNEEBERGER
C  CACM 9(5):381-387 (May, 1966)
C
C THIS ROUTINE CALCULATES THF GENERALIZED INVERSE OF A
C AND STORES THE TRANSPOSE OF IT IN A.
C MR -> FIRST DIMENSION OF A.
C NR -> NO. OF ROWS IN A
C NC -> NO, OF COLUMNS IN A
C U IS THE BOOKKEEPING MATRIX.
C AFLAG AND ATEMP ARE TEMPORRARY WORKING STORAGE.
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      INTEGER MR,NR,NC,I,J,K,L
      DIMENSION A(MR,NC) ,U(NC,NC) ,AFLAG(NC) ,ATEMP(NC)
      DO 10 I = 1,NC
        DO 5 J = 1,NC
5       U (I,J) = 0.0
10    U(I,I) = 1.0
      FAC = DOT(MR,NR,A,1,1,NC)
      FAC= 1.0/SQRT(FAC)
      DO 15 I = 1,NR
15    A(I,1) = A(I,1) * FAC
      DO 20 I = 1,NC
20    U (I,1) = U(I,1)*FAC
      AFLAG(1) = 1.0
C
C DEPENDENT COL TOLERANCE FOR N BIT FLOATING POINT FRACTION
C
      N = 27
      TOL = (10. * 0.5**N)**2
      WRITE(*,*) 'TOL=',TOL
      DO 100 J = 2,NC
        DOT1 = DOT(MR,NR,A,J,J,NC)
        JM1=J-1
         DO 50 L=1,2
          DO 30 K=1,JM1
30        ATEMP(K) = DOT(MR,NR,A,J,K,NC)
          DO 45 K=1,JM1
            DO 35 I = 1,NR
35          A(I,J) = A(I,J)-ATEMP(K)*A(I,K)*AFLAG(K)
            DO 40 I = 1,NC
40          U(I,J) = U(I,J)-ATEMP(K)*U(I,K)
45        CONTINUE
50       CONTINUE
        DOT2 = DOT(MR,NR,A,J,J,NC)
        IF((DOT2/DOT1) - TOL) 55,55,70
55      DO 60 I=I,JM1
          ATEMP (I)=0.0
          DO 60 K=1,I
60      ATEMP(I) = ATEMP(I) + U(K,I)*U(K,J)
        DO 65 I = 1,NR
          A( I,J)=0.0
          DO 65 K=I,JM1
65      A(I,J) = A(I,J) - A (I,K)*ATEMP(K)*AFLAG(K)
        AFLAG(J) = 0.0
        FAC = DOT(NC,NC,U,J,J,NC)
        FAC= 1.0/SQRT(FAC)
        GO TO 75
70      AFLAG(J) = 1.0
        FAC=1.0/SQRT(DOT2)
75      DO 80 I = 1,NR
80      A(I,J) = A(I,J)*FAC
        DO 85 I = 1,NC
85      U(I,J) = U(I,J)*FAC
100   CONTINUE
      DO 130 J=1,NC
      DO 130 I=1,NR
        FAC = 0.0
        DO 120 K = J,NC
120     FAC=FAC+A(I,K)*U(J,K)
130   A(I,J) = FAC
      RETURN
      END
      FUNCTION DOT(MR,NR,A,JC,KC,NC)
C  COMPUTES THE INNER RPODUCT OF COLUMNS JC AND KC
C  OF MATRIX A.
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      DIMENSION A(MR,NC)
      DOT=0.0
      DO 5 I = 1,NR
5     DOT = DOT + A(I,JC)*A(I,KC)
      RETURN
      END