C
C*    Group  Linear Solver subroutines (Code DECCON/SOLCON)
C
      SUBROUTINE BVDECC (A,NROW,NCOL,MCON,M,N,IRANK,COND,D,
     1                                            PIVOT,KRED,AH,V)
C----------------------------------------------------------------------
C
C     CONSTRAINED QR-DECOMPOSITION OF (M,N)-MATRIX A
C     FIRST MCON ROWS BELONG TO EQUALITY CONSTRAINTS
C
C
C  REFERENCES:
C     1. P.DEUFLHARD, V.APOSTOLESCU:
C        AN UNDERRELAXED GAUSS-NEWTON METHOD FOR EQUALITY CONSTRAINED
C        NONLINEAR LEAST SQUARES PROBLEMS.
C        LECTURE NOTES CONTROL INFORM. SCI. VOL. 7, P. 22-32 (1978)
C
C     2. P.DEUFLHARD, W.SAUTTER:
C        ON RANK-DEFICIENT PSEUDOINVERSES.
C        J. LIN. ALG. APPL. VOL. 29, P. 91-111 (1980)
C
C*********************************************************************
C
C     TO BE USED IN CONNECTION WITH SUBROUTINE SOLCON (BVSOLC)
C
C *********** DATE OF LATEST CHANGE ****** MAY 10,82 ******************
C
C     INPUT PARAMETERS (* MARKS INOUT PARAMETERS)
C     -----------------------------------------------
C
C
C      * A(NROW,NCOL)  INPUT MATRIX
C                      A(M,N) CONTAINS ACTUAL INPUT
C        NROW          DECLARED NUMBER OF ROWS OF A AND AH
C        NCOL          DECLARED NUMBER OF COLUMNS OF A AND AH
C     (*)MCON          NUMBER OF EQUALITY CONSTRAINTS (MCON<=N)
C                      INTERNALLY REDUCED IF EQUALITY CONSTRAINTS
C                      ARE LINEARLY DEPENDENT
C        M             TREATED NUMBER OF ROWS OF MATRIX A (M>=N)
C        N             TREATED NUMBER OF COLUMNS OF MATRIX A
C     (*)IRANK         PSEUDO-RANK OF MATRIX A
C        COND          PERMITTED UPPER BOUND OF DABS(D(1)/D(IRANKC))
C                      AND OF DABS(D(IRANKC+1))/D(IRANK))
C                      (SUB-CONDITION NUMBERS OF A)
C        KRED          >=0    HOUSEHOLDER TRIANGULARIZATION
C                             (BUILD UP OF PSEUDO-INVERSE,IF IRANK<N )
C                      < 0    REDUCTION OF PSEUDO-RANK OF MATRIX A
C                             SKIPPING HOUSEHOLDER TRIANGULARIZATION
C                             BUILD-UP OF NEW PSEUDO-INVERSE
C        V(N)          REAL WORK ARRAY
C
C     OUTPUT PARAMETERS
C     -----------------
C
C        A(M,N)        OUTPUT MATRIX UPDATING PRODUCT OF HOUSEHOLDER
C                      TRANSFORMATIONS AND UPPER TRIANGULAR MATRIX
C        MCON          PSEUDO-RANK OF CONSTRAINED PART OF MATRIX A
C        IRANK         PSEUDO-RANK OF TOTAL MATRIX A
C        D(1),..,D(IRANK) DIAGONAL ELEMENTS OF UPPER TRIANGULAR MATRIX
C        PIVOT(N)      INDEX VECTOR STORING PERMUTATION OF COLUMNS
C                      DUE TO PIVOTING
C        AH(N,N)       UPDATING MATRIX FOR PART OF PSEUDO INVERSE,
C                      USED ONLY,IF IRANK < N
C
C----------------------------------------------------------------------
C
      INTEGER  IRANK, KRED, MCON, M, N, PIVOT(N)
      INTEGER  I, II, IRK1, I1, J, JD, JJ, K, K1, MH, ISUB
      DOUBLE PRECISION    A(NROW,NCOL), AH(NCOL,NCOL), D(N), V(N)
      DOUBLE PRECISION    COND, ONE , DD, DSQRT
      DOUBLE PRECISION    H, HMAX, S, T, SMALL, ZERO, EPMACH
C
      DATA  ZERO/0.D0/ , ONE/1.D0/
C
C  RELATIVE MACHINE PRECISION
      CALL ZIBCONST(EPMACH,SMALL)
