SUBROUTINE PEDECC(A,NROW,NCOL,MCON,M,N,IRANK,COND,D,PIVOT, *KRED,AH,V) C* Begin Prologue DECCON INTEGER IRANK,MCON INTEGER M,N,NROW,NCOL,KRED INTEGER PIVOT(NCOL) DOUBLE PRECISION COND DOUBLE PRECISION A(NROW,NCOL),AH(NCOL,NCOL) DOUBLE PRECISION D(NCOL),V(NCOL) C ------------------------------------------------------------ C C* Title C C* Deccon - Constrained Least Squares QR-Decomposition C C* Written by P. Deuflhard, L.Weimann C* Purpose Solution of least squares problems, optionally C with equality constraints. C* Method Constrained Least Squares QR-Decomposition C (see references below) C* Category D9b1. - Singular, overdetermined or C underdetermined systems of linear C equations, generalized inverses. C Constrained Least Squares solution C* Keywords Linear Least Square Problems, constrained, C QR-decomposition, pseudo inverse. C* Version 1.0 C* Revision August 1987 C* Latest Change January 1991 C* Library CodeLib C* Code Fortran 77, Double Precision C* Environment Standard Fortran 77 environment on PC's, C workstations and hosts. C* Copyright (c) Konrad-Zuse-Zentrum fuer C Informationstechnik Berlin (ZIB) C Takustrasse 7, D-14195 Berlin-Dahlem C phone : + 49/30/84185-0 C fax : + 49/30/84185-125 C* Contact Bodo Erdmann C ZIB, Division Scientific Computing, C Department Numerical Analysis and Modelling C phone : + 49/30/84185-185 C fax : + 49/30/84185-107 C e-mail: erdmann@zib.de C C* References: C =========== C C /1/ P.Deuflhard, V.Apostolescu: C An underrelaxed Gauss-Newton method for equality C constrained nonlinear least squares problems. C Lecture Notes Control Inform. Sci. vol. 7, p. C 22-32 (1978) C /2/ P.Deuflhard, W.Sautter: C On rank-deficient pseudoinverses. C J. Lin. Alg. Appl. vol. 29, p. 91-111 (1980) C C* Related Programs: SOLCON (PESOLC) C C --------------------------------------------------------------- C C* Licence C You may use or modify this code for your own non commercial C purposes for an unlimited time. C In any case you should not deliver this code without a special C permission of ZIB. C In case you intend to use the code commercially, we oblige you C to sign an according licence agreement with ZIB. C C* Warranty C This code has been tested up to a certain level. Defects and C weaknesses, which may be included in the code, do not establish C any warranties by ZIB. ZIB does not take over any liabilities C which may follow from aquisition or application of this code. C C* Software status C This code is under partial care of ZIB and belongs to ZIB software C class 2. C C ------------------------------------------------------------ C C* Summary: C ======== C Constrained QR-decomposition of (M,N)-system with C computation of pseudoinverse in case of rank-defeciency . C First MCON rows belong to equality constraints. C C ------------------------------------------------------------ C C* Parameters: C =========== C C* Input parameters (* marks inout parameters) 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 C AH C * MCON number of equality constraints C MCON.LE.N C internally reduced if equality C constraints C are linearly dependent C M treated number of rows of matrix A 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 C (IRANK))and of DABS(D(IRANK+1)/D( C IRANK)) C ( sub - condition numbers of A ) C KRED Type of operation C >= 0 Householder triangularization C ( build up of pseudo - C inverse,if IRANK.LT.N ) C < 0 reduction of pseudo-rank of C matrix A , skipping Householder C triangularization, build-up of C new pseudo-inverse C V(N) real work array C C* Output parameters C ================= C C A(M,N) Output matrix updating product of C Householder transformations and C upper triangular matrix C MCON Pseudo-rank of constrained part of C matrix A C IRANK Pseudo-rank of total matrix A C D(IRANK) Diagonal elements of upper C triangular matrix C PIVOT(N) Index vector storing permutation of C columns due to pivoting C COND Sub-condition number of A C (in case of rank reduction: C sub-condition number which led to C rank reduction) C AH(N,N) Updating matrix for part of pseudo C inverse C C* Subroutines called: none C C* Machine constants used C ====================== C C EPMACH = relative machine precision C DOUBLE PRECISION EPMACH C C ____________________________________________________________ C* End Prologue INTRINSIC DABS,DSQRT DOUBLE PRECISION EPMACH,SMALL DOUBLE PRECISION ZERO PARAMETER (ZERO=0.0D0) DOUBLE PRECISION ONE PARAMETER (ONE=1.0D0) INTEGER L1 DOUBLE PRECISION S1 INTEGER I,II,IRK1,ISUB,I1,J,JD,JJ,K,K1,LEVEL,MH DOUBLE PRECISION DD,H,HMAX,S,SH,SMALLS,T C* Begin C ______________________________________________________________ C 1 Initialization CALL ZIBCONST(EPMACH,SMALL) SMALLS = DSQRT(EPMACH*10.0D0) IF(IRANK.GT.N) IRANK = N IF(IRANK.GT.M) IRANK = M C ______________________________________________________________ C 1.1 Special case M=1 and N=1 IF(M.EQ.1.AND.N.EQ.1)THEN PIVOT(1)=1 D(1)=A(1,1) COND = ONE RETURN ENDIF IF(KRED.GE.0)THEN C ____________________________________________________________ C 1.1 Initialize pivot-array DO 11 J=1,N PIVOT(J)=J 11 CONTINUE C ____________________________________________________________ C 2. Constrained Householder triangularization JD = 1 ISUB = 1 MH = MCON IF(MH.EQ.0) MH = M K1 = 1 LEVEL = 1 C DO (Until) 2 CONTINUE K = K1 IF(K.NE.N)THEN K1 = K+1 C DO (Until) 20 CONTINUE IF(JD.NE.0)THEN DO 201 J=K,N S = ZERO DO 2011 L1=K,MH S = S+A(L1,J)**2 2011 CONTINUE D(J)=S 201 CONTINUE ENDIF C ______________________________________________________ C 2.1 Column pivoting S1 = D(K) JJ = K DO 21 L1=K,N IF(D(L1).GT.S1) THEN S1=D(L1) JJ = L1 ENDIF 21 CONTINUE H = D(JJ) IF(JD.EQ.1) HMAX = H*SMALLS JD = 0 IF(H.LT.HMAX) JD = 1 IF(.NOT.(H.GE.HMAX)) GOTO 20 C UNTIL ( expression - negated above) IF(JJ.NE.K)THEN C ______________________________________________________ C 2.2 Column interchange I = PIVOT(K) PIVOT(K)=PIVOT(JJ) PIVOT(JJ)=I D(JJ)=D(K) DO 221 L1=1,M S1=A(L1,JJ) A(L1,JJ)=A(L1,K) A(L1,K)=S1 221 CONTINUE ENDIF ENDIF H = ZERO DO 222 L1=K,MH H = H+A(L1,K)**2 222 CONTINUE T = DSQRT(H) C __________________________________________________________ C 2.3.0 A-priori test on pseudo-rank IF(ISUB.GT.0) DD = T/COND ISUB = 0 IF(T.GT.DD.OR.K.LE.MCON)THEN IF(T.LE.DD)THEN C ______________________________________________________ C 2.3.1 Rank reduction MCON = K-1 K1 = K MH = M JD = 1 ISUB = 1 ELSE C ______________________________________________________ C 2.4 Householder step S = A(K,K) IF(S.GT.ZERO) T =-T D(K)=T A(K,K)=S-T IF(K.NE.N)THEN T = ONE/(H-S*T) DO 24 J=K1,N S = ZERO DO 241 L1=K,MH S = S+A(L1,K)*A(L1,J) 241 CONTINUE S = S*T S1 =-S DO 242 L1=K,M A(L1,J) = A(L1,J)+A(L1,K)*S1 242 CONTINUE D(J)=D(J)-A(K,J)**2 24 CONTINUE IF(K.NE.IRANK)THEN IF(K.EQ.MCON)THEN MH = M JD = 1 ISUB = 1 ENDIF ELSE LEVEL = 3 ENDIF ELSE LEVEL = 4 ENDIF ENDIF ELSE IRANK = K-1 IF(IRANK.EQ.0)THEN LEVEL = 4 ELSE LEVEL = 3 ENDIF ENDIF IF(.NOT.(LEVEL.NE.1)) GOTO 2 C UNTIL ( expression - negated above) ELSE LEVEL = 3 ENDIF IF(LEVEL.EQ.3)THEN C ______________________________________________________________ C 3 Rank-deficient pseudo-inverse IRK1 = IRANK+1 DO 3 J=IRK1,N DO 31 II=1,IRANK I = IRK1-II S = A(I,J) IF(II.NE.1)THEN SH = ZERO DO 311 L1=I1,IRANK SH=SH+A(I,L1)*V(L1) 311 CONTINUE S = S-SH ENDIF I1 = I V(I)=S/D(I) AH(I,J)=V(I) 31 CONTINUE DO 32 I=IRK1,J S = ZERO DO 321 L1=1,I-1 S = S+AH(L1,I)*V(L1) 321 CONTINUE IF(I.NE.J)THEN V(I)=-S/D(I) AH(I,J)=-V(I) ENDIF 32 CONTINUE D(J)=DSQRT(S+ONE) 3 CONTINUE ENDIF C ______________________________________________________________ C 9 Exit IF(K.EQ.IRANK)THEN T = D(IRANK) IF(T.NE.ZERO) COND = DABS(D(1)/T) ENDIF RETURN END SUBROUTINE PESOLC(A,NROW,NCOL,MCON,M,N,X,B,IRANK,D,PIVOT, *KRED,AH,V) C* Begin Prologue SOLCON DOUBLE PRECISION A(NROW,NCOL),AH(NCOL,NCOL) DOUBLE PRECISION X(NCOL),B(NROW),D(NCOL),V(NCOL) INTEGER NROW,NCOL,MCON,M,N,IRANK,KRED INTEGER PIVOT(NCOL) C ____________________________________________________________ C C* Summary C ======= C C Best constrained linear least squares solution of (M,N)- C system . First MCON rows comprise MCON equality constraints. C To be used in connection with subroutine DECCON (PEDECC) C References: See DECCON C Related Programs: DECCON C C* Parameters: C =========== 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 C 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 C system ) C C ____________________________________________________________ C* End Prologue DOUBLE PRECISION ZERO PARAMETER (ZERO=0.0D0) INTEGER L1,L2 DOUBLE PRECISION S1 INTEGER I,II,I1,IRK1,J,JJ,J1,MH DOUBLE PRECISION S,SH C* Begin C ____________________________________________________________ C 1 Solution for pseudo-rank zero IF(IRANK.EQ.0)THEN S1 = ZERO DO 1 L1=1,N X(L1)=S1 1 CONTINUE RETURN ENDIF IF(KRED.GE.0.AND.(M.NE.1.OR.N.NE.1))THEN C __________________________________________________________ C 2 Constrained householder transformations of right-hand side MH = MCON IF(MH.EQ.0) MH = M DO 2 J=1,IRANK S = ZERO DO 21 L1=J,MH S = S+A(L1,J)*B(L1) 21 CONTINUE S = S/(D(J)*A(J,J)) S1 = S DO 22 L1=J,M B(L1)=B(L1)+A(L1,J)*S1 22 CONTINUE IF(J.EQ.MCON) MH = M 2 CONTINUE ENDIF C ____________________________________________________________ C 3 Solution of upper triangular system IRK1 = IRANK+1 DO 31 II=1,IRANK I = IRK1-II I1 = I+1 S = B(I) IF(I1.LE.IRANK)THEN SH = ZERO DO 311 L1=I1,IRANK SH=SH+A(I,L1)*V(L1) 311 CONTINUE S = S-SH ENDIF V(I)=S/D(I) 31 CONTINUE IF(IRK1.LE.N)THEN C __________________________________________________________ C 3.2 Computation of the best constrained least squares- C solution DO 321 J=IRK1,N S = ZERO DO 3211 L1=1,J-1 S = S+AH(L1,J)*V(L1) 3211 CONTINUE V(J)=-S/D(J) 321 CONTINUE DO 322 JJ=1,N J = N-JJ+1 S = ZERO IF(JJ.NE.1)THEN DO 3221 L1=J1,N S=S+AH(J,L1)*V(L1) 3221 CONTINUE ENDIF IF(JJ.NE.1.AND.J.LE.IRANK)THEN V(J)=V(J)-S ELSE J1 = J V(J)=-(S+V(J))/D(J) ENDIF 322 CONTINUE ENDIF C ____________________________________________________________ C 4 Back-permutation of solution components DO 4 L1=1,N L2 = PIVOT(L1) X(L2) = V(L1) 4 CONTINUE RETURN END