C
C      ________________________________________________________
C     |                                                        |
C     |COMPUTE SINGULAR VALUE DECOMPOSITION OF A GENERAL MATRIX|
C     |      A = Q TIMES DIAGONAL MATRIX TIMES P TRANSPOSE     |
C     |                                                        |
C     |    INPUT:                                              |
C     |                                                        |
C     |         S     --ARRAY WITH AT LEAST N ELEMENTS         |
C     |                                                        |
C     |         LQ    --LEADING (ROW) DIMENSION OF ARRAY Q     |
C     |                                                        |
C     |         IQ    --AN INTEGER WHICH INDICATES WHICH COL-  |
C     |                 UMNS OF Q TO COMPUTE (= 0 MEANS NONE,  |
C     |                 = 1 MEANS FIRST L, = 2 MEANS LAST M-L, |
C     |                 = 3 MEANS ALL M WHERE L = MIN(M,N))    |
C     |                                                        |
C     |         LP    --LEADING (ROW) DIMENSION OF ARRAY P     |
C     |                                                        |
C     |         IP    --AN INTEGER (LIKE IQ) WHICH INDICATES   |
C     |                 WHICH COLUMNS OF P TO COMPUTE          |
C     |                                                        |
C     |         A     --ARRAY CONTAINING COEFFICIENT MATRIX    |
C     |                                                        |
C     |         LA    --LEADING (ROW) DIMENSION OF ARRAY A     |
C     |                                                        |
C     |         M     --ROW DIMENSION OF MATRIX STORED IN A    |
C     |                                                        |
C     |         N     --COLUMN DIMENSION OF MATRIX STORED IN A |
C     |                                                        |
C     |         W     --WORK ARRAY(LENGTH AT LEAST MAX(M,3L-1))|
C     |                                                        |
C     |    OUTPUT:                                             |
C     |                                                        |
C     |         Q     --Q FACTOR IN THE SINGULAR VALUE DECOMP. |
C     |                                                        |
C     |         S     --SINGULAR VALUES IN DECREASING ORDER    |
C     |                                                        |
C     |         P     --P FACTOR IN THE SINGULAR VALUE DECOMP. |
C     |                                                        |
C     |         A     --THE HOUSEHOLDER VECTORS USED IN THE    |
C     |                 REDUCTION PROCESS                      |
C     |                                                        |
C     |    NOTE:                                               |
C     |                                                        |
C     |         EITHER P OR Q CAN BE IDENTIFIED WITH A BUT NOT |
C     |         BOTH. WHEN EITHER P OR Q ARE IDENTIFIED WITH A,|
C     |         THE HOUSEHOLDER VECTORS IN A ARE DESTROYED. IF |
C     |         IQ = 2, Q MUST HAVE M COLUMNS EVEN THOUGH THE  |
C     |         OUTPUT ARRAY HAS JUST M-L COLUMNS. SIMILARLY IF|
C     |         IP = 2, P MUST HAVE N COLUMNS EVEN THOUGH THE  |
C     |         OUTPUT ARRAY HAS JUST N-L COLUMNS.             |
C     |                                                        |
C     |    BUILTIN FUNCTIONS: MIN0                             |
C     |    PACKAGE SUBROUTINES: BIDAG2,EIG3,FGV,HSR1-HSR5,SCL, |
C     |                         SFT,SINGB,SNG0,SNG1,SORT2      |
C     |________________________________________________________|
C
      SUBROUTINE SING(Q,LQ,IQ,S,P,LP,IP,A,LA,M,N,W)
      INTEGER I,IP,IQ,IU,JL,JR,L,LA,LP,LQ,M,N
      REAL A(LA,1),S(1),Q(LQ,1),P(LP,1),W(1)
      IF ( IQ .GE. 0 ) GOTO 20
10    WRITE(6,*) 'ERROR: INPUT PARAMETER IQ FOR SUBROUTINE SING'
      WRITE(6,*) 'EITHER LESS THAN 0 OR GREATER THAN 3'
      STOP
20    IF ( IQ .GT. 3 ) GOTO 10
      JL = 0
      IF ( IQ .EQ. 0 ) GOTO 30
      IF ( IQ .EQ. 2 ) GOTO 30
      JL = 1
30    IF ( IP .GE. 0 ) GOTO 50
40    WRITE(6,*) 'ERROR: INPUT PARAMETER IP FOR SUBROUTINE SING'
      WRITE(6,*) 'EITHER LESS THAN 0 OR GREATER THAN 3'
      STOP
50    IF ( IP .GT. 3 ) GOTO 40
      JR = 0
      IF ( IP .EQ. 0 ) GOTO 60
      IF ( IP .EQ. 2 ) GOTO 60
      JR = 1
60    CALL BIDAG2(S,W,Q,LQ,IQ,P,LP,IP,A,LA,M,N)
      L = MIN0(M,N)
      IF ( L .GT. 1 ) GOTO 80
      IF ( S(1) .GE. 0. ) RETURN
      S(1) = -S(1)
      IF ( JL .EQ. 0. ) RETURN
      DO 70 I = 1,M
70         Q(I,1) = -Q(I,1)
      RETURN
80    IU = 0
      IF ( M .GE. N ) GOTO 90
      IU = 1
90    CALL SINGB(S,L,W,IU,Q,LQ,M,JL,P,LP,N,JR,W(L),W(L+L))
      RETURN
      END