C
      IF(IRANK.GT.N) IRANK=N
C
      IF(M.GT.1) GOTO 100
C  SPECIAL CASE M = 1
      PIVOT(1)=1
      D(1)=A(1,1)
      RETURN
C
100   IF  (KRED.LT.0)  GO TO  3
C
C  CONSTRAINED HOUSEHOLDER TRIANGULARIZATION
C
      DO 1 J=1,N
      PIVOT(J) = J
 1    CONTINUE
C
      JD = 1
      ISUB = 1
      MH = MCON
      IF (MH.EQ.0) MH=M
      K1 = 1
201   K = K1
      IF (K.EQ.N)  GO TO 22
      K1 = K+1
      IF (JD.EQ.0)  GO TO 211
21    DO  210  J=K,N
      S = ZERO
      DO 2101  I=K,MH
2101  S = S+A(I,J)*A(I,J)
210   D(J) = S
C
C  COLUMN PIVOTING
211   H = D(K)
      JJ = K
      DO   212  J=K1,N
      IF (D(J).LE.H)  GO TO 212
      H = D(J)
      JJ = J
212   CONTINUE
      IF (JD.EQ.1)  HMAX = H * SMALL
      JD = 0
      IF (H.GE.HMAX)  GO TO 213
      JD = 1
      GO TO 21
 213  IF (JJ.EQ.K)  GO TO 22
C
C  COLUMN INTERCHANGE
      I = PIVOT(K)
      PIVOT(K) = PIVOT(JJ)
      PIVOT(JJ) = I
      D(JJ) = D(K)
      DO  215  I=1,M
      T = A(I,K)
      A(I,K) = A(I,JJ)
215   A(I,JJ) = T
C
22    H = ZERO
      DO  221  I=K,MH
221   H = H+A(I,K)*A(I,K)
      T = DSQRT(H)
C
C  A PRIORI TEST ON PSEUDO-RANK
      IF (ISUB.GT.0) DD = T/COND
      ISUB = 0
      IF (T.GT.DD) GOTO 23
C  RANK REDUCTION
      IF (K.GT.MCON) GOTO 222
C  CONSTRAINTS ARE LINEARLY DEPENDENT
C****
      DO 2210 L=K,N
2210  A(K,L)=ZERO
C****
      MCON = K-1
      K1 = K
      MH = M
      JD = 1
      ISUB = 1
      GO TO 201
C
222   IRANK = K - 1
      IF (IRANK.EQ.0)  RETURN
      GO TO 3
C
23    S = A(K,K)
      IF (S.GT.ZERO) T = -T
      D(K) = T
      A(K,K) = S-T
      IF (K.EQ.N)  RETURN
C
      T = ONE/(H-S*T)
      DO  24  J=K1,N
      S = ZERO
      DO  241  I=K,MH
241   S = S+A(I,K)*A(I,J)
      S = S*T
      DO  242  I=K,M
242   A(I,J) = A(I,J)-A(I,K)*S
24    D(J) = D(J)-A(K,J)*A(K,J)
C
      IF (K.EQ.IRANK) GOTO 3
      IF (K.NE.MCON) GOTO 201
      MH = M
      JD = 1
      ISUB = 1
      GOTO 201
C
C  RANK-DEFICIENT PSEUDO-INVERSE
C
3     IRK1 = IRANK+1
      DO  30  J=IRK1,N
      DO  31  II=1,IRANK
      I = IRK1-II
      S = A(I,J)
      IF (II.EQ.1)  GO TO 310
      DO  3111  JJ=I1,IRANK
3111  S = S-A(I,JJ)*V(JJ)
310   I1 = I
      V(I) = S/D(I)
31    AH(I,J) = V(I)
      DO  32  I=IRK1,J
      S = ZERO
      I1 = I-1
      DO  321  JJ=1,I1
321   S = S+AH(JJ,I)*V(JJ)
      IF (I.EQ.J)  GO TO 32
      V(I) = -S/D(I)
      AH(I,J) = -V(I)
32    CONTINUE
30    D(J) = DSQRT(S+ONE)
      RETURN
C
C     **********  LAST CARD OF DECCON  **********
C
      END
      SUBROUTINE BVSOLC (A,NROW,NCOL,MCON,M,N,X,B,IRANK,D,
     @                   PIVOT,KRED,AH,V)
C
      INTEGER  IRANK, KRED, M, MCON, N, PIVOT(N)
      DOUBLE PRECISION A(NROW,NCOL), AH(NCOL,NCOL)
      DOUBLE PRECISION B(M), D(N), V(N), X(N), S, ZERO
