C....*...1.........2.........3.........4.........5.........6.........7.*.......8
C     DAPLUS   1/10/72
C
C     PURPOSE
C     COMPUTE THE MOORE-PENROSE G-INVERSE OF A MATRIX - A+
C     OBTAIN THE SINGULAR VALUE DECOMPOSITION OF A MATRIX
C
C     USAGE
C     CALL DAPLUS(A,N,M,AMP,U,S,V,IR,EPS,D1,D2,D3)
C
C     SUBROUTINES CALLED
C     DSVD
C
C     ARGUMENTS
C     A    - AN N BY M MATRIX STORED COLUMNWISE (STORAGE MODE OF 0)
C            ELEMENTS OF A ARE REAL*8
C     AMP  - MOORE-PENROSE G-INVERSE OF A
C            AN M BY N MATRIX STORED COLUMNWISE (STORAGE MODE OF 0)
C            ELEMENTS ARE REAL*8
C     N    - NUMBER OF ROWS IN A  (ON OUTPUT THE NO. OF COLS. OF A+)
C     M    - NUMBER OF COLUMNS IN A (ON OUTPUT THE NO. OF ROWS OF A+)
C            N MUST BE GREATER THAN OR EQUAL TO M
C     U    - AN N BY M MATRIX STORED COLUMNWISE (STORAGE MODE OF 0)
C            ELEMENTS OF U ARE REAL*8
C     S    - A DIAGONAL MATRIX STORED AS AN M-VECTOR (STORAGE MODE OF 2)
C            ELEMENTS OF S ARE REAL*8
C     V    - AN M BY M MATRIX STORED COLUMNWISE (STORAGE MODE OF 0)
C            ELEMENTS OF V ARE REAL*8
C     IR   - COMPUTED RANK OF A
C            INTEGER
C     EPS  - VALUE USED TO TEST THE SINGULAR VALUES OF A AND DETERMINE
C            THE RANK OF A.  A REASONABLE VALUE IS EPS = 1.D-13
C            REAL*8
C     D1   - AN M VECTOR USED FOR WORKSPACE.
C     D2   - AN M VECTOR USED FOR WORKSPACE.
C     D3   - AN M VECTOR USED FOR WORKSPACE.
C            ELEMENTS OF D1,D2,D3 ARE REAL*8
C
C     REMARKS
C     A = U*S*(V-TRANSPOSE)
C     (U-TRANSPOSE)*U = (V-TRANSPOSE)*V = V*(V-TRANSPOSE) = I
C     A+ = V*(S+)*(U-TRANSPOSE)
C     S+(I) = 1/S(I) IF S(I).GT.0 AND S+(I)=0 OTHERWISE
C     IF A CAN BE DESTROYED THE FOLLOWING USAGE WILL CONSERVE CORE
C     CALL DAPLUS(A,N,M,AMP,A,S,V,IR,EPS,D1,D2,D3)
C     THE MATRIX A WILL CONTAIN U ON RETURN
C
C     REFERENCE
C     BUSINGER,P.A. AND GOLUB,G.H., SINGULAR VALUE DECOMPOSITION OF A
C     COMPLEX MATRIX. COMMUNICATIONS OF THE ACM 12: 564-565 (OCT. 1969)
C
      SUBROUTINE DAPLUS(A,N,M,AMP,U,S,V,IR,EPS,D1,D2,D3)
      IMPLICIT REAL*8 (A-H,O-Z)
      save
      REAL*8 A(N,M),AMP(M,N),U(N,M),S(M),V(M,M),D1(M),D2(M),D3(M)
      CALL Z1DAPLUS(A,N,M,AMP)
      CALL DSVD(AMP,N,M,N,M,0,M,M,S,U,V,D1,D2,D3)
C
C     RANK DETERMINATION
      IR=0
      TEST=S(1)*EPS
      DO 10 J=1,M
      IF(S(J).GT.TEST) IR=IR+1
      D1(J)=0.D0
10    IF(S(J).GT.TEST) D1(J)=1.D0/S(J)
C
      CALL Z2DAPLUS(AMP,M,N,V,D1,U,D2)
      RETURN
      END
      SUBROUTINE Z1DAPLUS(A,N,M,AMP)
      implicit real*8 (a-h,o-z)
      save
      REAL*8 A(1),AMP(1)
      L=M*N
      DO 10 I=1,L
10    AMP(I)=A(I)
      RETURN
      END
      SUBROUTINE Z2DAPLUS(A,M,N,V,D1,U,D2)
      implicit real*8 (a-h,o-z)
      save
      REAL*8 A(M,N), V(M,M),D1(M), U(N,M), D2(M)
      DO 30 I=1,M
      DO 10 K=1,M
10    D2(K)=V(I,K)*D1(K)
      DO 20 J=1,N
      A(I,J)=0.D0
      DO 20 K=1,M
20    A(I,J)=A(I,J)+D2(K)*U(J,K)
30    CONTINUE
      RETURN
      END