SUBROUTINE MPINV (SIZE,N,M,MINNM,A,AINV,S,E,U,V, 1 WORK,IRANK,IERR) C C FUNCTION: CF CF Calculates the Moore-Penrose pseudo-inverse of the N by M CF matrix A, and stores the result in AINV. The singular CF values of A are returned in S. The left singular vectors CF are returned in U, and the right singular vectors are CF returned in V. CF C USAGE: CU CU CALL MPINV (SIZE,N,M,MINNM,A,AINV,S,E,U,V,WORK,IERR) CU C INPUTS: CI CI SIZE integer - leading dimension of all 2-dim. arrays. CI SIZE must be greater than or equal to max(N,M). CI CI N integer - first dimension of matrix A. CI CI M integer - second dimension of matrix A. CI CI MINNM integer - MINNM = min(N,M). CI CI A double precision - array of dimension SIZE by M CI containing the N by M matrix A. CI C OUTPUTS: CO CO AINV double precision - array of dimension SIZE by N CO containing the M by N matrix AINV which is the CO pseudo-inverse of A. CO CO S double precision - vector of dimension MINNM containing CO the singular values of A if IERR = 0; otherwise, see CO the description of LINPACK routine DSVDC. CO CO E double precision - vector of dimension MINNM used by CO LINPACK routine DSVDC; normally zero. CO CO U double precision - array of dimension SIZE by N CO containing the left singular vectors of A. CO CO V double precision - array of dimension size by M CO containing the right singular vectors of A. CO CO WORK double precision - vector of dimension N used by CO LINPACK routine DSVDC. CO CO IRANK integer - the calculated rank of the matrix A used CO in the construction of the pseudo-inverse. co CO IERR integer - error condition returned from dsvdc. If CO IERR is not zero, then AINV is not correct. CO C ALGORITHM: CA CA Computes the singular value decomposition of A and forms CA the sum of the outer products of the left and right singular CA vectors of the matrix A divided by the corresponding singular CA value, for all non-zero singular values. CA C MACHINE DEPENDENCIES: CM CM none CM C HISTORY: CH CH written by: Birdwell & Laub CH date: May 17, 1985 CH current version: 2.0 CH modifications: modifications so that it works CH under standard F77:bb:8-86. CH add minnm argument:bb:8-86. C ROUTINES CALLED: CC CC LINPACK: DSVDC, DROT, DAXPY, DDOT, DSCAL, DSWAP, DNRM2, DROTG CC C COMMON MEMORY USED: CM CM none CM C---------------------------------------------------------------------- C written for: The CASCADE Project C Oak Ridge National Laboratory C U.S. Department of Energy C Contract number DE-AC05-840R21400 C Subcontract number 37B-7685 S13 C Organization: The University of Tennessee C---------------------------------------------------------------------- C THIS SOFTWARE IS IN THE PUBLIC DOMAIN C NO RESTRICTIONS ON ITS USE ARE IMPLIED C---------------------------------------------------------------------- C C--global variables: C INTEGER SIZE INTEGER N INTEGER M INTEGER IRANK INTEGER IERR C DOUBLE PRECISION A(SIZE,M) DOUBLE PRECISION AINV(SIZE,N) DOUBLE PRECISION S(MINNM) DOUBLE PRECISION E(MINNM) DOUBLE PRECISION U(SIZE,N) DOUBLE PRECISION V(SIZE,M) DOUBLE PRECISION WORK(N) C C--local variables: C INTEGER I INTEGER J C DOUBLE PRECISION EPS C C--code: C C Call LINPACK routine to compute SVD. C CALL DSVDC (A,SIZE,N,M,S,E,U,SIZE,V,SIZE,WORK,21,IERR) IF (IERR .NE. 0) RETURN C C Compute machine epsilon. C EPS = 1.0D0 10 IF (EPS + 1.D0 .NE. 1.D0) THEN EPS = EPS / 2.D0 GO TO 10 END IF EPS = EPS * 2.D0 C EPS = EPS * S(1) C DO 20, J = 1, N DO 20, I = 1, M AINV(I,J) = 0.D0 20 CONTINUE C IRANK = 1 30 IF (S(IRANK) .GT. EPS .AND. IRANK .LE. MIN(N,M)) THEN DO 40, J = 1, N DO 40, I = 1, M AINV(I,J) = AINV(I,J) 1 + V(I,IRANK)*U(J,IRANK)/S(IRANK) 40 CONTINUE IRANK = IRANK + 1 GO TO 30 END IF C IRANK = IRANK - 1 C RETURN END