SUBROUTINE IB01RD( JOB, N, M, L, NSMP, A, LDA, B, LDB, C, LDC, D, $ LDD, U, LDU, Y, LDY, X0, TOL, IWORK, DWORK, $ LDWORK, IWARN, INFO ) C C RELEASE 4.0, WGS COPYRIGHT 2000. C C PURPOSE C C To estimate the initial state of a linear time-invariant (LTI) C discrete-time system, given the system matrices (A,B,C,D) and C the input and output trajectories of the system. The model C structure is : C C x(k+1) = Ax(k) + Bu(k), k >= 0, C y(k) = Cx(k) + Du(k), C C where x(k) is the n-dimensional state vector (at time k), C u(k) is the m-dimensional input vector, C y(k) is the l-dimensional output vector, C and A, B, C, and D are real matrices of appropriate dimensions. C Matrix A is assumed to be in a real Schur form. C C ARGUMENTS C C Mode Parameters C C JOB CHARACTER*1 C Specifies whether or not the matrix D is zero, as follows: C = 'Z': the matrix D is zero; C = 'N': the matrix D is not zero. C C Input/Output Parameters C C N (input) INTEGER C The order of the system. N >= 0. C C M (input) INTEGER C The number of system inputs. M >= 0. C C L (input) INTEGER C The number of system outputs. L > 0. C C NSMP (input) INTEGER C The number of rows of matrices U and Y (number of C samples used, t). NSMP >= N. C C A (input) DOUBLE PRECISION array, dimension (LDA,N) C The leading N-by-N part of this array must contain the C system state matrix A in a real Schur form. C C LDA INTEGER C The leading dimension of the array A. LDA >= MAX(1,N). C C B (input) DOUBLE PRECISION array, dimension (LDB,M) C The leading N-by-M part of this array must contain the C system input matrix B (corresponding to the real Schur C form of A). C If N = 0 or M = 0, this array is not referenced. C C LDB INTEGER C The leading dimension of the array B. C LDB >= N, if N > 0 and M > 0; C LDB >= 1, if N = 0 or M = 0. C C C (input) DOUBLE PRECISION array, dimension (LDC,N) C The leading L-by-N part of this array must contain the C system output matrix C (corresponding to the real Schur C form of A). C C LDC INTEGER C The leading dimension of the array C. LDC >= L. C C D (input) DOUBLE PRECISION array, dimension (LDD,M) C The leading L-by-M part of this array must contain the C system input-output matrix. C If M = 0 or JOB = 'Z', this array is not referenced. C C LDD INTEGER C The leading dimension of the array D. C LDD >= L, if M > 0 and JOB = 'N'; C LDD >= 1, if M = 0 or JOB = 'Z'. C C U (input) DOUBLE PRECISION array, dimension (LDU,M) C If M > 0, the leading NSMP-by-M part of this array must C contain the t-by-m input-data sequence matrix U, C U = [u_1 u_2 ... u_m]. Column j of U contains the C NSMP values of the j-th input component for consecutive C time increments. C If M = 0, this array is not referenced. C C LDU INTEGER C The leading dimension of the array U. C LDU >= MAX(1,NSMP), if M > 0; C LDU >= 1, if M = 0. C C Y (input) DOUBLE PRECISION array, dimension (LDY,L) C The leading NSMP-by-L part of this array must contain the C t-by-l output-data sequence matrix Y, C Y = [y_1 y_2 ... y_l]. Column j of Y contains the C NSMP values of the j-th output component for consecutive C time increments. C C LDY INTEGER C The leading dimension of the array Y. LDY >= MAX(1,NSMP). C C X0 (output) DOUBLE PRECISION array, dimension (N) C The estimated initial state of the system, x(0). C C Tolerances C C TOL DOUBLE PRECISION C The tolerance to be used for estimating the rank of C matrices. If the user sets TOL > 0, then the given value C of TOL is used as a lower bound for the reciprocal C condition number; a matrix whose estimated condition C number is less than 1/TOL is considered to be of full C rank. If the user sets TOL <= 0, then EPS is used C instead, where EPS is the relative machine precision C (see LAPACK Library routine DLAMCH). TOL <= 1. C C Workspace C C IWORK INTEGER array, dimension (N) C C DWORK DOUBLE PRECISION array, dimension (LDWORK) C On exit, if INFO = 0, DWORK(1) returns the optimal value C of LDWORK and DWORK(2) contains the reciprocal condition C number of the triangular factor of the QR factorization of C the matrix Gamma (see METHOD). C On exit, if INFO = -22, DWORK(1) returns the minimum C value of LDWORK. C C LDWORK INTEGER C The length of the array DWORK. C LDWORK >= max( 2, min( LDW1, LDW2 ) ), where C LDW1 = t*L*(N + 1) + 2*N + max( 2*N*N, 4*N ), C LDW2 = N*(N + 1) + 2*N + C max( q*(N + 1) + 2*N*N + L*N, 4*N ), C q = N*L. C For good performance, LDWORK should be larger. C If LDWORK >= LDW1, then standard QR factorization of C the matrix Gamma (see METHOD) is used. Otherwise, the C QR factorization is computed sequentially by performing C NCYCLE cycles, each cycle (except possibly the last one) C processing s samples, where s is chosen by equating C LDWORK to LDW2, for q replaced by s*L. C The computational effort may increase and the accuracy may C decrease with the decrease of s. Recommended value is C LDRWRK = LDW1, assuming a large enough cache size, to C also accommodate A, B, C, D, U, and Y. C C Warning Indicator C C IWARN INTEGER C = 0: no warning; C = 4: the least squares problem to be solved has a C rank-deficient coefficient matrix. C C Error Indicator C C INFO INTEGER C = 0: successful exit; C < 0: if INFO = -i, the i-th argument had an illegal C value; C = 2: the singular value decomposition (SVD) algorithm did C not converge. C C METHOD C C An extension and refinement of the method in [1] is used. C Specifically, the output y0(k) of the system for zero initial C state is computed for k = 0, 1, ..., t-1 using the given model. C Then the following least squares problem is solved for x(0) C C ( C ) ( y(0) - y0(0) ) C ( C*A ) ( y(1) - y0(1) ) C Gamma * x(0) = ( : ) * x(0) = ( : ). C ( : ) ( : ) C ( C*A^(t-1) ) ( y(t-1) - y0(t-1) ) C C The coefficient matrix Gamma is evaluated using powers of A with C exponents 2^k. The QR decomposition of this matrix is computed. C If its triangular factor R is too ill conditioned, then singular C value decomposition of R is used. C C If the coefficient matrix cannot be stored in the workspace (i.e., C LDWORK < LDW1), the QR decomposition is computed sequentially. C C REFERENCES C C [1] Verhaegen M., and Varga, A. C Some Experience with the MOESP Class of Subspace Model C Identification Methods in Identifying the BO105 Helicopter. C Report TR R165-94, DLR Oberpfaffenhofen, 1994. C C NUMERICAL ASPECTS C C The implemented method is numerically stable. C C CONTRIBUTOR C C V. Sima, Research Institute for Informatics, Bucharest, Apr. 2000. C C REVISIONS C C - C C KEYWORDS C C Identification methods; least squares solutions; multivariable C systems; QR decomposition; singular value decomposition. C C ****************************************************************** C C .. Parameters .. DOUBLE PRECISION ZERO, ONE, TWO, THREE PARAMETER ( ZERO = 0.0D0, ONE = 1.0D0, TWO = 2.0D0, $ THREE = 3.0D0 ) C IBLOCK is a threshold value for switching to a block algorithm C for U (to avoid row by row passing through U). INTEGER IBLOCK PARAMETER ( IBLOCK = 16384 ) C .. Scalar Arguments .. DOUBLE PRECISION TOL INTEGER INFO, IWARN, L, LDA, LDB, LDC, LDD, LDU, $ LDWORK, LDY, M, N, NSMP CHARACTER JOB C .. Array Arguments .. DOUBLE PRECISION A(LDA, *), B(LDB, *), C(LDC, *), D(LDD, *), $ DWORK(*), U(LDU, *), X0(*), Y(LDY, *) INTEGER IWORK(*) C .. Local Scalars .. DOUBLE PRECISION RCOND, TOLL INTEGER I2, IA, IAS, IC, ICYCLE, IE, IERR, IEXPON, $ IG, INIGAM, INIH, INIR, INIT, IQ, IREM, IRHS, $ ISIZE, ISV, ITAU, IU, IUPNT, IUT, IUTRAN, IX, $ IXINIT, IY, IYPNT, J, JWORK, K, LDDW, LDR, $ LDW1, LDW2, MAXWRK, MINSMP, MINWLS, MINWRK, NC, $ NCP1, NCYCLE, NN, NOBS, NRBL, NROW, NSMPL, RANK LOGICAL BLOCK, FIRST, NCYC, POWER2, SWITCH, WITHD C .. Local Arrays .. DOUBLE PRECISION DUM(1) C .. External Functions .. LOGICAL LSAME INTEGER ILAENV DOUBLE PRECISION DLAMCH EXTERNAL DLAMCH, ILAENV, LSAME C .. External Subroutines .. EXTERNAL DAXPY, DCOPY, DGELSS, DGEMV, DGEQRF, DLACPY, $ DLASET, DORMQR, DTRCON, DTRMM, DTRMV, DTRSV, $ MA02AD, MB01TD, MB04OD, XERBLA C .. Intrinsic Functions .. INTRINSIC DBLE, INT, LOG, MAX, MIN, MOD C .. Executable Statements .. C C Check the input parameters. C WITHD = LSAME( JOB, 'N' ) IWARN = 0 INFO = 0 NN = N*N MINSMP = N C IF( .NOT.( LSAME( JOB, 'Z' ) .OR. WITHD ) ) THEN INFO = -1 ELSE IF( N.LT.0 ) THEN INFO = -2 ELSE IF( M.LT.0 ) THEN INFO = -3 ELSE IF( L.LE.0 ) THEN INFO = -4 ELSE IF( NSMP.LT.MINSMP ) THEN INFO = -5 ELSE IF( LDA.LT.MAX( 1, N ) ) THEN INFO = -7 ELSE IF( LDB.LT.1 .OR. ( LDB.LT.N .AND. M.GT.0 ) ) THEN INFO = -9 ELSE IF( LDC.LT.L ) THEN INFO = -11 ELSE IF( LDD.LT.1 .OR. ( WITHD .AND. LDD.LT.L .AND. M.GT.0 ) ) $ THEN INFO = -13 ELSE IF( LDU.LT.1 .OR. ( M.GT.0 .AND. LDU.LT.NSMP ) ) THEN INFO = -15 ELSE IF( LDY.LT.MAX( 1, NSMP ) ) THEN INFO = -17 ELSE IF( TOL.GT.ONE ) THEN INFO = -19 END IF C C Compute workspace. C (Note: Comments in the code beginning "Workspace:" describe the C minimal amount of workspace needed at that point in the code, C as well as the preferred amount for good performance. C NB refers to the optimal block size for the immediately C following subroutine, as returned by ILAENV.) C NSMPL = NSMP*L IQ = MINSMP*L NCP1 = N + 1 ISIZE = NSMPL*NCP1 IC = 2*NN MINWLS = MINSMP*NCP1 ITAU = IC + L*N LDW1 = ISIZE + 2*N + MAX( IC, 4*N ) LDW2 = MINWLS + 2*N + MAX( IQ*NCP1 + ITAU, 4*N ) MINWRK = MAX( MIN( LDW1, LDW2 ), 2 ) IF ( INFO.EQ.0 .AND. LDWORK.GE.MINWRK ) THEN MAXWRK = ISIZE + 2*N + MAX( N*ILAENV( 1, 'DGEQRF', ' ', NSMPL, $ N, -1, -1 ), $ ILAENV( 1, 'DORMQR', 'LT', NSMPL, $ 1, N, -1 ) ) MAXWRK = MAX( MAXWRK, MINWRK ) END IF C IF ( INFO.EQ.0 .AND. LDWORK.LT.MINWRK ) THEN INFO = -22 DWORK(1) = MINWRK END IF C C Return if there are illegal arguments. C IF( INFO.NE.0 ) THEN CALL XERBLA( 'IB01RD', -INFO ) RETURN END IF C C Quick return if possible. C IF ( N.EQ.0 ) THEN DWORK(1) = TWO DWORK(2) = ONE RETURN END IF C C Set up the least squares problem, either directly, if enough C workspace, or sequentially, otherwise. C IYPNT = 1 IUPNT = 1 INIR = 1 IF ( LDWORK.GE.LDW1 ) THEN C C Enough workspace for solving the problem directly. C NCYCLE = 1 NOBS = NSMP LDDW = NSMPL INIGAM = 1 ELSE C C NCYCLE > 1 cycles are needed for solving the problem C sequentially, taking NOBS samples in each cycle (or the C remaining samples in the last cycle). C JWORK = LDWORK - MINWLS - 2*N - ITAU LDDW = JWORK/NCP1 NOBS = LDDW/L LDDW = L*NOBS NCYCLE = NSMP/NOBS IF ( MOD( NSMP, NOBS ).NE.0 ) $ NCYCLE = NCYCLE + 1 INIH = INIR + NN INIGAM = INIH + N END IF C NCYC = NCYCLE.GT.1 IRHS = INIGAM + LDDW*N IXINIT = IRHS + LDDW IC = IXINIT + N IF( NCYC ) THEN IA = IC + L*N LDR = N IE = INIGAM ELSE INIH = IRHS IA = IC LDR = LDDW IE = IXINIT END IF IUTRAN = IA IAS = IA + NN ITAU = IA DUM(1) = ZERO C C Set block parameters for passing through the array U. C BLOCK = M.GT.1 .AND. NSMP*M.GE.IBLOCK IF ( BLOCK ) THEN NRBL = ( LDWORK - IUTRAN + 1 )/M NC = NOBS/NRBL IF ( MOD( NOBS, NRBL ).NE.0 ) $ NC = NC + 1 INIT = ( NC - 1 )*NRBL BLOCK = BLOCK .AND. NRBL.GT.1 END IF C C Perform direct of sequential compression of the matrix Gamma. C DO 150 ICYCLE = 1, NCYCLE FIRST = ICYCLE.EQ.1 IF ( .NOT.FIRST ) THEN IF ( ICYCLE.EQ.NCYCLE ) THEN NOBS = NSMP - ( NCYCLE - 1 )*NOBS LDDW = L*NOBS IF ( BLOCK ) THEN NC = NOBS/NRBL IF ( MOD( NOBS, NRBL ).NE.0 ) $ NC = NC + 1 INIT = ( NC - 1 )*NRBL END IF END IF END IF C C Compute the extended observability matrix Gamma. C Workspace: need s*L*(N + 1) + 2*N*N + 2*N + a + w, C where s = NOBS, C a = 0, w = 0, if NCYCLE = 1, C a = L*N, w = N*(N + 1), if NCYCLE > 1; C prefer as above, with s = t, a = w = 0. C JWORK = IAS + NN IEXPON = INT( LOG( DBLE( NOBS ) )/LOG( TWO ) ) IREM = L*( NOBS - 2**IEXPON ) POWER2 = IREM.EQ.0 IF ( .NOT.POWER2 ) $ IEXPON = IEXPON + 1 C IF ( FIRST ) THEN CALL DLACPY( 'Full', L, N, C, LDC, DWORK(INIGAM), LDDW ) ELSE CALL DLACPY( 'Full', L, N, DWORK(IC), L, DWORK(INIGAM), $ LDDW ) END IF C p C Use powers of the matrix A: A , p = 2**(J-1). C CALL DLACPY( 'Upper', N, N, A, LDA, DWORK(IA), N ) CALL DCOPY( N-1, A(2,1), LDA+1, DWORK(IA+1), N+1 ) I2 = L NROW = 0 C DO 20 J = 1, IEXPON IG = INIGAM IF ( J.LT.IEXPON .OR. POWER2 ) THEN NROW = I2 ELSE NROW = IREM END IF C CALL DLACPY( 'Full', NROW, N, DWORK(IG), LDDW, DWORK(IG+I2), $ LDDW ) CALL DTRMM( 'Right', 'Upper', 'No Transpose', 'Non Unit', $ NROW, N, ONE, DWORK(IA), N, DWORK(IG+I2), $ LDDW ) C p C Compute the contribution of the subdiagonal of A to the C product. C DO 10 IX = 1, N - 1 CALL DAXPY( NROW, DWORK(IA+(IX-1)*N+IX), DWORK(IG+LDDW), $ 1, DWORK(IG+I2), 1 ) IG = IG + LDDW 10 CONTINUE C IF ( J.LT.IEXPON ) THEN CALL DLACPY( 'Upper', N, N, DWORK(IA), N, DWORK(IAS), N ) CALL DCOPY( N-1, DWORK(IA+1), N+1, DWORK(IAS+1), N+1 ) CALL MB01TD( N, DWORK(IAS), N, DWORK(IA), N, $ DWORK(JWORK), IERR ) I2 = I2*2 END IF 20 CONTINUE C IF ( NCYC ) THEN IG = INIGAM + I2 + NROW - L CALL DLACPY( 'Full', L, N, DWORK(IG), LDDW, DWORK(IC), L ) CALL DTRMM( 'Right', 'Upper', 'No Transpose', 'Non Unit', L, $ N, ONE, A, LDA, DWORK(IC), L ) C C Compute the contribution of the subdiagonal of A to the C product. C DO 30 IX = 1, N - 1 CALL DAXPY( L, A(IX+1,IX), DWORK(IG+LDDW), 1, $ DWORK(IC+(IX-1)*L), 1 ) IG = IG + LDDW 30 CONTINUE C END IF C C Setup (part of) the right hand side of the least squares C problem starting from DWORK(IRHS); use the estimated output C trajectory for zero initial state, or for the saved final state C value of the previous cycle. C A specialization of SLICOT Library routine TF01ND is used. C For large input sets (NSMP*M >= IBLOCK), chunks of U are C transposed, to reduce the number of row-wise passes. C Workspace: need s*L*(N + 1) + N + w; C prefer as above, with s = t, w = 0. C IF ( FIRST ) $ CALL DCOPY( N, DUM, 0, DWORK(IXINIT), 1 ) CALL DCOPY( N, DWORK(IXINIT), 1, X0, 1 ) IY = IRHS C DO 40 J = 1, L CALL DCOPY( NOBS, Y(IYPNT,J), 1, DWORK(IY), L ) IY = IY + 1 40 CONTINUE C IY = IRHS IU = IUPNT IF ( M.GT.0 ) THEN IF ( WITHD ) THEN C IF ( BLOCK ) THEN SWITCH = .TRUE. NROW = NRBL C DO 60 K = 1, NOBS IF ( MOD( K-1, NROW ).EQ.0 .AND. SWITCH ) THEN IUT = IUTRAN IF ( K.GT.INIT ) THEN NROW = NOBS - INIT SWITCH = .FALSE. END IF CALL MA02AD( 'Full', NROW, M, U(IU,1), LDU, $ DWORK(IUT), M ) IU = IU + NROW END IF CALL DGEMV( 'No transpose', L, N, -ONE, C, LDC, X0, $ 1, ONE, DWORK(IY), 1 ) CALL DGEMV( 'No transpose', L, M, -ONE, D, LDD, $ DWORK(IUT), 1, ONE, DWORK(IY), 1 ) CALL DTRMV( 'Upper', 'No transpose', 'Non-unit', N, $ A, LDA, X0, 1 ) C DO 50 IX = 2, N X0(IX) = X0(IX) + A(IX,IX-1)*DWORK(IXINIT+IX-2) 50 CONTINUE C CALL DGEMV( 'No transpose', N, M, ONE, B, LDB, $ DWORK(IUT), 1, ONE, X0, 1 ) CALL DCOPY( N, X0, 1, DWORK(IXINIT), 1 ) IY = IY + L IUT = IUT + M 60 CONTINUE C ELSE C DO 80 K = 1, NOBS CALL DGEMV( 'No transpose', L, N, -ONE, C, LDC, X0, $ 1, ONE, DWORK(IY), 1 ) CALL DGEMV( 'No transpose', L, M, -ONE, D, LDD, $ U(IU,1), LDU, ONE, DWORK(IY), 1 ) CALL DTRMV( 'Upper', 'No transpose', 'Non-unit', N, $ A, LDA, X0, 1 ) C DO 70 IX = 2, N X0(IX) = X0(IX) + A(IX,IX-1)*DWORK(IXINIT+IX-2) 70 CONTINUE C CALL DGEMV( 'No transpose', N, M, ONE, B, LDB, $ U(IU,1), LDU, ONE, X0, 1 ) CALL DCOPY( N, X0, 1, DWORK(IXINIT), 1 ) IY = IY + L IU = IU + 1 80 CONTINUE C END IF C ELSE C IF ( BLOCK ) THEN SWITCH = .TRUE. NROW = NRBL C DO 100 K = 1, NOBS IF ( MOD( K-1, NROW ).EQ.0 .AND. SWITCH ) THEN IUT = IUTRAN IF ( K.GT.INIT ) THEN NROW = NOBS - INIT SWITCH = .FALSE. END IF CALL MA02AD( 'Full', NROW, M, U(IU,1), LDU, $ DWORK(IUT), M ) IU = IU + NROW END IF CALL DGEMV( 'No transpose', L, N, -ONE, C, LDC, X0, $ 1, ONE, DWORK(IY), 1 ) CALL DTRMV( 'Upper', 'No transpose', 'Non-unit', N, $ A, LDA, X0, 1 ) C DO 90 IX = 2, N X0(IX) = X0(IX) + A(IX,IX-1)*DWORK(IXINIT+IX-2) 90 CONTINUE C CALL DGEMV( 'No transpose', N, M, ONE, B, LDB, $ DWORK(IUT), 1, ONE, X0, 1 ) CALL DCOPY( N, X0, 1, DWORK(IXINIT), 1 ) IY = IY + L IUT = IUT + M 100 CONTINUE C ELSE C DO 120 K = 1, NOBS CALL DGEMV( 'No transpose', L, N, -ONE, C, LDC, X0, $ 1, ONE, DWORK(IY), 1 ) CALL DTRMV( 'Upper', 'No transpose', 'Non-unit', N, $ A, LDA, X0, 1 ) C DO 110 IX = 2, N X0(IX) = X0(IX) + A(IX,IX-1)*DWORK(IXINIT+IX-2) 110 CONTINUE C CALL DGEMV( 'No transpose', N, M, ONE, B, LDB, $ U(IU,1), LDU, ONE, X0, 1 ) CALL DCOPY( N, X0, 1, DWORK(IXINIT), 1 ) IY = IY + L IU = IU + 1 120 CONTINUE C END IF C END IF C ELSE C DO 140 K = 1, NOBS CALL DGEMV( 'No transpose', L, N, -ONE, C, LDC, X0, 1, $ ONE, DWORK(IY), 1 ) CALL DTRMV( 'Upper', 'No transpose', 'Non-unit', N, A, $ LDA, X0, 1 ) C DO 130 IX = 2, N X0(IX) = X0(IX) + A(IX,IX-1)*DWORK(IXINIT+IX-2) 130 CONTINUE C CALL DCOPY( N, X0, 1, DWORK(IXINIT), 1 ) IY = IY + L 140 CONTINUE C END IF C C Compress the data using (sequential) QR factorization. C Workspace: need v + 2*N; C where v = s*L*(N + 1) + N + a + w. C JWORK = ITAU + N IF ( FIRST ) THEN C C Compress the first data segment of Gamma. C Workspace: need v + 2*N, C prefer v + N + N*NB. C CALL DGEQRF( LDDW, N, DWORK(INIGAM), LDDW, DWORK(ITAU), $ DWORK(JWORK), LDWORK-JWORK+1, IERR ) C C Apply the transformation to the right hand side part. C Workspace: need v + N + 1, C prefer v + N + NB. C CALL DORMQR( 'Left', 'Transpose', LDDW, 1, N, DWORK(INIGAM), $ LDDW, DWORK(ITAU), DWORK(IRHS), LDDW, $ DWORK(JWORK), LDWORK-JWORK+1, IERR ) C IF ( NCYC ) THEN C C Save the triangular factor of Gamma and the C corresponding right hand side. C CALL DLACPY( 'Upper', N, NCP1, DWORK(INIGAM), LDDW, $ DWORK(INIR), LDR ) END IF ELSE C C Compress the current (but not the first) data segment of C Gamma. C Workspace: need v + N - 1. C CALL MB04OD( 'Full', N, 1, LDDW, DWORK(INIR), LDR, $ DWORK(INIGAM), LDDW, DWORK(INIH), LDR, $ DWORK(IRHS), LDDW, DWORK(ITAU), DWORK(JWORK) ) END IF C IUPNT = IUPNT + NOBS IYPNT = IYPNT + NOBS 150 CONTINUE C C Estimate the reciprocal condition number of the triangular factor C of the QR decomposition. C Workspace: need u + 3*N, where C u = t*L*(N + 1), if NCYCLE = 1; C u = w, if NCYCLE > 1. C CALL DTRCON( '1-norm', 'Upper', 'No Transpose', N, DWORK(INIR), $ LDR, RCOND, DWORK(IE), IWORK, IERR ) C TOLL = TOL IF ( TOLL.LE.ZERO ) $ TOLL = DLAMCH( 'Precision' ) IF ( RCOND.LE.TOLL**( TWO/THREE ) ) THEN IWARN = 4 C C The least squares problem is ill-conditioned. C Use SVD to solve it. C Workspace: need u + 6*N; C prefer larger. C CALL DLASET( 'Lower', N-1, N-1, ZERO, ZERO, DWORK(INIR+1), $ LDR ) ISV = IE JWORK = ISV + N CALL DGELSS( N, N, 1, DWORK(INIR), LDR, DWORK(INIH), LDR, $ DWORK(ISV), TOLL, RANK, DWORK(JWORK), $ LDWORK-JWORK+1, IERR ) IF ( IERR.GT.0 ) THEN C C Return if SVD algorithm did not converge. C INFO = 2 RETURN END IF MAXWRK = MAX( MAXWRK, INT( DWORK(JWORK) ) - JWORK + 1 ) ELSE C C Find the least squares solution using QR decomposition only. C CALL DTRSV( 'Upper', 'No Transpose', 'Non Unit', N, $ DWORK(INIR), LDR, DWORK(INIH), 1 ) END IF C C Return the estimated initial state of the system x0. C CALL DCOPY( N, DWORK(INIH), 1, X0, 1 ) C DWORK(1) = MAXWRK DWORK(2) = RCOND C RETURN C C *** End of IB01RD *** END