C
C      ________________________________________________________
C     |                                                        |
C     | GIVEN THE SINGULAR VALUE DECOMPOSITION QDP TRANSPOSE OF|
C     | A GENERAL MATRIX, COMPUTE ITS REGULARIZED PSEUDOINVERSE|
C     |                                                        |
C     |    INPUT:                                              |
C     |                                                        |
C     |         LA    --LEADING (ROW) DIMENSION OF ARRAY A     |
C     |                                                        |
C     |         Q     --FIRST FACTOR IN THE SVD DECOMPOSITION  |
C     |                                                        |
C     |         LQ    --LEADING (ROW) DIMENSION OF ARRAY Q     |
C     |                                                        |
C     |         MQ    --ROW DIMENSION OF MATRIX STORED IN Q    |
C     |                                                        |
C     |         D     --SINGULAR VALUES                        |
C     |                                                        |
C     |         P     --LAST FACTOR IN THE SVD DECOMPOSITION   |
C     |                                                        |
C     |         LP    --LEADING (ROW) DIMENSION OF ARRAY P     |
C     |                                                        |
C     |         MP    --ROW DIMENSION OF MATRIX STORED IN P    |
C     |                                                        |
C     |         R     --REGULARIZATION PARAMETER               |
C     |                                                        |
C     |    OUTPUT:                                             |
C     |                                                        |
C     |         A     --THE REGULARIZED PSEUDOINVERSE          |
C     |                                                        |
C     |    BUILTIN FUNCTIONS: ABS,MIN0,SQRT                    |
C     |________________________________________________________|
C
      SUBROUTINE PSEUDO(A,LA,Q,LQ,MQ,D,P,LP,MP,R)
      INTEGER I,J,K,LA,LP,LQ,MP,MQ,N
      REAL A(LA,1),D(1),Q(LQ,1),P(LP,1),R,S,T,Y,Z
      N = MIN0(MQ,MP)
      Z = ABS(R)
      Y = SQRT(Z)
      DO 10 J = 1,MQ
           DO 10 I = 1,MP
10              A(I,J) = 0.
      DO 50 K = 1,N
           IF ( D(K) .EQ. 0. ) GOTO 50
           IF ( D(K) .GT. Y ) GOTO 20
           S = D(K)/(D(K)**2+Z)
           GOTO 30
20         S = 1./(D(K)+Z/D(K))
30         DO 40 J = 1,MQ
                T = S*Q(J,K)
                DO 40 I = 1,MP
40                   A(I,J) = A(I,J) + T*P(I,K)
50    CONTINUE
      RETURN
      END