C
C
C     BEST CONSTRAINED LINEAR LEAST SQUARES SOLUTION OF (M,N)-SYSTEM
C     FIRST MCON ROWS COMPRISE MCON EQUALITY CONSTRAINTS
C
C *********************************************************************
C
C     TO BE USED IN CONNECTION WITH SUBROUTINE DECCON (BVDECC)
C
C********** DATE OF LATEST CHANGE ************ MAY 10, 82 **********
C
C     INPUT PARAMETERS (* MARKS INOUT PARAMETERS)
C     -----------------------------------------------
C
C        A(M,N)      SEE OUTPUT OF DECCON
C        NROW        SEE OUTPUT OF DECCON
C        NCOL        SEE OUTPUT OF DECCON
C        M           SEE OUTPUT OF DECCON
C        N           SEE OUTPUT OF DECCON
C        MCON        SEE OUTPUT OF DECCON
C        IRANK       SEE OUTPUT OF DECCON
C        D(N)        SEE OUTPUT OF DECCON
C        PIVOT(N)    SEE OUTPUT OF DECCON
C        AH(N,N)     SEE OUTPUT OF DECCON
C        KRED        SEE OUTPUT OF DECCON
C      * B(M)        RIGHT-HAND SIDE OF LINEAR SYSTEM, IF (KRED.GE.0)
C                    RIGHT-HAND SIDE OF UPPER LINEAR SYSTEM,
C                                                      IF (KRED.LT.0)
C        V(N)        REAL WORK ARRAY
C
C     OUTPUT PARAMETERS
C     -----------------
C
C        X(N)        BEST LSQ-SOLUTION OF LINEAR SYSTEM
C        B(M)        RIGHT-HAND OF UPPER TRIGULAR SYSTEM
C                    (TRANSFORMED RIGHT-HAND SIDE OF LINEAR SYSTEM)
C
C
      INTEGER  I, II, I1, IRK1, J, JJ, J1, MH
C
C
      DATA  ZERO/0.D0/
C
C
      IF (IRANK.GT.0)  GO TO 110
C
C  SOLUTION FOR PSEUDO-RANK ZERO
C
      DO  1  I=1,N
1     X(I) = ZERO
      RETURN
C
110   IF (KRED.LT.0 .OR. M.EQ.1) GOTO 4
C
C  CONSTRAINED HOUSEHOLDER TRANSFORMATIONS OF RIGHT-HAND SIDE
C
3     MH = MCON
      IF (MH.EQ.0)  MH = M
      DO  31  J=1,IRANK
      S = ZERO
      DO  311  I=J,MH
311   S = S+A(I,J)*B(I)
      S = S/(D(J)*A(J,J))
      DO  312  I=J,M
312   B(I) = B(I)+A(I,J)*S
      IF (J.EQ.MCON)  MH = M
 31   CONTINUE
C
C  SOLUTION OF UPPER TRIANGULAR SYSTEM
C
4     IRK1 = IRANK+1
      DO  41  II=1,IRANK
      I = IRK1-II
      I1 = I + 1
      S = B(I)
      IF (I1.GT.IRANK)  GO TO 41
      DO  4111  JJ=I1,IRANK
4111  S = S-A(I,JJ)*V(JJ)
41    V(I) = S/D(I)
      IF (IRANK.EQ.N) GOTO 5
C
C  COMPUTATION OF THE BEST CONSTRAINED LSQ-SOLUTION
C
      DO  421  J=IRK1,N
      S = ZERO
      J1 = J-1
      DO  4211  I=1,J1
4211  S = S+AH(I,J)*V(I)
421   V(J) = -S/D(J)
      DO  422  JJ=1,N
      J = N-JJ+1
      S = ZERO
      IF (JJ.EQ.1) GOTO 4222
      DO  4221  I=J1,N
4221  S = S+AH(J,I)*V(I)
      IF (J.LE.IRANK) GOTO 4223
4222  J1=J
      V(J)=-(V(J)+S)/D(J)
      GOTO 422
4223  V(J) = V(J)-S
422   CONTINUE
C
C BACK-PERMUTATION OF SOLUTION COMPONENTS
C
5     DO  50  J=1,N
      IH=PIVOT(J)
50    X(IH) = V(J)
      RETURN
C
C     **********  LAST CARD OF SOLCON  **********
C
      END