1 SUBROUTINE IB01PD( METH, JOB, JOBCV, NOBR, N, M, L, NSMPL, R, 2 $ LDR, A, LDA, C, LDC, B, LDB, D, LDD, Q, LDQ, 3 $ RY, LDRY, S, LDS, O, LDO, TOL, IWORK, DWORK, 4 $ LDWORK, IWARN, INFO ) 5C 6C RELEASE 4.0, WGS COPYRIGHT 2000. 7C 8C PURPOSE 9C 10C To estimate the matrices A, C, B, and D of a linear time-invariant 11C (LTI) state space model, using the singular value decomposition 12C information provided by other routines. Optionally, the system and 13C noise covariance matrices, needed for the Kalman gain, are also 14C determined. 15C 16C ARGUMENTS 17C 18C Mode Parameters 19C 20C METH CHARACTER*1 21C Specifies the subspace identification method to be used, 22C as follows: 23C = 'M': MOESP algorithm with past inputs and outputs; 24C = 'N': N4SID algorithm. 25C 26C JOB CHARACTER*1 27C Specifies which matrices should be computed, as follows: 28C = 'A': compute all system matrices, A, B, C, and D; 29C = 'C': compute the matrices A and C only; 30C = 'B': compute the matrix B only; 31C = 'D': compute the matrices B and D only. 32C 33C JOBCV CHARACTER*1 34C Specifies whether or not the covariance matrices are to 35C be computed, as follows: 36C = 'C': the covariance matrices should be computed; 37C = 'N': the covariance matrices should not be computed. 38C 39C Input/Output Parameters 40C 41C NOBR (input) INTEGER 42C The number of block rows, s, in the input and output 43C Hankel matrices processed by other routines. NOBR > 1. 44C 45C N (input) INTEGER 46C The order of the system. NOBR > N > 0. 47C 48C M (input) INTEGER 49C The number of system inputs. M >= 0. 50C 51C L (input) INTEGER 52C The number of system outputs. L > 0. 53C 54C NSMPL (input) INTEGER 55C If JOBCV = 'C', the total number of samples used for 56C calculating the covariance matrices. 57C NSMPL >= 2*(M+L)*NOBR. 58C This parameter is not meaningful if JOBCV = 'N'. 59C 60C R (input/workspace) DOUBLE PRECISION array, dimension 61C ( LDR,2*(M+L)*NOBR ) 62C On entry, the leading 2*(M+L)*NOBR-by-2*(M+L)*NOBR part 63C of this array must contain the relevant data for the MOESP 64C or N4SID algorithms, as constructed by SLICOT Library 65C routines IB01AD or IB01ND. Let R_ij, i,j = 1:4, be the 66C ij submatrix of R (denoted S in IB01AD and IB01ND), 67C partitioned by M*NOBR, L*NOBR, M*NOBR, and L*NOBR 68C rows and columns. The submatrix R_22 contains the matrix 69C of left singular vectors used. Also needed, for 70C METH = 'N' or JOBCV = 'C', are the submatrices R_11, 71C R_14 : R_44, and, for METH = 'M' and JOB <> 'C', the 72C submatrices R_31 and R_12, containing the processed 73C matrices R_1c and R_2c, respectively, as returned by 74C SLICOT Library routines IB01AD or IB01ND. 75C Moreover, if METH = 'N' and JOB = 'A' or 'C', the 76C block-row R_41 : R_43 must contain the transpose of the 77C block-column R_14 : R_34 as returned by SLICOT Library 78C routines IB01AD or IB01ND. 79C The remaining part of R is used as workspace. 80C On exit, part of this array is overwritten. Specifically, 81C if METH = 'M', R_22 and R_31 are overwritten if 82C JOB = 'B' or 'D', and R_12, R_22, R_14 : R_34, 83C and possibly R_11 are overwritten if JOBCV = 'C'; 84C if METH = 'N', all needed submatrices are overwritten. 85C 86C LDR INTEGER 87C The leading dimension of the array R. 88C LDR >= 2*(M+L)*NOBR. 89C 90C A (input or output) DOUBLE PRECISION array, dimension 91C (LDA,N) 92C On entry, if METH = 'N' and JOB = 'B' or 'D', the 93C leading N-by-N part of this array must contain the system 94C state matrix. 95C If METH = 'M' or (METH = 'N' and JOB = 'A' or 'C'), 96C this array need not be set on input. 97C On exit, if JOB = 'A' or 'C' and INFO = 0, the 98C leading N-by-N part of this array contains the system 99C state matrix. 100C 101C LDA INTEGER 102C The leading dimension of the array A. 103C LDA >= N, if JOB = 'A' or 'C', or METH = 'N' and 104C JOB = 'B' or 'D'; 105C LDA >= 1, otherwise. 106C 107C C (input or output) DOUBLE PRECISION array, dimension 108C (LDC,N) 109C On entry, if METH = 'N' and JOB = 'B' or 'D', the 110C leading L-by-N part of this array must contain the system 111C output matrix. 112C If METH = 'M' or (METH = 'N' and JOB = 'A' or 'C'), 113C this array need not be set on input. 114C On exit, if JOB = 'A' or 'C' and INFO = 0, or 115C INFO = 3 (or INFO >= 0, for METH = 'M'), the leading 116C L-by-N part of this array contains the system output 117C matrix. 118C 119C LDC INTEGER 120C The leading dimension of the array C. 121C LDC >= L, if JOB = 'A' or 'C', or METH = 'N' and 122C JOB = 'B' or 'D'; 123C LDC >= 1, otherwise. 124C 125C B (output) DOUBLE PRECISION array, dimension (LDB,M) 126C If M > 0, JOB = 'A', 'B', or 'D' and INFO = 0, the 127C leading N-by-M part of this array contains the system 128C input matrix. If M = 0 or JOB = 'C', this array is 129C not referenced. 130C 131C LDB INTEGER 132C The leading dimension of the array B. 133C LDB >= N, if M > 0 and JOB = 'A', 'B', or 'D'; 134C LDB >= 1, if M = 0 or JOB = 'C'. 135C 136C D (output) DOUBLE PRECISION array, dimension (LDD,M) 137C If M > 0, JOB = 'A' or 'D' and INFO = 0, the leading 138C L-by-M part of this array contains the system input-output 139C matrix. If M = 0 or JOB = 'C' or 'B', this array is 140C not referenced. 141C 142C LDD INTEGER 143C The leading dimension of the array D. 144C LDD >= L, if M > 0 and JOB = 'A' or 'D'; 145C LDD >= 1, if M = 0 or JOB = 'C' or 'B'. 146C 147C Q (output) DOUBLE PRECISION array, dimension (LDQ,N) 148C If JOBCV = 'C', the leading N-by-N part of this array 149C contains the positive semidefinite state covariance matrix 150C to be used as state weighting matrix when computing the 151C Kalman gain. 152C This parameter is not referenced if JOBCV = 'N'. 153C 154C LDQ INTEGER 155C The leading dimension of the array Q. 156C LDQ >= N, if JOBCV = 'C'; 157C LDQ >= 1, if JOBCV = 'N'. 158C 159C RY (output) DOUBLE PRECISION array, dimension (LDRY,L) 160C If JOBCV = 'C', the leading L-by-L part of this array 161C contains the positive (semi)definite output covariance 162C matrix to be used as output weighting matrix when 163C computing the Kalman gain. 164C This parameter is not referenced if JOBCV = 'N'. 165C 166C LDRY INTEGER 167C The leading dimension of the array RY. 168C LDRY >= L, if JOBCV = 'C'; 169C LDRY >= 1, if JOBCV = 'N'. 170C 171C S (output) DOUBLE PRECISION array, dimension (LDS,L) 172C If JOBCV = 'C', the leading N-by-L part of this array 173C contains the state-output cross-covariance matrix to be 174C used as cross-weighting matrix when computing the Kalman 175C gain. 176C This parameter is not referenced if JOBCV = 'N'. 177C 178C LDS INTEGER 179C The leading dimension of the array S. 180C LDS >= N, if JOBCV = 'C'; 181C LDS >= 1, if JOBCV = 'N'. 182C 183C O (output) DOUBLE PRECISION array, dimension ( LDO,N ) 184C If METH = 'M' and JOBCV = 'C', or METH = 'N', 185C the leading L*NOBR-by-N part of this array contains 186C the estimated extended observability matrix, i.e., the 187C first N columns of the relevant singular vectors. 188C If METH = 'M' and JOBCV = 'N', this array is not 189C referenced. 190C 191C LDO INTEGER 192C The leading dimension of the array O. 193C LDO >= L*NOBR, if JOBCV = 'C' or METH = 'N'; 194C LDO >= 1, otherwise. 195C 196C Tolerances 197C 198C TOL DOUBLE PRECISION 199C The tolerance to be used for estimating the rank of 200C matrices. If the user sets TOL > 0, then the given value 201C of TOL is used as a lower bound for the reciprocal 202C condition number; an m-by-n matrix whose estimated 203C condition number is less than 1/TOL is considered to 204C be of full rank. If the user sets TOL <= 0, then an 205C implicitly computed, default tolerance, defined by 206C TOLDEF = m*n*EPS, is used instead, where EPS is the 207C relative machine precision (see LAPACK Library routine 208C DLAMCH). 209C 210C Workspace 211C 212C IWORK INTEGER array, dimension (LIWORK) 213C LIWORK = N, if METH = 'M' and M = 0 214C or JOB = 'C' and JOBCV = 'N'; 215C LIWORK = M*NOBR+N, if METH = 'M', JOB = 'C', 216C and JOBCV = 'C'; 217C LIWORK = max(L*NOBR,M*NOBR), if METH = 'M', JOB <> 'C', 218C and JOBCV = 'N'; 219C LIWORK = max(L*NOBR,M*NOBR+N), if METH = 'M', JOB <> 'C', 220C and JOBCV = 'C'; 221C LIWORK = max(M*NOBR+N,M*(N+L)), if METH = 'N'. 222C 223C DWORK DOUBLE PRECISION array, dimension (LDWORK) 224C On exit, if INFO = 0, DWORK(1) returns the optimal value 225C of LDWORK, and DWORK(2), DWORK(3), DWORK(4), and 226C DWORK(5) contain the reciprocal condition numbers of the 227C triangular factors of the matrices, defined in the code, 228C GaL (GaL = Un(1:(s-1)*L,1:n)), R_1c (if METH = 'M'), 229C M (if JOBCV = 'C' or METH = 'N'), and Q or T (see 230C SLICOT Library routines IB01PY or IB01PX), respectively. 231C If METH = 'N', DWORK(3) is set to one without any 232C calculations. Similarly, if METH = 'M' and JOBCV = 'N', 233C DWORK(4) is set to one. If M = 0 or JOB = 'C', 234C DWORK(3) and DWORK(5) are set to one. 235C On exit, if INFO = -30, DWORK(1) returns the minimum 236C value of LDWORK. 237C 238C LDWORK INTEGER 239C The length of the array DWORK. 240C LDWORK >= max( LDW1,LDW2 ), where, if METH = 'M', 241C LDW1 >= max( 2*(L*NOBR-L)*N+2*N, (L*NOBR-L)*N+N*N+7*N ), 242C if JOB = 'C' or JOB = 'A' and M = 0; 243C LDW1 >= max( 2*(L*NOBR-L)*N+N*N+7*N, 244C (L*NOBR-L)*N+N+6*M*NOBR, (L*NOBR-L)*N+N+ 245C max( L+M*NOBR, L*NOBR + max( 3*L*NOBR, M ))) 246C if M > 0 and JOB = 'A', 'B', or 'D'; 247C LDW2 >= 0, if JOBCV = 'N'; 248C LDW2 >= max( (L*NOBR-L)*N+Aw+2*N+max(5*N,(2*M+L)*NOBR+L), 249C 4*(M*NOBR+N), M*NOBR+2*N+L ), if JOBCV = 'C', 250C where Aw = N+N*N, if M = 0 or JOB = 'C'; 251C Aw = 0, otherwise; 252C and, if METH = 'N', 253C LDW1 >= max( (L*NOBR-L)*N+2*N+(2*M+L)*NOBR+L, 254C 2*(L*NOBR-L)*N+N*N+8*N, N+4*(M*NOBR+N), 255C M*NOBR+3*N+L ); 256C LDW2 >= 0, if M = 0 or JOB = 'C'; 257C LDW2 >= M*NOBR*(N+L)*(M*(N+L)+1)+ 258C max( (N+L)**2, 4*M*(N+L)+1 ), 259C if M > 0 and JOB = 'A', 'B', or 'D'. 260C For good performance, LDWORK should be larger. 261C 262C Warning Indicator 263C 264C IWARN INTEGER 265C = 0: no warning; 266C = 4: a least squares problem to be solved has a 267C rank-deficient coefficient matrix; 268C = 5: the computed covariance matrices are too small. 269C The problem seems to be a deterministic one. 270C 271C Error Indicator 272C 273C INFO INTEGER 274C = 0: successful exit; 275C < 0: if INFO = -i, the i-th argument had an illegal 276C value; 277C = 2: the singular value decomposition (SVD) algorithm did 278C not converge; 279C = 3: a singular upper triangular matrix was found. 280C 281C METHOD 282C 283C In the MOESP approach, the matrices A and C are first 284C computed from an estimated extended observability matrix [1], 285C and then, the matrices B and D are obtained by solving an 286C extended linear system in a least squares sense. 287C In the N4SID approach, besides the estimated extended 288C observability matrix, the solutions of two least squares problems 289C are used to build another least squares problem, whose solution 290C is needed to compute the system matrices A, C, B, and D. The 291C solutions of the two least squares problems are also optionally 292C used by both approaches to find the covariance matrices. 293C 294C REFERENCES 295C 296C [1] Verhaegen M., and Dewilde, P. 297C Subspace Model Identification. Part 1: The output-error state- 298C space model identification class of algorithms. 299C Int. J. Control, 56, pp. 1187-1210, 1992. 300C 301C [2] Van Overschee, P., and De Moor, B. 302C N4SID: Two Subspace Algorithms for the Identification 303C of Combined Deterministic-Stochastic Systems. 304C Automatica, Vol.30, No.1, pp. 75-93, 1994. 305C 306C [3] Van Overschee, P. 307C Subspace Identification : Theory - Implementation - 308C Applications. 309C Ph. D. Thesis, Department of Electrical Engineering, 310C Katholieke Universiteit Leuven, Belgium, Feb. 1995. 311C 312C [4] Sima, V. 313C Subspace-based Algorithms for Multivariable System 314C Identification. 315C Studies in Informatics and Control, 5, pp. 335-344, 1996. 316C 317C NUMERICAL ASPECTS 318C 319C The implemented method is numerically stable. 320C 321C FURTHER COMMENTS 322C 323C In some applications, it is useful to compute the system matrices 324C using two calls to this routine, the first one with JOB = 'C', 325C and the second one with JOB = 'B' or 'D'. This is slightly less 326C efficient than using a single call with JOB = 'A', because some 327C calculations are repeated. If METH = 'N', all the calculations 328C at the first call are performed again at the second call; 329C moreover, it is required to save the needed submatrices of R 330C before the first call and restore them before the second call. 331C If the covariance matrices are desired, JOBCV should be set 332C to 'C' at the second call. If B and D are both needed, they 333C should be computed at once. 334C It is possible to compute the matrices A and C using the MOESP 335C algorithm (METH = 'M'), and the matrices B and D using the N4SID 336C algorithm (METH = 'N'). This combination could be slightly more 337C efficient than N4SID algorithm alone and it could be more accurate 338C than MOESP algorithm. No saving/restoring is needed in such a 339C combination, provided JOBCV is set to 'N' at the first call. 340C Recommended usage: either one call with JOB = 'A', or 341C first call with METH = 'M', JOB = 'C', JOBCV = 'N', 342C second call with METH = 'M', JOB = 'D', JOBCV = 'C', or 343C first call with METH = 'M', JOB = 'C', JOBCV = 'N', 344C second call with METH = 'N', JOB = 'D', JOBCV = 'C'. 345C 346C CONTRIBUTOR 347C 348C V. Sima, Research Institute for Informatics, Bucharest, Dec. 1999. 349C 350C REVISIONS 351C 352C March 2000, Feb. 2001. 353C 354C KEYWORDS 355C 356C Identification methods; least squares solutions; multivariable 357C systems; QR decomposition; singular value decomposition. 358C 359C ****************************************************************** 360C 361C .. Parameters .. 362 DOUBLE PRECISION ZERO, ONE, TWO, THREE 363 PARAMETER ( ZERO = 0.0D0, ONE = 1.0D0, TWO = 2.0D0, 364 $ THREE = 3.0D0 ) 365C .. Scalar Arguments .. 366 DOUBLE PRECISION TOL 367 INTEGER INFO, IWARN, L, LDA, LDB, LDC, LDD, LDO, LDQ, 368 $ LDR, LDRY, LDS, LDWORK, M, N, NOBR, NSMPL 369 CHARACTER JOB, JOBCV, METH 370C .. Array Arguments .. 371 DOUBLE PRECISION A(LDA, *), B(LDB, *), C(LDC, *), D(LDD, *), 372 $ DWORK(*), O(LDO, *), Q(LDQ, *), R(LDR, *), 373 $ RY(LDRY, *), S(LDS, *) 374 INTEGER IWORK( * ) 375C .. Local Scalars .. 376 DOUBLE PRECISION EPS, RCOND1, RCOND2, RCOND3, RCOND4, RNRM, 377 $ SVLMAX, THRESH, TOLL, TOLL1 378 INTEGER I, IAW, ID, IERR, IGAL, IHOUS, ISV, ITAU, 379 $ ITAU1, ITAU2, IU, IUN2, IWARNL, IX, JWORK, 380 $ LDUN2, LDUNN, LDW, LMMNOB, LMMNOL, LMNOBR, 381 $ LNOBR, LNOBRN, MAXWRK, MINWRK, MNOBR, MNOBRN, 382 $ N2, NCOL, NN, NPL, NR, NR2, NR3, NR4, NR4MN, 383 $ NR4PL, NROW, RANK, RANK11, RANKM 384 CHARACTER FACT, JOBP, JOBPY 385 LOGICAL FULLR, MOESP, N4SID, SHIFT, WITHAL, WITHB, 386 $ WITHC, WITHCO, WITHD 387C .. Local Array .. 388 DOUBLE PRECISION SVAL(3) 389C .. External Functions .. 390 LOGICAL LSAME 391 INTEGER ILAENV 392 DOUBLE PRECISION DLAMCH, DLANGE 393 EXTERNAL DLAMCH, DLANGE, ILAENV, LSAME 394C .. External Subroutines .. 395 EXTERNAL DCOPY, DGEMM, DGEQRF, DLACPY, DLASET, DORMQR, 396 $ DSYRK, DTRCON, DTRSM, DTRTRS, IB01PX, IB01PY, 397 $ MA02AD, MA02ED, MB02QY, MB02UD, MB03OD, XERBLA 398C .. Intrinsic Functions .. 399 INTRINSIC DBLE, MAX 400C .. Executable Statements .. 401C 402C Decode the scalar input parameters. 403C 404 MOESP = LSAME( METH, 'M' ) 405 N4SID = LSAME( METH, 'N' ) 406 WITHAL = LSAME( JOB, 'A' ) 407 WITHC = LSAME( JOB, 'C' ) .OR. WITHAL 408 WITHD = LSAME( JOB, 'D' ) .OR. WITHAL 409 WITHB = LSAME( JOB, 'B' ) .OR. WITHD 410 WITHCO = LSAME( JOBCV, 'C' ) 411 MNOBR = M*NOBR 412 LNOBR = L*NOBR 413 LMNOBR = LNOBR + MNOBR 414 LMMNOB = LNOBR + 2*MNOBR 415 MNOBRN = MNOBR + N 416 LNOBRN = LNOBR - N 417 LDUN2 = LNOBR - L 418 LDUNN = LDUN2*N 419 LMMNOL = LMMNOB + L 420 NR = LMNOBR + LMNOBR 421 NPL = N + L 422 N2 = N + N 423 NN = N*N 424 MINWRK = 1 425 IWARN = 0 426 INFO = 0 427C 428C Check the scalar input parameters. 429C 430 IF( .NOT.( MOESP .OR. N4SID ) ) THEN 431 INFO = -1 432 ELSE IF( .NOT.( WITHB .OR. WITHC ) ) THEN 433 INFO = -2 434 ELSE IF( .NOT.( WITHCO .OR. LSAME( JOBCV, 'N' ) ) ) THEN 435 INFO = -3 436 ELSE IF( NOBR.LE.1 ) THEN 437 INFO = -4 438 ELSE IF( N.LE.0 .OR. N.GE.NOBR ) THEN 439 INFO = -5 440 ELSE IF( M.LT.0 ) THEN 441 INFO = -6 442 ELSE IF( L.LE.0 ) THEN 443 INFO = -7 444 ELSE IF( WITHCO .AND. NSMPL.LT.NR ) THEN 445 INFO = -8 446 ELSE IF( LDR.LT.NR ) THEN 447 INFO = -10 448 ELSE IF( LDA.LT.1 .OR. ( ( WITHC .OR. ( WITHB .AND. N4SID ) ) 449 $ .AND. LDA.LT.N ) ) THEN 450 INFO = -12 451 ELSE IF( LDC.LT.1 .OR. ( ( WITHC .OR. ( WITHB .AND. N4SID ) ) 452 $ .AND. LDC.LT.L ) ) THEN 453 INFO = -14 454 ELSE IF( LDB.LT.1 .OR. ( WITHB .AND. LDB.LT.N .AND. M.GT.0 ) ) 455 $ THEN 456 INFO = -16 457 ELSE IF( LDD.LT.1 .OR. ( WITHD .AND. LDD.LT.L .AND. M.GT.0 ) ) 458 $ THEN 459 INFO = -18 460 ELSE IF( LDQ.LT.1 .OR. ( WITHCO .AND. LDQ.LT.N ) ) THEN 461 INFO = -20 462 ELSE IF( LDRY.LT.1 .OR. ( WITHCO .AND. LDRY.LT.L ) ) THEN 463 INFO = -22 464 ELSE IF( LDS.LT.1 .OR. ( WITHCO .AND. LDS.LT.N ) ) THEN 465 INFO = -24 466 ELSE IF( LDO.LT.1 .OR. ( ( WITHCO .OR. N4SID ) .AND. 467 $ LDO.LT.LNOBR ) ) THEN 468 INFO = -26 469 ELSE IF( LDWORK.GE.1 ) THEN 470C 471C Compute workspace. 472C (Note: Comments in the code beginning "Workspace:" describe the 473C minimal amount of workspace needed at that point in the code, 474C as well as the preferred amount for good performance. 475C NB refers to the optimal block size for the immediately 476C following subroutine, as returned by ILAENV.) 477C 478 IAW = 0 479 MINWRK = LDUNN + 4*N 480 MAXWRK = LDUNN + N + N*ILAENV( 1, 'DGEQRF', ' ', LDUN2, N, -1, 481 $ -1 ) 482 IF( MOESP ) THEN 483 ID = 0 484 IF( WITHC ) THEN 485 MINWRK = MAX( MINWRK, 2*LDUNN + N2, LDUNN + NN + 7*N ) 486 MAXWRK = MAX( MAXWRK, 2*LDUNN + N + N*ILAENV( 1, 487 $ 'DORMQR', 'LT', LDUN2, N, N, -1 ) ) 488 END IF 489 ELSE 490 ID = N 491 END IF 492C 493 IF( ( M.GT.0 .AND. WITHB ) .OR. N4SID ) THEN 494 MINWRK = MAX( MINWRK, 2*LDUNN + NN + ID + 7*N ) 495 IF ( MOESP ) 496 $ MINWRK = MAX( MINWRK, LDUNN + N + 6*MNOBR, LDUNN + N + 497 $ MAX( L + MNOBR, LNOBR + MAX( 3*LNOBR, M ) ) 498 $ ) 499 ELSE 500 IF( MOESP ) 501 $ IAW = N + NN 502 END IF 503C 504 IF( N4SID .OR. WITHCO ) THEN 505 MINWRK = MAX( MINWRK, LDUNN + IAW + N2 + MAX( 5*N, LMMNOL ), 506 $ ID + 4*MNOBRN, ID + MNOBRN + NPL ) 507 MAXWRK = MAX( MAXWRK, LDUNN + IAW + N2 + 508 $ MAX( N*ILAENV( 1, 'DGEQRF', ' ', LNOBR, N, -1, 509 $ -1 ), LMMNOB* 510 $ ILAENV( 1, 'DORMQR', 'LT', LNOBR, 511 $ LMMNOB, N, -1 ), LMMNOL* 512 $ ILAENV( 1, 'DORMQR', 'LT', LDUN2, 513 $ LMMNOL, N, -1 ) ), 514 $ ID + N + N*ILAENV( 1, 'DGEQRF', ' ', LMNOBR, 515 $ N, -1, -1 ), 516 $ ID + N + NPL*ILAENV( 1, 'DORMQR', 'LT', 517 $ LMNOBR, NPL, N, -1 ) ) 518 IF( N4SID .AND. ( M.GT.0 .AND. WITHB ) ) 519 $ MINWRK = MAX( MINWRK, MNOBR*NPL*( M*NPL + 1 ) + 520 $ MAX( NPL**2, 4*M*NPL + 1 ) ) 521 END IF 522 MAXWRK = MAX( MINWRK, MAXWRK ) 523C 524 IF ( LDWORK.LT.MINWRK ) THEN 525 INFO = -30 526 DWORK( 1 ) = MINWRK 527 END IF 528 END IF 529C 530C Return if there are illegal arguments. 531C 532 IF( INFO.NE.0 ) THEN 533 CALL XERBLA( 'IB01PD', -INFO ) 534 RETURN 535 END IF 536C 537 NR2 = MNOBR + 1 538 NR3 = LMNOBR + 1 539 NR4 = LMMNOB + 1 540C 541C Set the precision parameters. A threshold value EPS**(2/3) is 542C used for deciding to use pivoting or not, where EPS is the 543C relative machine precision (see LAPACK Library routine DLAMCH). 544C 545 EPS = DLAMCH( 'Precision' ) 546 THRESH = EPS**( TWO/THREE ) 547 SVLMAX = ZERO 548 RCOND4 = ONE 549C 550C Let Un be the matrix of left singular vectors (stored in R_22). 551C Copy un1 = GaL = Un(1:(s-1)*L,1:n) in the workspace. 552C 553 IGAL = 1 554 CALL DLACPY( 'Full', LDUN2, N, R(NR2,NR2), LDR, DWORK(IGAL), 555 $ LDUN2 ) 556C 557C Factor un1 = Q1*[r1' 0]' (' means transposition). 558C Workspace: need L*(NOBR-1)*N+2*N, 559C prefer L*(NOBR-1)*N+N+N*NB. 560C 561 ITAU1 = IGAL + LDUNN 562 JWORK = ITAU1 + N 563 LDW = JWORK 564 CALL DGEQRF( LDUN2, N, DWORK(IGAL), LDUN2, DWORK(ITAU1), 565 $ DWORK(JWORK), LDWORK-JWORK+1, IERR ) 566C 567C Compute the reciprocal of the condition number of r1. 568C Workspace: need L*(NOBR-1)*N+4*N. 569C 570 CALL DTRCON( '1-norm', 'Upper', 'NonUnit', N, DWORK(IGAL), LDUN2, 571 $ RCOND1, DWORK(JWORK), IWORK, INFO ) 572C 573 TOLL1 = TOL 574 IF( TOLL1.LE.ZERO ) 575 $ TOLL1 = NN*EPS 576C 577 IF ( ( M.GT.0 .AND. WITHB ) .OR. N4SID ) THEN 578 JOBP = 'P' 579 IF ( WITHAL ) THEN 580 JOBPY = 'D' 581 ELSE 582 JOBPY = JOB 583 END IF 584 ELSE 585 JOBP = 'N' 586 END IF 587C 588 IF ( MOESP ) THEN 589 NCOL = 0 590 IUN2 = JWORK 591 IF ( WITHC ) THEN 592C 593C Set C = Un(1:L,1:n) and then compute the system matrix A. 594C 595C Set un2 = Un(L+1:L*s,1:n) in DWORK(IUN2). 596C Workspace: need 2*L*(NOBR-1)*N+N. 597C 598 CALL DLACPY( 'Full', L, N, R(NR2,NR2), LDR, C, LDC ) 599 CALL DLACPY( 'Full', LDUN2, N, R(NR2+L,NR2), LDR, 600 $ DWORK(IUN2), LDUN2 ) 601C 602C Note that un1 has already been factored as 603C un1 = Q1*[r1' 0]' and usually (generically, assuming 604C observability) has full column rank. 605C Update un2 <-- Q1'*un2 in DWORK(IUN2) and save its 606C first n rows in A. 607C Workspace: need 2*L*(NOBR-1)*N+2*N; 608C prefer 2*L*(NOBR-1)*N+N+N*NB. 609C 610 JWORK = IUN2 + LDUNN 611 CALL DORMQR( 'Left', 'Transpose', LDUN2, N, N, DWORK(IGAL), 612 $ LDUN2, DWORK(ITAU1), DWORK(IUN2), LDUN2, 613 $ DWORK(JWORK), LDWORK-JWORK+1, IERR ) 614 CALL DLACPY( 'Full', N, N, DWORK(IUN2), LDUN2, A, LDA ) 615 NCOL = N 616 JWORK = IUN2 617 END IF 618C 619 IF ( RCOND1.GT.MAX( TOLL1, THRESH ) ) THEN 620C 621C The triangular factor r1 is considered to be of full rank. 622C Solve for A (if requested), r1*A = un2(1:n,:) in A. 623C 624 IF ( WITHC ) THEN 625 CALL DTRTRS( 'Upper', 'NoTranspose', 'NonUnit', N, N, 626 $ DWORK(IGAL), LDUN2, A, LDA, IERR ) 627 IF ( IERR.GT.0 ) THEN 628 INFO = 3 629 RETURN 630 END IF 631 END IF 632 RANK = N 633 ELSE 634C 635C Rank-deficient triangular factor r1. Use SVD of r1, 636C r1 = U*S*V', also for computing A (if requested) from 637C r1*A = un2(1:n,:). Matrix U is computed in DWORK(IU), 638C and V' overwrites r1. If B is requested, the 639C pseudoinverse of r1 and then of GaL are also computed 640C in R(NR3,NR2). 641C Workspace: need c*L*(NOBR-1)*N+N*N+7*N, 642C where c = 1 if B and D are not needed, 643C c = 2 if B and D are needed; 644C prefer larger. 645C 646 IU = IUN2 647 ISV = IU + NN 648 JWORK = ISV + N 649 IF ( M.GT.0 .AND. WITHB ) THEN 650C 651C Save the elementary reflectors used for computing r1, 652C if B, D are needed. 653C Workspace: need 2*L*(NOBR-1)*N+2*N+N*N. 654C 655 IHOUS = JWORK 656 JWORK = IHOUS + LDUNN 657 CALL DLACPY( 'Lower', LDUN2, N, DWORK(IGAL), LDUN2, 658 $ DWORK(IHOUS), LDUN2 ) 659 ELSE 660 IHOUS = IGAL 661 END IF 662C 663 CALL MB02UD( 'Not factored', 'Left', 'NoTranspose', JOBP, N, 664 $ NCOL, ONE, TOLL1, RANK, DWORK(IGAL), LDUN2, 665 $ DWORK(IU), N, DWORK(ISV), A, LDA, R(NR3,NR2), 666 $ LDR, DWORK(JWORK), LDWORK-JWORK+1, IERR ) 667 IF ( IERR.NE.0 ) THEN 668 INFO = 2 669 RETURN 670 END IF 671 MAXWRK = MAX( MAXWRK, INT( DWORK(JWORK) ) + JWORK - 1 ) 672C 673 IF ( RANK.EQ.0 ) THEN 674 JOBP = 'N' 675 ELSE IF ( M.GT.0 .AND. WITHB ) THEN 676C 677C Compute pinv(GaL) in R(NR3,NR2) if B, D are needed. 678C Workspace: need 2*L*(NOBR-1)*N+N*N+3*N; 679C prefer 2*L*(NOBR-1)*N+N*N+2*N+N*NB. 680C 681 CALL DLASET( 'Full', N, LDUN2-N, ZERO, ZERO, 682 $ R(NR3,NR2+N), LDR ) 683 CALL DORMQR( 'Right', 'Transpose', N, LDUN2, N, 684 $ DWORK(IHOUS), LDUN2, DWORK(ITAU1), 685 $ R(NR3,NR2), LDR, DWORK(JWORK), 686 $ LDWORK-JWORK+1, IERR ) 687 MAXWRK = MAX( MAXWRK, INT( DWORK(JWORK) ) + JWORK - 1 ) 688 IF ( WITHCO ) THEN 689C 690C Save pinv(GaL) in DWORK(IGAL). 691C 692 CALL DLACPY( 'Full', N, LDUN2, R(NR3,NR2), LDR, 693 $ DWORK(IGAL), N ) 694 END IF 695 JWORK = IUN2 696 END IF 697 LDW = JWORK 698 END IF 699C 700 IF ( M.GT.0 .AND. WITHB ) THEN 701C 702C Computation of B and D. 703C 704C Compute the reciprocal of the condition number of R_1c. 705C Workspace: need L*(NOBR-1)*N+N+3*M*NOBR. 706C 707 CALL DTRCON( '1-norm', 'Upper', 'NonUnit', MNOBR, R(NR3,1), 708 $ LDR, RCOND2, DWORK(JWORK), IWORK, IERR ) 709C 710 TOLL = TOL 711 IF( TOLL.LE.ZERO ) 712 $ TOLL = MNOBR*MNOBR*EPS 713C 714C Compute the right hand side and solve for K (in R_23), 715C K*R_1c' = u2'*R_2c', 716C where u2 = Un(:,n+1:L*s), and K is (Ls-n) x ms. 717C 718 CALL DGEMM( 'Transpose', 'Transpose', LNOBRN, MNOBR, LNOBR, 719 $ ONE, R(NR2,NR2+N), LDR, R(1,NR2), LDR, ZERO, 720 $ R(NR2,NR3), LDR ) 721C 722 IF ( RCOND2.GT.MAX( TOLL, THRESH ) ) THEN 723C 724C The triangular factor R_1c is considered to be of full 725C rank. Solve for K, K*R_1c' = u2'*R_2c'. 726C 727 CALL DTRSM( 'Right', 'Upper', 'Transpose', 'Non-unit', 728 $ LNOBRN, MNOBR, ONE, R(NR3,1), LDR, 729 $ R(NR2,NR3), LDR ) 730 ELSE 731C 732C Rank-deficient triangular factor R_1c. Use SVD of R_1c 733C for computing K from K*R_1c' = u2'*R_2c', where 734C R_1c = U1*S1*V1'. Matrix U1 is computed in R_33, 735C and V1' overwrites R_1c. 736C Workspace: need L*(NOBR-1)*N+N+6*M*NOBR; 737C prefer larger. 738C 739 ISV = LDW 740 JWORK = ISV + MNOBR 741 CALL MB02UD( 'Not factored', 'Right', 'Transpose', 742 $ 'No pinv', LNOBRN, MNOBR, ONE, TOLL, RANK11, 743 $ R(NR3,1), LDR, R(NR3,NR3), LDR, DWORK(ISV), 744 $ R(NR2,NR3), LDR, DWORK(JWORK), 1, 745 $ DWORK(JWORK), LDWORK-JWORK+1, IERR ) 746 IF ( IERR.NE.0 ) THEN 747 INFO = 2 748 RETURN 749 END IF 750 MAXWRK = MAX( MAXWRK, INT( DWORK(JWORK) ) + JWORK - 1 ) 751 JWORK = LDW 752 END IF 753C 754C Compute the triangular factor of the structured matrix Q 755C and apply the transformations to the matrix Kexpand, where 756C Q and Kexpand are defined in SLICOT Library routine 757C IB01PY. Compute also the matrices B, D. 758C Workspace: need L*(NOBR-1)*N+N+max(L+M*NOBR,L*NOBR+ 759C max(3*L*NOBR,M)); 760C prefer larger. 761C 762 IF ( WITHCO ) 763 $ CALL DLACPY( 'Full', LNOBR, N, R(NR2,NR2), LDR, O, LDO ) 764 CALL IB01PY( METH, JOBPY, NOBR, N, M, L, RANK, R(NR2,NR2), 765 $ LDR, DWORK(IGAL), LDUN2, DWORK(ITAU1), 766 $ R(NR3,NR2), LDR, R(NR2,NR3), LDR, R(NR4,NR2), 767 $ LDR, R(NR4,NR3), LDR, B, LDB, D, LDD, TOL, 768 $ IWORK, DWORK(JWORK), LDWORK-JWORK+1, IWARN, 769 $ INFO ) 770 IF ( INFO.NE.0 ) 771 $ RETURN 772 MAXWRK = MAX( MAXWRK, INT( DWORK(JWORK) ) + JWORK - 1 ) 773 RCOND4 = DWORK(JWORK+1) 774 IF ( WITHCO ) 775 $ CALL DLACPY( 'Full', LNOBR, N, O, LDO, R(NR2,1), LDR ) 776C 777 ELSE 778 RCOND2 = ONE 779 END IF 780C 781 IF ( .NOT.WITHCO ) THEN 782 RCOND3 = ONE 783 GO TO 30 784 END IF 785 ELSE 786C 787C For N4SID, set RCOND2 to one. 788C 789 RCOND2 = ONE 790 END IF 791C 792C If needed, save the first n columns, representing Gam, of the 793C matrix of left singular vectors, Un, in R_21 and in O. 794C 795 IF ( N4SID .OR. ( WITHC .AND. .NOT.WITHAL ) ) THEN 796 IF ( M.GT.0 ) 797 $ CALL DLACPY( 'Full', LNOBR, N, R(NR2,NR2), LDR, R(NR2,1), 798 $ LDR ) 799 CALL DLACPY( 'Full', LNOBR, N, R(NR2,NR2), LDR, O, LDO ) 800 END IF 801C 802C Computations for covariance matrices, and system matrices (N4SID). 803C Solve the least squares problems Gam*Y = R4(1:L*s,1:(2*m+L)*s), 804C GaL*X = R4(L+1:L*s,:), where 805C GaL = Gam(1:L*(s-1),:), Gam has full column rank, and 806C R4 = [ R_14' R_24' R_34' R_44L' ], R_44L = R_44(1:L,:), as 807C returned by SLICOT Library routine IB01ND. 808C First, find the QR factorization of Gam, Gam = Q*R. 809C Workspace: need L*(NOBR-1)*N+Aw+3*N; 810C prefer L*(NOBR-1)*N+Aw+2*N+N*NB, where 811C Aw = N+N*N, if (M = 0 or JOB = 'C'), rank(r1) < N, 812C and METH = 'M'; 813C Aw = 0, otherwise. 814C 815 ITAU2 = LDW 816 JWORK = ITAU2 + N 817 CALL DGEQRF( LNOBR, N, R(NR2,1), LDR, DWORK(ITAU2), 818 $ DWORK(JWORK), LDWORK-JWORK+1, IERR ) 819C 820C For METH = 'M' or when JOB = 'B' or 'D', transpose 821C [ R_14' R_24' R_34' ]' in the last block-row of R, obtaining Z, 822C and for METH = 'N' and JOB = 'A' or 'C', use the matrix Z 823C already available in the last block-row of R, and then apply 824C the transformations, Z <-- Q'*Z. 825C Workspace: need L*(NOBR-1)*N+Aw+2*N+(2*M+L)*NOBR; 826C prefer L*(NOBR-1)*N+Aw+2*N+(2*M+L)*NOBR*NB. 827C 828 IF ( MOESP .OR. ( WITHB .AND. .NOT. WITHAL ) ) 829 $ CALL MA02AD( 'Full', LMMNOB, LNOBR, R(1,NR4), LDR, R(NR4,1), 830 $ LDR ) 831 CALL DORMQR( 'Left', 'Transpose', LNOBR, LMMNOB, N, R(NR2,1), LDR, 832 $ DWORK(ITAU2), R(NR4,1), LDR, DWORK(JWORK), 833 $ LDWORK-JWORK+1, IERR ) 834C 835C Solve for Y, RY = Z in Z and save the transpose of the 836C solution Y in the second block-column of R. 837C 838 CALL DTRTRS( 'Upper', 'NoTranspose', 'NonUnit', N, LMMNOB, 839 $ R(NR2,1), LDR, R(NR4,1), LDR, IERR ) 840 IF ( IERR.GT.0 ) THEN 841 INFO = 3 842 RETURN 843 END IF 844 CALL MA02AD( 'Full', N, LMMNOB, R(NR4,1), LDR, R(1,NR2), LDR ) 845 NR4MN = NR4 - N 846 NR4PL = NR4 + L 847 NROW = LMMNOL 848C 849C SHIFT is .TRUE. if some columns of R_14 : R_44L should be 850C shifted to the right, to avoid overwriting useful information. 851C 852 SHIFT = M.EQ.0 .AND. LNOBR.LT.N2 853C 854 IF ( RCOND1.GT.MAX( TOLL, THRESH ) ) THEN 855C 856C The triangular factor r1 of GaL (GaL = Q1*r1) is 857C considered to be of full rank. 858C 859C Transpose [ R_14' R_24' R_34' R_44L' ]'(:,L+1:L*s) in the 860C last block-row of R (beginning with the (L+1)-th row), 861C obtaining Z1, and then apply the transformations, 862C Z1 <-- Q1'*Z1. 863C Workspace: need L*(NOBR-1)*N+Aw+2*N+ (2*M+L)*NOBR + L; 864C prefer L*(NOBR-1)*N+Aw+2*N+((2*M+L)*NOBR + L)*NB. 865C 866 CALL MA02AD( 'Full', LMMNOL, LDUN2, R(1,NR4PL), LDR, 867 $ R(NR4PL,1), LDR ) 868 CALL DORMQR( 'Left', 'Transpose', LDUN2, LMMNOL, N, 869 $ DWORK(IGAL), LDUN2, DWORK(ITAU1), R(NR4PL,1), LDR, 870 $ DWORK(JWORK), LDWORK-JWORK+1, IERR ) 871C 872C Solve for X, r1*X = Z1 in Z1, and copy the transpose of X 873C into the last part of the third block-column of R. 874C 875 CALL DTRTRS( 'Upper', 'NoTranspose', 'NonUnit', N, LMMNOL, 876 $ DWORK(IGAL), LDUN2, R(NR4PL,1), LDR, IERR ) 877 IF ( IERR.GT.0 ) THEN 878 INFO = 3 879 RETURN 880 END IF 881C 882 IF ( SHIFT ) THEN 883 NR4MN = NR4 884C 885 DO 10 I = L - 1, 0, -1 886 CALL DCOPY( LMMNOL, R(1,NR4+I), 1, R(1,NR4+N+I), 1 ) 887 10 CONTINUE 888C 889 END IF 890 CALL MA02AD( 'Full', N, LMMNOL, R(NR4PL,1), LDR, R(1,NR4MN), 891 $ LDR ) 892 NROW = 0 893 END IF 894C 895 IF ( N4SID .OR. NROW.GT.0 ) THEN 896C 897C METH = 'N' or rank-deficient triangular factor r1. 898C For METH = 'N', use SVD of r1, r1 = U*S*V', for computing 899C X' from X'*GaL' = Z1', if rank(r1) < N. Matrix U is 900C computed in DWORK(IU) and V' overwrites r1. Then, the 901C pseudoinverse of GaL is determined in R(NR4+L,NR2). 902C For METH = 'M', the pseudoinverse of GaL is already available 903C if M > 0 and B is requested; otherwise, the SVD of r1 is 904C available in DWORK(IU), DWORK(ISV), and DWORK(IGAL). 905C Workspace for N4SID: need 2*L*(NOBR-1)*N+N*N+8*N; 906C prefer larger. 907C 908 IF ( MOESP ) THEN 909 FACT = 'F' 910 IF ( M.GT.0 .AND. WITHB ) 911 $ CALL DLACPY( 'Full', N, LDUN2, DWORK(IGAL), N, 912 $ R(NR4PL,NR2), LDR ) 913 ELSE 914C 915C Save the elementary reflectors used for computing r1. 916C 917 IHOUS = JWORK 918 CALL DLACPY( 'Lower', LDUN2, N, DWORK(IGAL), LDUN2, 919 $ DWORK(IHOUS), LDUN2 ) 920 FACT = 'N' 921 IU = IHOUS + LDUNN 922 ISV = IU + NN 923 JWORK = ISV + N 924 END IF 925C 926 CALL MB02UD( FACT, 'Right', 'Transpose', JOBP, NROW, N, ONE, 927 $ TOLL1, RANK, DWORK(IGAL), LDUN2, DWORK(IU), N, 928 $ DWORK(ISV), R(1,NR4PL), LDR, R(NR4PL,NR2), LDR, 929 $ DWORK(JWORK), LDWORK-JWORK+1, IERR ) 930 IF ( NROW.GT.0 ) THEN 931 IF ( SHIFT ) THEN 932 NR4MN = NR4 933 CALL DLACPY( 'Full', LMMNOL, L, R(1,NR4), LDR, 934 $ R(1,NR4-L), LDR ) 935 CALL DLACPY( 'Full', LMMNOL, N, R(1,NR4PL), LDR, 936 $ R(1,NR4MN), LDR ) 937 CALL DLACPY( 'Full', LMMNOL, L, R(1,NR4-L), LDR, 938 $ R(1,NR4+N), LDR ) 939 ELSE 940 CALL DLACPY( 'Full', LMMNOL, N, R(1,NR4PL), LDR, 941 $ R(1,NR4MN), LDR ) 942 END IF 943 END IF 944C 945 IF ( N4SID ) THEN 946 IF ( IERR.NE.0 ) THEN 947 INFO = 2 948 RETURN 949 END IF 950 MAXWRK = MAX( MAXWRK, INT( DWORK(JWORK) ) + JWORK - 1 ) 951C 952C Compute pinv(GaL) in R(NR4+L,NR2). 953C Workspace: need 2*L*(NOBR-1)*N+3*N; 954C prefer 2*L*(NOBR-1)*N+2*N+N*NB. 955C 956 JWORK = IU 957 CALL DLASET( 'Full', N, LDUN2-N, ZERO, ZERO, R(NR4PL,NR2+N), 958 $ LDR ) 959 CALL DORMQR( 'Right', 'Transpose', N, LDUN2, N, 960 $ DWORK(IHOUS), LDUN2, DWORK(ITAU1), 961 $ R(NR4PL,NR2), LDR, DWORK(JWORK), 962 $ LDWORK-JWORK+1, IERR ) 963 MAXWRK = MAX( MAXWRK, INT( DWORK(JWORK) ) + JWORK - 1 ) 964 END IF 965 END IF 966C 967C For METH = 'N', find part of the solution (corresponding to A 968C and C) and, optionally, for both METH = 'M', or METH = 'N', 969C find the residual of the least squares problem that gives the 970C covariances, M*V = N, where 971C ( R_11 ) 972C M = ( Y' ), N = ( X' R4'(:,1:L) ), V = V(n+m*s, n+L), 973C ( 0 0 ) 974C with M((2*m+L)*s+L, n+m*s), N((2*m+L)*s+L, n+L), R4' being 975C stored in the last block-column of R. The last L rows of M 976C are not explicitly considered. Note that, for efficiency, the 977C last m*s columns of M are in the first positions of arrray R. 978C This permutation does not affect the residual, only the 979C solution. (The solution is not needed for METH = 'M'.) 980C Note that R_11 corresponds to the future outputs for both 981C METH = 'M', or METH = 'N' approaches. (For METH = 'N', the 982C first two block-columns have been interchanged.) 983C For METH = 'N', A and C are obtained as follows: 984C [ A' C' ] = V(m*s+1:m*s+n,:). 985C 986C First, find the QR factorization of Y'(m*s+1:(2*m+L)*s,:) 987C and apply the transformations to the corresponding part of N. 988C Compress the workspace for N4SID by moving the scalar reflectors 989C corresponding to Q. 990C Workspace: need d*N+2*N; 991C prefer d*N+N+N*NB; 992C where d = 0, for MOESP, and d = 1, for N4SID. 993C 994 IF ( MOESP ) THEN 995 ITAU = 1 996 ELSE 997 CALL DCOPY( N, DWORK(ITAU2), 1, DWORK, 1 ) 998 ITAU = N + 1 999 END IF 1000C 1001 JWORK = ITAU + N 1002 CALL DGEQRF( LMNOBR, N, R(NR2,NR2), LDR, DWORK(ITAU), 1003 $ DWORK(JWORK), LDWORK-JWORK+1, IERR ) 1004C 1005C Workspace: need d*N+N+(N+L); 1006C prefer d*N+N+(N+L)*NB. 1007C 1008 CALL DORMQR( 'Left', 'Transpose', LMNOBR, NPL, N, R(NR2,NR2), LDR, 1009 $ DWORK(ITAU), R(NR2,NR4MN), LDR, DWORK(JWORK), 1010 $ LDWORK-JWORK+1, IERR ) 1011C 1012 CALL DLASET( 'Lower', L-1, L-1, ZERO, ZERO, R(NR4+1,NR4), LDR ) 1013C 1014C Now, matrix M with permuted block-columns has been 1015C triangularized. 1016C Compute the reciprocal of the condition number of its 1017C triangular factor in R(1:m*s+n,1:m*s+n). 1018C Workspace: need d*N+3*(M*NOBR+N). 1019C 1020 JWORK = ITAU 1021 CALL DTRCON( '1-norm', 'Upper', 'NonUnit', MNOBRN, R, LDR, RCOND3, 1022 $ DWORK(JWORK), IWORK, INFO ) 1023C 1024 TOLL = TOL 1025 IF( TOLL.LE.ZERO ) 1026 $ TOLL = MNOBRN*MNOBRN*EPS 1027 IF ( RCOND3.GT.MAX( TOLL, THRESH ) ) THEN 1028C 1029C The triangular factor is considered to be of full rank. 1030C Solve for V(m*s+1:m*s+n,:), giving [ A' C' ]. 1031C 1032 FULLR = .TRUE. 1033 RANKM = MNOBRN 1034 IF ( N4SID ) 1035 $ CALL DTRSM( 'Left', 'Upper', 'NoTranspose', 'NonUnit', N, 1036 $ NPL, ONE, R(NR2,NR2), LDR, R(NR2,NR4MN), LDR ) 1037 ELSE 1038 FULLR = .FALSE. 1039C 1040C Use QR factorization (with pivoting). For METH = 'N', save 1041C (and then restore) information about the QR factorization of 1042C Gam, for later use. Note that R_11 could be modified by 1043C MB03OD, but the corresponding part of N is also modified 1044C accordingly. 1045C Workspace: need d*N+4*(M*NOBR+N). 1046C 1047 DO 20 I = 1, MNOBRN 1048 IWORK(I) = 0 1049 20 CONTINUE 1050C 1051 IF ( N4SID .AND. ( M.GT.0 .AND. WITHB ) ) 1052 $ CALL DLACPY( 'Full', LNOBR, N, R(NR2,1), LDR, R(NR4,1), 1053 $ LDR ) 1054 JWORK = ITAU + MNOBRN 1055 CALL DLASET( 'Lower', MNOBRN-1, MNOBRN, ZERO, ZERO, R(2,1), 1056 $ LDR ) 1057 CALL MB03OD( 'QR', MNOBRN, MNOBRN, R, LDR, IWORK, TOLL, 1058 $ SVLMAX, DWORK(ITAU), RANKM, SVAL, DWORK(JWORK), 1059 $ IERR ) 1060C 1061C Workspace: need d*N+M*NOBR+N+N+L; 1062C prefer d*N+M*NOBR+N+(N+L)*NB. 1063C 1064 CALL DORMQR( 'Left', 'Transpose', MNOBRN, NPL, MNOBRN, 1065 $ R, LDR, DWORK(ITAU), R(1,NR4MN), LDR, 1066 $ DWORK(JWORK), LDWORK-JWORK+1, IERR ) 1067 MAXWRK = MAX( MAXWRK, INT( DWORK(JWORK) ) + JWORK - 1 ) 1068 END IF 1069C 1070 IF ( WITHCO ) THEN 1071C 1072C The residual (transposed) of the least squares solution 1073C (multiplied by a matrix with orthogonal columns) is stored 1074C in the rows RANKM+1:(2*m+L)*s+L of V, and it should be 1075C squared-up for getting the covariance matrices. (Generically, 1076C RANKM = m*s+n.) 1077C 1078 RNRM = ONE/DBLE( NSMPL ) 1079 IF ( MOESP ) THEN 1080 CALL DSYRK( 'Upper', 'Transpose', NPL, LMMNOL-RANKM, RNRM, 1081 $ R(RANKM+1,NR4MN), LDR, ZERO, R, LDR ) 1082 CALL DLACPY( 'Upper', N, N, R, LDR, Q, LDQ ) 1083 CALL DLACPY( 'Full', N, L, R(1,N+1), LDR, S, LDS ) 1084 CALL DLACPY( 'Upper', L, L, R(N+1,N+1), LDR, RY, LDRY ) 1085 ELSE 1086 CALL DSYRK( 'Upper', 'Transpose', NPL, LMMNOL-RANKM, RNRM, 1087 $ R(RANKM+1,NR4MN), LDR, ZERO, DWORK(JWORK), NPL ) 1088 CALL DLACPY( 'Upper', N, N, DWORK(JWORK), NPL, Q, LDQ ) 1089 CALL DLACPY( 'Full', N, L, DWORK(JWORK+N*NPL), NPL, S, 1090 $ LDS ) 1091 CALL DLACPY( 'Upper', L, L, DWORK(JWORK+N*(NPL+1)), NPL, RY, 1092 $ LDRY ) 1093 END IF 1094 CALL MA02ED( 'Upper', N, Q, LDQ ) 1095 CALL MA02ED( 'Upper', L, RY, LDRY ) 1096C 1097C Check the magnitude of the residual. 1098C 1099 RNRM = DLANGE( '1-norm', LMMNOL-RANKM, NPL, R(RANKM+1,NR4MN), 1100 $ LDR, DWORK(JWORK) ) 1101 IF ( RNRM.LT.THRESH ) 1102 $ IWARN = 5 1103 END IF 1104C 1105 IF ( N4SID ) THEN 1106 IF ( .NOT.FULLR ) THEN 1107 IWARN = 4 1108C 1109C Compute part of the solution of the least squares problem, 1110C M*V = N, for the rank-deficient problem. 1111C Remark: this computation should not be performed before the 1112C symmetric updating operation above. 1113C Workspace: need M*NOBR+3*N+L; 1114C prefer larger. 1115C 1116 CALL MB03OD( 'No QR', N, N, R(NR2,NR2), LDR, IWORK, TOLL, 1117 $ SVLMAX, DWORK(ITAU), RANKM, SVAL, DWORK(JWORK), 1118 $ IERR ) 1119 CALL MB02QY( N, N, NPL, RANKM, R(NR2,NR2), LDR, IWORK, 1120 $ R(NR2,NR4MN), LDR, DWORK(ITAU+MNOBR), 1121 $ DWORK(JWORK), LDWORK-JWORK+1, INFO ) 1122 MAXWRK = MAX( MAXWRK, INT( DWORK(JWORK) ) + JWORK - 1 ) 1123 JWORK = ITAU 1124 IF ( M.GT.0 .AND. WITHB ) 1125 $ CALL DLACPY( 'Full', LNOBR, N, R(NR4,1), LDR, R(NR2,1), 1126 $ LDR ) 1127 END IF 1128C 1129 IF ( WITHC ) THEN 1130C 1131C Obtain A and C, noting that block-permutations have been 1132C implicitly used. 1133C 1134 CALL MA02AD( 'Full', N, N, R(NR2,NR4MN), LDR, A, LDA ) 1135 CALL MA02AD( 'Full', N, L, R(NR2,NR4MN+N), LDR, C, LDC ) 1136 ELSE 1137C 1138C Use the given A and C. 1139C 1140 CALL MA02AD( 'Full', N, N, A, LDA, R(NR2,NR4MN), LDR ) 1141 CALL MA02AD( 'Full', L, N, C, LDC, R(NR2,NR4MN+N), LDR ) 1142 END IF 1143C 1144 IF ( M.GT.0 .AND. WITHB ) THEN 1145C 1146C Obtain B and D. 1147C First, compute the transpose of the matrix K as 1148C N(1:m*s,:) - M(1:m*s,m*s+1:m*s+n)*[A' C'], in the first 1149C m*s rows of R(1,NR4MN). 1150C 1151 CALL DGEMM ( 'NoTranspose', 'NoTranspose', MNOBR, NPL, N, 1152 $ -ONE, R(1,NR2), LDR, R(NR2,NR4MN), LDR, ONE, 1153 $ R(1,NR4MN), LDR ) 1154C 1155C Denote M = pinv(GaL) and construct 1156C 1157C [ [ A ] -1 ] [ R ] 1158C and L = [ [ ] R 0 ] Q', where Gam = Q * [ ]. 1159C [ [ C ] ] [ 0 ] 1160C 1161C Then, solve the least squares problem. 1162C 1163 CALL DLACPY( 'Full', N, N, A, LDA, R(NR2,NR4), LDR ) 1164 CALL DLACPY( 'Full', L, N, C, LDC, R(NR2+N,NR4), LDR ) 1165 CALL DTRSM( 'Right', 'Upper', 'NoTranspose', 'NonUnit', 1166 $ NPL, N, ONE, R(NR2,1), LDR, R(NR2,NR4), LDR ) 1167 CALL DLASET( 'Full', NPL, LNOBRN, ZERO, ZERO, R(NR2,NR4+N), 1168 $ LDR ) 1169C 1170C Workspace: need 2*N+L; prefer N + (N+L)*NB. 1171C 1172 CALL DORMQR( 'Right', 'Transpose', NPL, LNOBR, N, R(NR2,1), 1173 $ LDR, DWORK, R(NR2,NR4), LDR, DWORK(JWORK), 1174 $ LDWORK-JWORK+1, IERR ) 1175C 1176C Obtain the matrix K by transposition, and find B and D. 1177C Workspace: need NOBR*(M*(N+L))**2+M*NOBR*(N+L)+ 1178C max((N+L)**2,4*M*(N+L)+1); 1179C prefer larger. 1180C 1181 CALL MA02AD( 'Full', MNOBR, NPL, R(1,NR4MN), LDR, 1182 $ R(NR2,NR3), LDR ) 1183 IX = MNOBR*NPL**2*M + 1 1184 JWORK = IX + MNOBR*NPL 1185 CALL IB01PX( JOBPY, NOBR, N, M, L, R, LDR, O, LDO, 1186 $ R(NR2,NR4), LDR, R(NR4PL,NR2), LDR, R(NR2,NR3), 1187 $ LDR, DWORK, MNOBR*NPL, DWORK(IX), B, LDB, D, 1188 $ LDD, TOL, IWORK, DWORK(JWORK), LDWORK-JWORK+1, 1189 $ IWARNL, INFO ) 1190 IF ( INFO.NE.0 ) 1191 $ RETURN 1192 IWARN = MAX( IWARN, IWARNL ) 1193 MAXWRK = MAX( MAXWRK, INT( DWORK(JWORK) ) + JWORK - 1 ) 1194 RCOND4 = DWORK(JWORK+1) 1195C 1196 END IF 1197 END IF 1198C 1199 30 CONTINUE 1200C 1201C Return optimal workspace in DWORK(1) and reciprocal condition 1202C numbers in the next locations. 1203C 1204 DWORK(1) = MAXWRK 1205 DWORK(2) = RCOND1 1206 DWORK(3) = RCOND2 1207 DWORK(4) = RCOND3 1208 DWORK(5) = RCOND4 1209 RETURN 1210C 1211C *** Last line of IB01PD *** 1212 END 1213