1 SUBROUTINE AB13MD( FACT, N, Z, LDZ, M, NBLOCK, ITYPE, X, BOUND, D, 2 $ G, IWORK, DWORK, LDWORK, ZWORK, LZWORK, INFO ) 3C 4C RELEASE 4.0, WGS COPYRIGHT 2000. 5C 6C PURPOSE 7C 8C To compute an upper bound on the structured singular value for a 9C given square complex matrix and a given block structure of the 10C uncertainty. 11C 12C ARGUMENTS 13C 14C Mode Parameters 15C 16C FACT CHARACTER*1 17C Specifies whether or not an information from the 18C previous call is supplied in the vector X. 19C = 'F': On entry, X contains information from the 20C previous call. 21C = 'N': On entry, X does not contain an information from 22C the previous call. 23C 24C Input/Output Parameters 25C 26C N (input) INTEGER 27C The order of the matrix Z. N >= 0. 28C 29C Z (input) COMPLEX*16 array, dimension (LDZ,N) 30C The leading N-by-N part of this array must contain the 31C complex matrix Z for which the upper bound on the 32C structured singular value is to be computed. 33C 34C LDZ INTEGER 35C The leading dimension of the array Z. LDZ >= max(1,N). 36C 37C M (input) INTEGER 38C The number of diagonal blocks in the block structure of 39C the uncertainty. M >= 1. 40C 41C NBLOCK (input) INTEGER array, dimension (M) 42C The vector of length M containing the block structure 43C of the uncertainty. NBLOCK(I), I = 1:M, is the size of 44C each block. 45C 46C ITYPE (input) INTEGER array, dimension (M) 47C The vector of length M indicating the type of each block. 48C For I = 1:M, 49C ITYPE(I) = 1 indicates that the corresponding block is a 50C real block, and 51C ITYPE(I) = 2 indicates that the corresponding block is a 52C complex block. 53C NBLOCK(I) must be equal to 1 if ITYPE(I) is equal to 1. 54C 55C X (input/output) DOUBLE PRECISION array, dimension 56C ( M + MR - 1 ), where MR is the number of the real blocks. 57C On entry, if FACT = 'F' and NBLOCK(1) < N, this array 58C must contain information from the previous call to AB13MD. 59C If NBLOCK(1) = N, this array is not used. 60C On exit, if NBLOCK(1) < N, this array contains information 61C that can be used in the next call to AB13MD for a matrix 62C close to Z. 63C 64C BOUND (output) DOUBLE PRECISION 65C The upper bound on the structured singular value. 66C 67C D, G (output) DOUBLE PRECISION arrays, dimension (N) 68C The vectors of length N containing the diagonal entries 69C of the diagonal N-by-N matrices D and G, respectively, 70C such that the matrix 71C Z'*D^2*Z + sqrt(-1)*(G*Z-Z'*G) - BOUND^2*D^2 72C is negative semidefinite. 73C 74C Workspace 75C 76C IWORK INTEGER array, dimension MAX(4*M-2,N) 77C 78C DWORK DOUBLE PRECISION array, dimension (LDWORK) 79C On exit, if INFO = 0, DWORK(1) contains the optimal value 80C of LDWORK. 81C 82C LDWORK INTEGER 83C The dimension of the array DWORK. 84C LDWORK >= 2*N*N*M - N*N + 9*M*M + N*M + 11*N + 33*M - 11. 85C For best performance 86C LDWORK >= 2*N*N*M - N*N + 9*M*M + N*M + 6*N + 33*M - 11 + 87C MAX( 5*N,2*N*NB ) 88C where NB is the optimal blocksize returned by ILAENV. 89C 90C ZWORK COMPLEX*16 array, dimension (LZWORK) 91C On exit, if INFO = 0, ZWORK(1) contains the optimal value 92C of LZWORK. 93C 94C LZWORK INTEGER 95C The dimension of the array ZWORK. 96C LZWORK >= 6*N*N*M + 12*N*N + 6*M + 6*N - 3. 97C For best performance 98C LZWORK >= 6*N*N*M + 12*N*N + 6*M + 3*N - 3 + 99C MAX( 3*N,N*NB ) 100C where NB is the optimal blocksize returned by ILAENV. 101C 102C Error Indicator 103C 104C INFO INTEGER 105C = 0: successful exit; 106C < 0: if INFO = -i, the i-th argument had an illegal 107C value; 108C = 1: the block sizes must be positive integers; 109C = 2: the sum of block sizes must be equal to N; 110C = 3: the size of a real block must be equal to 1; 111C = 4: the block type must be either 1 or 2; 112C = 5: errors in solving linear equations or in matrix 113C inversion; 114C = 6: errors in computing eigenvalues or singular values. 115C 116C METHOD 117C 118C The routine computes the upper bound proposed in [1]. 119C 120C REFERENCES 121C 122C [1] Fan, M.K.H., Tits, A.L., and Doyle, J.C. 123C Robustness in the presence of mixed parametric uncertainty 124C and unmodeled dynamics. 125C IEEE Trans. Automatic Control, vol. AC-36, 1991, pp. 25-38. 126C 127C NUMERICAL ASPECTS 128C 129C The accuracy and speed of computation depend on the value of 130C the internal threshold TOL. 131C 132C CONTRIBUTORS 133C 134C P.Hr. Petkov, F. Delebecque, D.W. Gu, M.M. Konstantinov and 135C S. Steer with the assistance of V. Sima, September 2000. 136C 137C REVISIONS 138C 139C V. Sima, Katholieke Universiteit Leuven, February 2001. 140C 141C KEYWORDS 142C 143C H-infinity optimal control, Robust control, Structured singular 144C value. 145C 146C ****************************************************************** 147C 148C .. Parameters .. 149 COMPLEX*16 CZERO, CONE, CIMAG 150 PARAMETER ( CZERO = ( 0.0D+0, 0.0D+0 ), 151 $ CONE = ( 1.0D+0, 0.0D+0 ), 152 $ CIMAG = ( 0.0D+0, 1.0D+0 ) ) 153 DOUBLE PRECISION ZERO, ONE, TWO, FOUR, FIVE, EIGHT, TEN, FORTY, 154 $ FIFTY 155 PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0, TWO = 2.0D+0, 156 $ FOUR = 4.0D+0, FIVE = 5.0D+0, EIGHT = 8.0D+0, 157 $ TEN = 1.0D+1, FORTY = 4.0D+1, FIFTY = 5.0D+1 158 $ ) 159 DOUBLE PRECISION ALPHA, BETA, THETA 160 PARAMETER ( ALPHA = 100.0D+0, BETA = 1.0D-2, 161 $ THETA = 1.0D-2 ) 162 DOUBLE PRECISION C1, C2, C3, C4, C5, C6, C7, C8, C9 163 PARAMETER ( C1 = 1.0D-3, C2 = 1.0D-2, C3 = 0.25D+0, 164 $ C4 = 0.9D+0, C5 = 1.5D+0, C6 = 1.0D+1, 165 $ C7 = 1.0D+2, C8 = 1.0D+3, C9 = 1.0D+4 ) 166C .. 167C .. Scalar Arguments .. 168 CHARACTER FACT 169 INTEGER INFO, LDWORK, LDZ, LZWORK, M, N 170 DOUBLE PRECISION BOUND 171C .. 172C .. Array Arguments .. 173 INTEGER ITYPE( * ), IWORK( * ), NBLOCK( * ) 174 COMPLEX*16 Z( LDZ, * ), ZWORK( * ) 175 DOUBLE PRECISION D( * ), DWORK( * ), G( * ), X( * ) 176C .. 177C .. Local Scalars .. 178 INTEGER I, INFO2, ISUM, ITER, IW2, IW3, IW4, IW5, IW6, 179 $ IW7, IW8, IW9, IW10, IW11, IW12, IW13, IW14, 180 $ IW15, IW16, IW17, IW18, IW19, IW20, IW21, IW22, 181 $ IW23, IW24, IW25, IW26, IW27, IW28, IW29, IW30, 182 $ IW31, IW32, IW33, IWRK, IZ2, IZ3, IZ4, IZ5, 183 $ IZ6, IZ7, IZ8, IZ9, IZ10, IZ11, IZ12, IZ13, 184 $ IZ14, IZ15, IZ16, IZ17, IZ18, IZ19, IZ20, IZ21, 185 $ IZ22, IZ23, IZ24, IZWRK, J, K, L, LWA, LWAMAX, 186 $ LZA, LZAMAX, MINWRK, MINZRK, MR, MT, NSUM, SDIM 187 COMPLEX*16 DETF, TEMPIJ, TEMPJI 188 DOUBLE PRECISION C, COLSUM, DELTA, DLAMBD, E, EMAX, EMIN, EPS, 189 $ HN, HNORM, HNORM1, PHI, PP, PROD, RAT, RCOND, 190 $ REGPAR, ROWSUM, SCALE, SNORM, STSIZE, SVLAM, 191 $ T1, T2, T3, TAU, TEMP, TOL, TOL2, TOL3, TOL4, 192 $ TOL5, YNORM1, YNORM2, ZNORM, ZNORM2 193 LOGICAL GTEST, POS, SELECT, XFACT 194C .. 195C .. Local Arrays .. 196 LOGICAL BWORK( 1 ) 197C .. 198C .. External Functions 199 DOUBLE PRECISION DDOT, DLAMCH, DLANGE, ZLANGE 200 LOGICAL LSAME 201 EXTERNAL DDOT, DLAMCH, DLANGE, LSAME, ZLANGE 202C .. 203C .. External Subroutines .. 204 EXTERNAL DCOPY, DGEMV, DLACPY, DLASET, DSCAL, DSYCON, 205 $ DSYSV, DSYTRF, DSYTRS, XERBLA, ZCOPY, ZGEES, 206 $ ZGEMM, ZGEMV, ZGESVD, ZGETRF, ZGETRI, ZLACPY, 207 $ ZLASCL 208C .. 209C .. Intrinsic Functions .. 210 INTRINSIC ABS, DCMPLX, DCONJG, DFLOAT, DREAL, INT, LOG, 211 $ MAX, SQRT 212C .. 213C .. Executable Statements .. 214C 215C Compute workspace. 216C 217 MINWRK = 2*N*N*M - N*N + 9*M*M + N*M + 11*N + 33*M - 11 218 MINZRK = 6*N*N*M + 12*N*N + 6*M + 6*N - 3 219C 220C Decode and Test input parameters. 221C 222 INFO = 0 223 XFACT = LSAME( FACT, 'F' ) 224 IF( .NOT.XFACT .AND. .NOT.LSAME( FACT, 'N' ) ) THEN 225 INFO = -1 226 ELSE IF( N.LT.0 ) THEN 227 INFO = -2 228 ELSE IF( LDZ.LT.MAX( 1, N ) ) THEN 229 INFO = -4 230 ELSE IF( M.LT.1 ) THEN 231 INFO = -5 232 ELSE IF( LDWORK.LT.MINWRK ) THEN 233 INFO = -14 234 ELSE IF( LZWORK.LT.MINZRK ) THEN 235 INFO = -16 236 END IF 237 IF( INFO.NE.0 ) THEN 238 CALL XERBLA( 'AB13MD', -INFO ) 239 RETURN 240 END IF 241C 242 NSUM = 0 243 ISUM = 0 244 MR = 0 245 DO 10 I = 1, M 246 IF( NBLOCK( I ).LT.1 ) THEN 247 INFO = 1 248 RETURN 249 END IF 250 IF( ITYPE( I ).EQ.1 .AND. NBLOCK( I ).GT.1 ) THEN 251 INFO = 3 252 RETURN 253 END IF 254 NSUM = NSUM + NBLOCK( I ) 255 IF( ITYPE( I ).EQ.1 ) MR = MR + 1 256 IF( ITYPE( I ).EQ.1 .OR. ITYPE( I ).EQ.2 ) ISUM = ISUM + 1 257 10 CONTINUE 258 IF( NSUM.NE.N ) THEN 259 INFO = 2 260 RETURN 261 END IF 262 IF( ISUM.NE.M ) THEN 263 INFO = 4 264 RETURN 265 END IF 266 MT = M + MR - 1 267C 268 LWAMAX = 0 269 LZAMAX = 0 270C 271C Set D = In, G = 0. 272C 273 CALL DLASET( 'Full', N, 1, ONE, ONE, D, N ) 274 CALL DLASET( 'Full', N, 1, ZERO, ZERO, G, N ) 275C 276C Quick return if possible. 277C 278 ZNORM = ZLANGE( 'F', N, N, Z, LDZ, DWORK ) 279 IF( ZNORM.EQ.ZERO ) THEN 280 BOUND = ZERO 281 DWORK( 1 ) = ONE 282 ZWORK( 1 ) = CONE 283 RETURN 284 END IF 285C 286C Copy Z into ZWORK. 287C 288 CALL ZLACPY( 'Full', N, N, Z, LDZ, ZWORK, N ) 289C 290C Exact bound for the case NBLOCK( 1 ) = N. 291C 292 IF( NBLOCK( 1 ).EQ.N ) THEN 293 IF( ITYPE( 1 ).EQ.1 ) THEN 294C 295C 1-by-1 real block. 296C 297 BOUND = ZERO 298 DWORK( 1 ) = ONE 299 ZWORK( 1 ) = CONE 300 ELSE 301C 302C N-by-N complex block. 303C 304 CALL ZGESVD( 'N', 'N', N, N, ZWORK, N, DWORK, ZWORK, 1, 305 $ ZWORK, 1, ZWORK( N*N+1 ), LZWORK, 306 $ DWORK( N+1 ), INFO2 ) 307 IF( INFO2.GT.0 ) THEN 308 INFO = 6 309 RETURN 310 END IF 311 BOUND = DWORK( 1 ) 312 LZA = N*N + INT( ZWORK( N*N+1 ) ) 313 DWORK( 1 ) = 5*N 314 ZWORK( 1 ) = DCMPLX( LZA ) 315 END IF 316 RETURN 317 END IF 318C 319C Get machine precision. 320C 321 EPS = DLAMCH( 'P' ) 322C 323C Set tolerances. 324C 325 TOL = C7*SQRT( EPS ) 326 TOL2 = C9*EPS 327 TOL3 = C6*EPS 328 TOL4 = C1 329 TOL5 = C1 330 REGPAR = C8*EPS 331C 332C Real workspace usage. 333C 334 IW2 = M*M 335 IW3 = IW2 + M 336 IW4 = IW3 + N 337 IW5 = IW4 + M 338 IW6 = IW5 + M 339 IW7 = IW6 + N 340 IW8 = IW7 + N 341 IW9 = IW8 + N*( M - 1 ) 342 IW10 = IW9 + N*N*MT 343 IW11 = IW10 + MT 344 IW12 = IW11 + MT*MT 345 IW13 = IW12 + N 346 IW14 = IW13 + MT + 1 347 IW15 = IW14 + MT + 1 348 IW16 = IW15 + MT + 1 349 IW17 = IW16 + MT + 1 350 IW18 = IW17 + MT + 1 351 IW19 = IW18 + MT 352 IW20 = IW19 + MT 353 IW21 = IW20 + MT 354 IW22 = IW21 + N 355 IW23 = IW22 + M - 1 356 IW24 = IW23 + MR 357 IW25 = IW24 + N 358 IW26 = IW25 + 2*MT 359 IW27 = IW26 + MT 360 IW28 = IW27 + MT 361 IW29 = IW28 + M - 1 362 IW30 = IW29 + MR 363 IW31 = IW30 + N + 2*MT 364 IW32 = IW31 + MT*MT 365 IW33 = IW32 + MT 366 IWRK = IW33 + MT + 1 367C 368C Double complex workspace usage. 369C 370 IZ2 = N*N 371 IZ3 = IZ2 + N*N 372 IZ4 = IZ3 + N*N 373 IZ5 = IZ4 + N*N 374 IZ6 = IZ5 + N*N 375 IZ7 = IZ6 + N*N*MT 376 IZ8 = IZ7 + N*N 377 IZ9 = IZ8 + N*N 378 IZ10 = IZ9 + N*N 379 IZ11 = IZ10 + MT 380 IZ12 = IZ11 + N*N 381 IZ13 = IZ12 + N 382 IZ14 = IZ13 + N*N 383 IZ15 = IZ14 + N 384 IZ16 = IZ15 + N*N 385 IZ17 = IZ16 + N 386 IZ18 = IZ17 + N*N 387 IZ19 = IZ18 + N*N*MT 388 IZ20 = IZ19 + MT 389 IZ21 = IZ20 + N*N*MT 390 IZ22 = IZ21 + N*N 391 IZ23 = IZ22 + N*N 392 IZ24 = IZ23 + N*N 393 IZWRK = IZ24 + MT 394C 395C Compute the cummulative sums of blocks dimensions. 396C 397 IWORK( 1 ) = 0 398 DO 20 I = 2, M+1 399 IWORK( I ) = IWORK( I - 1 ) + NBLOCK( I - 1 ) 400 20 CONTINUE 401C 402C Find Osborne scaling if initial scaling is not given. 403C 404 IF( .NOT.XFACT ) THEN 405 CALL DLASET( 'Full', M, M, ZERO, ZERO, DWORK, M ) 406 CALL DLASET( 'Full', M, 1, ONE, ONE, DWORK( IW2+1 ), M ) 407 ZNORM = ZLANGE( 'F', N, N, ZWORK, N, DWORK ) 408 DO 40 J = 1, M 409 DO 30 I = 1, M 410 IF( I.NE.J ) THEN 411 CALL ZLACPY( 'Full', IWORK( I+1 )-IWORK( I ), 412 $ IWORK( J+1 )-IWORK( J ), 413 $ Z( IWORK( I )+1, IWORK( J )+1 ), LDZ, 414 $ ZWORK( IZ2+1 ), N ) 415 CALL ZGESVD( 'N', 'N', IWORK( I+1 )-IWORK( I ), 416 $ IWORK( J+1 )-IWORK( J ), ZWORK( IZ2+1 ), 417 $ N, DWORK( IW3+1 ), ZWORK, 1, ZWORK, 1, 418 $ ZWORK( IZWRK+1 ), LZWORK-IZWRK, 419 $ DWORK( IWRK+1 ), INFO2 ) 420 IF( INFO2.GT.0 ) THEN 421 INFO = 6 422 RETURN 423 END IF 424 LZA = INT( ZWORK( IZWRK+1 ) ) 425 LZAMAX = MAX( LZA, LZAMAX ) 426 ZNORM2 = DWORK( IW3+1 ) 427 DWORK( I+(J-1)*M ) = ZNORM2 + ZNORM*TOL2 428 END IF 429 30 CONTINUE 430 40 CONTINUE 431 CALL DLASET( 'Full', M, 1, ZERO, ZERO, DWORK( IW4+1 ), M ) 432 50 DO 60 I = 1, M 433 DWORK( IW5+I ) = DWORK( IW4+I ) - ONE 434 60 CONTINUE 435 HNORM = DLANGE( 'F', M, 1, DWORK( IW5+1 ), M, DWORK ) 436 IF( HNORM.LE.TOL2 ) GO TO 120 437 DO 110 K = 1, M 438 COLSUM = ZERO 439 DO 70 I = 1, M 440 COLSUM = COLSUM + DWORK( I+(K-1)*M ) 441 70 CONTINUE 442 ROWSUM = ZERO 443 DO 80 J = 1, M 444 ROWSUM = ROWSUM + DWORK( K+(J-1)*M ) 445 80 CONTINUE 446 RAT = SQRT( COLSUM / ROWSUM ) 447 DWORK( IW4+K ) = RAT 448 DO 90 I = 1, M 449 DWORK( I+(K-1)*M ) = DWORK( I+(K-1)*M ) / RAT 450 90 CONTINUE 451 DO 100 J = 1, M 452 DWORK( K+(J-1)*M ) = DWORK( K+(J-1)*M )*RAT 453 100 CONTINUE 454 DWORK( IW2+K ) = DWORK( IW2+K )*RAT 455 110 CONTINUE 456 GO TO 50 457 120 SCALE = ONE / DWORK( IW2+1 ) 458 CALL DSCAL( M, SCALE, DWORK( IW2+1 ), 1 ) 459 ELSE 460 DWORK( IW2+1 ) = ONE 461 DO 130 I = 2, M 462 DWORK( IW2+I ) = SQRT( X( I-1 ) ) 463 130 CONTINUE 464 END IF 465 DO 150 J = 1, M 466 DO 140 I = 1, M 467 IF( I.NE.J ) THEN 468 CALL ZLASCL( 'G', M, M, DWORK( IW2+J ), DWORK( IW2+I ), 469 $ IWORK( I+1 )-IWORK( I ), 470 $ IWORK( J+1 )-IWORK( J ), 471 $ ZWORK( IWORK( I )+1+IWORK( J )*N ), N, 472 $ INFO2 ) 473 END IF 474 140 CONTINUE 475 150 CONTINUE 476C 477C Scale Z by its 2-norm. 478C 479 CALL ZLACPY( 'Full', N, N, ZWORK, N, ZWORK( IZ2+1 ), N ) 480 CALL ZGESVD( 'N', 'N', N, N, ZWORK( IZ2+1 ), N, DWORK( IW3+1 ), 481 $ ZWORK, 1, ZWORK, 1, ZWORK( IZWRK+1 ), LZWORK-IZWRK, 482 $ DWORK( IWRK+1 ), INFO2 ) 483 IF( INFO2.GT.0 ) THEN 484 INFO = 6 485 RETURN 486 END IF 487 LZA = INT( ZWORK( IZWRK+1 ) ) 488 LZAMAX = MAX( LZA, LZAMAX ) 489 ZNORM = DWORK( IW3+1 ) 490 CALL ZLASCL( 'G', M, M, ZNORM, ONE, N, N, ZWORK, N, INFO2 ) 491C 492C Set BB. 493C 494 CALL DLASET( 'Full', N*N, MT, ZERO, ZERO, DWORK( IW9+1 ), N*N ) 495C 496C Set P. 497C 498 DO 160 I = 1, NBLOCK( 1 ) 499 DWORK( IW6+I ) = ONE 500 160 CONTINUE 501 DO 170 I = NBLOCK( 1 )+1, N 502 DWORK( IW6+I ) = ZERO 503 170 CONTINUE 504C 505C Compute P*Z. 506C 507 DO 190 J = 1, N 508 DO 180 I = 1, N 509 ZWORK( IZ3+I+(J-1)*N ) = DCMPLX( DWORK( IW6+I ) )* 510 $ ZWORK( I+(J-1)*N ) 511 180 CONTINUE 512 190 CONTINUE 513C 514C Compute Z'*P*Z. 515C 516 CALL ZGEMM( 'C', 'N', N, N, N, CONE, ZWORK, N, ZWORK( IZ3+1 ), N, 517 $ CZERO, ZWORK( IZ4+1 ), N ) 518C 519C Copy Z'*P*Z into A0. 520C 521 CALL ZLACPY( 'Full', N, N, ZWORK( IZ4+1 ), N, ZWORK( IZ5+1 ), N ) 522C 523C Copy diag(P) into B0d. 524C 525 CALL DCOPY( N, DWORK( IW6+1 ), 1, DWORK( IW7+1 ), 1 ) 526C 527 DO 270 K = 2, M 528C 529C Set P. 530C 531 DO 200 I = 1, IWORK( K ) 532 DWORK( IW6+I ) = ZERO 533 200 CONTINUE 534 DO 210 I = IWORK( K )+1, IWORK( K )+NBLOCK( K ) 535 DWORK( IW6+I ) = ONE 536 210 CONTINUE 537 IF( K.LT.M ) THEN 538 DO 220 I = IWORK( K+1 )+1, N 539 DWORK( IW6+I ) = ZERO 540 220 CONTINUE 541 END IF 542C 543C Compute P*Z. 544C 545 DO 240 J = 1, N 546 DO 230 I = 1, N 547 ZWORK( IZ3+I+(J-1)*N ) = DCMPLX( DWORK( IW6+I ) )* 548 $ ZWORK( I+(J-1)*N ) 549 230 CONTINUE 550 240 CONTINUE 551C 552C Compute t = Z'*P*Z. 553C 554 CALL ZGEMM( 'C', 'N', N, N, N, CONE, ZWORK, N, ZWORK( IZ3+1 ), 555 $ N, CZERO, ZWORK( IZ4+1 ), N ) 556C 557C Copy t(:) into the (k-1)-th column of AA. 558C 559 CALL ZCOPY( N*N, ZWORK( IZ4+1 ), 1, ZWORK( IZ6+1+(K-2)*N*N ), 560 $ 1 ) 561C 562C Copy diag(P) into the (k-1)-th column of BBd. 563C 564 CALL DCOPY( N, DWORK( IW6+1 ), 1, DWORK( IW8+1+(K-2)*N ), 1 ) 565C 566C Copy P(:) into the (k-1)-th column of BB. 567C 568 DO 260 I = 1, N 569 DWORK( IW9+I+(I-1)*N+(K-2)*N*N ) = DWORK( IW6+I ) 570 260 CONTINUE 571 270 CONTINUE 572C 573 L = 0 574C 575 DO 350 K = 1, M 576 IF( ITYPE( K ).EQ.1 ) THEN 577 L = L + 1 578C 579C Set P. 580C 581 DO 280 I = 1, IWORK( K ) 582 DWORK( IW6+I ) = ZERO 583 280 CONTINUE 584 DO 290 I = IWORK( K )+1, IWORK( K )+NBLOCK( K ) 585 DWORK( IW6+I ) = ONE 586 290 CONTINUE 587 IF( K.LT.M ) THEN 588 DO 300 I = IWORK( K+1 )+1, N 589 DWORK( IW6+I ) = ZERO 590 300 CONTINUE 591 END IF 592C 593C Compute P*Z. 594C 595 DO 320 J = 1, N 596 DO 310 I = 1, N 597 ZWORK( IZ3+I+(J-1)*N ) = DCMPLX( DWORK( IW6+I ) )* 598 $ ZWORK( I+(J-1)*N ) 599 310 CONTINUE 600 320 CONTINUE 601C 602C Compute t = sqrt(-1)*( P*Z - Z'*P ). 603C 604 DO 340 J = 1, N 605 DO 330 I = 1, J 606 TEMPIJ = ZWORK( IZ3+I+(J-1)*N ) 607 TEMPJI = ZWORK( IZ3+J+(I-1)*N ) 608 ZWORK( IZ4+I+(J-1)*N ) = CIMAG*( TEMPIJ - 609 $ DCONJG( TEMPJI ) ) 610 ZWORK( IZ4+J+(I-1)*N ) = CIMAG*( TEMPJI - 611 $ DCONJG( TEMPIJ ) ) 612 330 CONTINUE 613 340 CONTINUE 614C 615C Copy t(:) into the (m-1+l)-th column of AA. 616C 617 CALL ZCOPY( N*N, ZWORK( IZ4+1 ), 1, 618 $ ZWORK( IZ6+1+(M-2+L)*N*N ), 1 ) 619 END IF 620 350 CONTINUE 621C 622C Set initial X. 623C 624 DO 360 I = 1, M - 1 625 X( I ) = ONE 626 360 CONTINUE 627 IF( MR.GT.0 ) THEN 628 IF( .NOT.XFACT ) THEN 629 DO 370 I = 1, MR 630 X( M-1+I ) = ZERO 631 370 CONTINUE 632 ELSE 633 L = 0 634 DO 380 K = 1, M 635 IF( ITYPE( K ).EQ.1 ) THEN 636 L = L + 1 637 X( M-1+L ) = X( M-1+L ) / DWORK( IW2+K )**2 638 END IF 639 380 CONTINUE 640 END IF 641 END IF 642C 643C Set constants. 644C 645 SVLAM = ONE / EPS 646 C = ONE 647C 648C Set H. 649C 650 CALL DLASET( 'Full', MT, MT, ZERO, ONE, DWORK( IW11+1 ), MT ) 651C 652 ITER = -1 653C 654C Main iteration loop. 655C 656 390 ITER = ITER + 1 657C 658C Compute A(:) = A0 + AA*x. 659C 660 DO 400 I = 1, MT 661 ZWORK( IZ10+I ) = DCMPLX( X( I ) ) 662 400 CONTINUE 663 CALL ZCOPY( N*N, ZWORK( IZ5+1 ), 1, ZWORK( IZ7+1 ), 1 ) 664 CALL ZGEMV( 'N', N*N, MT, CONE, ZWORK( IZ6+1 ), N*N, 665 $ ZWORK( IZ10+1 ), 1, CONE, ZWORK( IZ7+1 ), 1 ) 666C 667C Compute diag( Binv ). 668C 669 CALL DCOPY( N, DWORK( IW7+1 ), 1, DWORK( IW12+1 ), 1 ) 670 CALL DGEMV( 'N', N, M-1, ONE, DWORK( IW8+1 ), N, X, 1, ONE, 671 $ DWORK( IW12+1 ), 1 ) 672 DO 410 I = 1, N 673 DWORK( IW12+I ) = ONE / DWORK( IW12+I ) 674 410 CONTINUE 675C 676C Compute Binv*A. 677C 678 DO 430 J = 1, N 679 DO 420 I = 1, N 680 ZWORK( IZ11+I+(J-1)*N ) = DCMPLX( DWORK( IW12+I ) )* 681 $ ZWORK( IZ7+I+(J-1)*N ) 682 420 CONTINUE 683 430 CONTINUE 684C 685C Compute eig( Binv*A ). 686C 687 CALL ZGEES( 'N', 'N', SELECT, N, ZWORK( IZ11+1 ), N, SDIM, 688 $ ZWORK( IZ12+1 ), ZWORK, N, ZWORK( IZWRK+1 ), 689 $ LZWORK-IZWRK, DWORK( IWRK+1 ), BWORK, INFO2 ) 690 IF( INFO2.GT.0 ) THEN 691 INFO = 6 692 RETURN 693 END IF 694 LZA = INT( ZWORK( IZWRK+1 ) ) 695 LZAMAX = MAX( LZA, LZAMAX ) 696 E = DREAL( ZWORK( IZ12+1 ) ) 697 IF( N.GT.1 ) THEN 698 DO 440 I = 2, N 699 IF( DREAL( ZWORK( IZ12+I ) ).GT.E ) 700 $ E = DREAL( ZWORK( IZ12+I ) ) 701 440 CONTINUE 702 END IF 703C 704C Set tau. 705C 706 IF( MR.GT.0 ) THEN 707 SNORM = ABS( X( M ) ) 708 IF( MR.GT.1 ) THEN 709 DO 450 I = M+1, MT 710 IF( ABS( X( I ) ).GT.SNORM ) SNORM = ABS( X( I ) ) 711 450 CONTINUE 712 END IF 713 IF( SNORM.GT.FORTY ) THEN 714 TAU = C7 715 ELSE IF( SNORM.GT.EIGHT ) THEN 716 TAU = FIFTY 717 ELSE IF( SNORM.GT.FOUR ) THEN 718 TAU = TEN 719 ELSE IF( SNORM.GT.ONE ) THEN 720 TAU = FIVE 721 ELSE 722 TAU = TWO 723 END IF 724 END IF 725 IF( ITER.EQ.0 ) THEN 726 DLAMBD = E + C1 727 ELSE 728 DWORK( IW13+1 ) = E 729 CALL DCOPY( MT, X, 1, DWORK( IW13+2 ), 1 ) 730 DLAMBD = ( ONE - THETA )*DWORK( IW13+1 ) + 731 $ THETA*DWORK( IW14+1 ) 732 CALL DCOPY( MT, DWORK( IW13+2 ), 1, DWORK( IW18+1 ), 1 ) 733 CALL DCOPY( MT, DWORK( IW14+2 ), 1, DWORK( IW19+1 ), 1 ) 734 L = 0 735 460 DO 470 I = 1, MT 736 X( I ) = ( ONE - THETA / TWO**L )*DWORK( IW18+I ) + 737 $ ( THETA / TWO**L )*DWORK( IW19+I ) 738 470 CONTINUE 739C 740C Compute At(:) = A0 + AA*x. 741C 742 DO 480 I = 1, MT 743 ZWORK( IZ10+I ) = DCMPLX( X( I ) ) 744 480 CONTINUE 745 CALL ZCOPY( N*N, ZWORK( IZ5+1 ), 1, ZWORK( IZ9+1 ), 1 ) 746 CALL ZGEMV( 'N', N*N, MT, CONE, ZWORK( IZ6+1 ), N*N, 747 $ ZWORK( IZ10+1 ), 1, CONE, ZWORK( IZ9+1 ), 1 ) 748C 749C Compute diag(Bt). 750C 751 CALL DCOPY( N, DWORK( IW7+1 ), 1, DWORK( IW21+1 ), 1 ) 752 CALL DGEMV( 'N', N, M-1, ONE, DWORK( IW8+1 ), N, X, 1, ONE, 753 $ DWORK( IW21+1 ), 1 ) 754C 755C Compute W. 756C 757 DO 500 J = 1, N 758 DO 490 I = 1, N 759 IF( I.EQ.J ) THEN 760 ZWORK( IZ13+I+(I-1)*N ) = DCMPLX( THETA*BETA* 761 $ ( DWORK( IW14+1 ) - DWORK( IW13+1 ) ) /TWO - 762 $ DLAMBD*DWORK( IW21+I ) ) + 763 $ ZWORK( IZ9+I+(I-1)*N ) 764 ELSE 765 ZWORK( IZ13+I+(J-1)*N ) = ZWORK( IZ9+I+(J-1)*N ) 766 END IF 767 490 CONTINUE 768 500 CONTINUE 769C 770C Compute eig( W ). 771C 772 CALL ZGEES( 'N', 'N', SELECT, N, ZWORK( IZ13+1 ), N, SDIM, 773 $ ZWORK( IZ14+1 ), ZWORK, N, ZWORK( IZWRK+1 ), 774 $ LZWORK-IZWRK, DWORK( IWRK+1 ), BWORK, INFO2 ) 775 IF( INFO2.GT.0 ) THEN 776 INFO = 6 777 RETURN 778 END IF 779 LZA = INT( ZWORK( IZWRK+1 ) ) 780 LZAMAX = MAX( LZA, LZAMAX ) 781 EMAX = DREAL( ZWORK( IZ14+1 ) ) 782 IF( N.GT.1 ) THEN 783 DO 510 I = 2, N 784 IF( DREAL( ZWORK( IZ14+I ) ).GT.EMAX ) 785 $ EMAX = DREAL( ZWORK( IZ14+I ) ) 786 510 CONTINUE 787 END IF 788 IF( EMAX.LE.ZERO ) THEN 789 GO TO 515 790 ELSE 791 L = L + 1 792 GO TO 460 793 END IF 794 END IF 795C 796C Set y. 797C 798 515 DWORK( IW13+1 ) = DLAMBD 799 CALL DCOPY( MT, X, 1, DWORK( IW13+2 ), 1 ) 800C 801 IF( ( SVLAM - DLAMBD ).LT.TOL ) THEN 802 BOUND = SQRT( MAX( E, ZERO ) )*ZNORM 803 DO 520 I = 1, M - 1 804 X( I ) = X( I )*DWORK( IW2+I+1 )**2 805 520 CONTINUE 806C 807C Compute sqrt( x ). 808C 809 DO 530 I = 1, M-1 810 DWORK( IW20+I ) = SQRT( X( I ) ) 811 530 CONTINUE 812C 813C Compute diag( D ). 814C 815 CALL DCOPY( N, DWORK( IW7+1 ), 1, D, 1 ) 816 CALL DGEMV( 'N', N, M-1, ONE, DWORK( IW8+1 ), N, 817 $ DWORK( IW20+1 ), 1, ONE, D, 1 ) 818C 819C Compute diag( G ). 820C 821 J = 0 822 L = 0 823 DO 540 K = 1, M 824 J = J + NBLOCK( K ) 825 IF( ITYPE( K ).EQ.1 ) THEN 826 L = L + 1 827 X( M-1+L ) = X( M-1+L )*DWORK( IW2+K )**2 828 G( J ) = X( M-1+L ) 829 END IF 830 540 CONTINUE 831 CALL DSCAL( N, ZNORM, G, 1 ) 832 DWORK( 1 ) = DFLOAT( MINWRK - 5*N + LWAMAX ) 833 ZWORK( 1 ) = DCMPLX( MINZRK - 3*N + LZAMAX ) 834 RETURN 835 END IF 836 SVLAM = DLAMBD 837 DO 800 K = 1, M 838C 839C Store xD. 840C 841 CALL DCOPY( M-1, X, 1, DWORK( IW22+1 ), 1 ) 842 IF( MR.GT.0 ) THEN 843C 844C Store xG. 845C 846 CALL DCOPY( MR, X( M ), 1, DWORK( IW23+1 ), 1 ) 847 END IF 848C 849C Compute A(:) = A0 + AA*x. 850C 851 DO 550 I = 1, MT 852 ZWORK( IZ10+I ) = DCMPLX( X( I ) ) 853 550 CONTINUE 854 CALL ZCOPY( N*N, ZWORK( IZ5+1 ), 1, ZWORK( IZ7+1 ), 1 ) 855 CALL ZGEMV( 'N', N*N, MT, CONE, ZWORK( IZ6+1 ), N*N, 856 $ ZWORK( IZ10+1 ), 1, CONE, ZWORK( IZ7+1 ), 1 ) 857C 858C Compute B = B0d + BBd*xD. 859C 860 CALL DCOPY( N, DWORK( IW7+1 ), 1, DWORK( IW24+1 ), 1 ) 861 CALL DGEMV( 'N', N, M-1, ONE, DWORK( IW8+1 ), N, 862 $ DWORK( IW22+1 ), 1, ONE, DWORK( IW24+1 ), 1 ) 863C 864C Compute F. 865C 866 DO 556 J = 1, N 867 DO 555 I = 1, N 868 IF( I.EQ.J ) THEN 869 ZWORK( IZ15+I+(I-1)*N ) = DCMPLX( DLAMBD* 870 $ DWORK( IW24+I ) ) - ZWORK( IZ7+I+(I-1)*N ) 871 ELSE 872 ZWORK( IZ15+I+(J-1)*N ) = -ZWORK( IZ7+I+(J-1)*N ) 873 END IF 874 555 CONTINUE 875 556 CONTINUE 876 CALL ZLACPY( 'Full', N, N, ZWORK( IZ15+1 ), N, 877 $ ZWORK( IZ17+1 ), N ) 878C 879C Compute det( F ). 880C 881 CALL ZGEES( 'N', 'N', SELECT, N, ZWORK( IZ15+1 ), N, SDIM, 882 $ ZWORK( IZ16+1 ), ZWORK, N, ZWORK( IZWRK+1 ), 883 $ LZWORK-IZWRK, DWORK( IWRK+1 ), BWORK, INFO2 ) 884 IF( INFO2.GT.0 ) THEN 885 INFO = 6 886 RETURN 887 END IF 888 LZA = INT( ZWORK( IZWRK+1 ) ) 889 LZAMAX = MAX( LZA, LZAMAX ) 890 DETF = CONE 891 DO 560 I = 1, N 892 DETF = DETF*ZWORK( IZ16+I ) 893 560 CONTINUE 894C 895C Compute Finv. 896C 897 CALL ZGETRF( N, N, ZWORK( IZ17+1 ), N, IWORK, INFO2 ) 898 IF( INFO2.GT.0 ) THEN 899 INFO = 5 900 RETURN 901 END IF 902 CALL ZGETRI( N, ZWORK( IZ17+1 ), N, IWORK, ZWORK( IZWRK+1 ), 903 $ LDWORK-IWRK, INFO2 ) 904 LZA = INT( ZWORK( IZWRK+1 ) ) 905 LZAMAX = MAX( LZA, LZAMAX ) 906C 907C Compute phi. 908C 909 DO 570 I = 1, M-1 910 DWORK( IW25+I ) = DWORK( IW22+I ) - BETA 911 DWORK( IW25+M-1+I ) = ALPHA - DWORK( IW22+I ) 912 570 CONTINUE 913 IF( MR.GT.0 ) THEN 914 DO 580 I = 1, MR 915 DWORK( IW25+2*(M-1)+I ) = DWORK( IW23+I ) + TAU 916 DWORK( IW25+2*(M-1)+MR+I ) = TAU - DWORK( IW23+I ) 917 580 CONTINUE 918 END IF 919 PROD = ONE 920 DO 590 I = 1, 2*MT 921 PROD = PROD*DWORK( IW25+I ) 922 590 CONTINUE 923 TEMP = DREAL( DETF ) 924 IF( TEMP.LT.EPS ) TEMP = EPS 925 PHI = -LOG( TEMP ) - LOG( PROD ) 926C 927C Compute g. 928C 929 DO 610 J = 1, MT 930 DO 600 I = 1, N*N 931 ZWORK( IZ18+I+(J-1)*N*N ) = DCMPLX( DLAMBD* 932 $ DWORK( IW9+I+(J-1)*N*N ) ) - ZWORK( IZ6+I+(J-1)*N*N ) 933 600 CONTINUE 934 610 CONTINUE 935 CALL ZGEMV( 'C', N*N, MT, CONE, ZWORK( IZ18+1 ), N*N, 936 $ ZWORK( IZ17+1 ), 1, CZERO, ZWORK( IZ19+1 ), 1 ) 937 DO 620 I = 1, M-1 938 DWORK( IW26+I ) = ONE / ( DWORK( IW22+I ) - BETA ) - 939 $ ONE / ( ALPHA - DWORK( IW22+I ) ) 940 620 CONTINUE 941 IF( MR.GT.0 ) THEN 942 DO 630 I = 1, MR 943 DWORK( IW26+M-1+I ) = ONE / ( DWORK( IW23+I ) + TAU ) 944 $ -ONE / ( TAU - DWORK( IW23+I ) ) 945 630 CONTINUE 946 END IF 947 DO 640 I = 1, MT 948 DWORK( IW26+I ) = -DREAL( ZWORK( IZ19+I ) ) - 949 $ DWORK( IW26+I ) 950 640 CONTINUE 951C 952C Compute h. 953C 954 CALL DLACPY( 'Full', MT, MT, DWORK( IW11+1 ), MT, 955 $ DWORK( IW31+1 ), MT ) 956 CALL DCOPY( MT, DWORK( IW26+1 ), 1, DWORK( IW27+1 ), 1 ) 957 CALL DSYSV( 'U', MT, 1, DWORK( IW31+1 ), MT, IWORK, 958 $ DWORK( IW27+1 ), MT, DWORK( IWRK+1 ), 959 $ LDWORK-IWRK, INFO2 ) 960 IF( INFO2.GT.0 ) THEN 961 INFO = 5 962 RETURN 963 END IF 964 LWA = INT( DWORK( IWRK+1 ) ) 965 LWAMAX = MAX( LWA, LWAMAX ) 966 STSIZE = ONE 967C 968C Store hD. 969C 970 CALL DCOPY( M-1, DWORK( IW27+1 ), 1, DWORK( IW28+1 ), 1 ) 971C 972C Determine stepsize. 973C 974 L = 0 975 DO 650 I = 1, M-1 976 IF( DWORK( IW28+I ).GT.ZERO ) THEN 977 L = L + 1 978 IF( L.EQ.1 ) THEN 979 TEMP = ( DWORK( IW22+I ) - BETA ) / DWORK( IW28+I ) 980 ELSE 981 TEMP = MIN( TEMP, ( DWORK( IW22+I ) - BETA ) / 982 $ DWORK( IW28+I ) ) 983 END IF 984 END IF 985 650 CONTINUE 986 IF( L.GT.0 ) STSIZE = MIN( STSIZE, TEMP ) 987 L = 0 988 DO 660 I = 1, M-1 989 IF( DWORK( IW28+I ).LT.ZERO ) THEN 990 L = L + 1 991 IF( L.EQ.1 ) THEN 992 TEMP = ( ALPHA - DWORK( IW22+I ) ) / 993 $ ( -DWORK( IW28+I ) ) 994 ELSE 995 TEMP = MIN( TEMP, ( ALPHA - DWORK( IW22+I ) ) / 996 $ ( -DWORK( IW28+I ) ) ) 997 END IF 998 END IF 999 660 CONTINUE 1000 IF( L.GT.0 ) STSIZE = MIN( STSIZE, TEMP ) 1001 IF( MR.GT.0 ) THEN 1002C 1003C Store hG. 1004C 1005 CALL DCOPY( MR, DWORK( IW27+M ), 1, DWORK( IW29+1 ), 1 ) 1006C 1007C Determine stepsize. 1008C 1009 L = 0 1010 DO 670 I = 1, MR 1011 IF( DWORK( IW29+I ).GT.ZERO ) THEN 1012 L = L + 1 1013 IF( L.EQ.1 ) THEN 1014 TEMP = ( DWORK( IW23+I ) + TAU ) / 1015 $ DWORK( IW29+I ) 1016 ELSE 1017 TEMP = MIN( TEMP, ( DWORK( IW23+I ) + TAU ) / 1018 $ DWORK( IW29+I ) ) 1019 END IF 1020 END IF 1021 670 CONTINUE 1022 IF( L.GT.0 ) STSIZE = MIN( STSIZE, TEMP ) 1023 L = 0 1024 DO 680 I = 1, MR 1025 IF( DWORK( IW29+I ).LT.ZERO ) THEN 1026 L = L + 1 1027 IF( L.EQ.1 ) THEN 1028 TEMP = ( TAU - DWORK( IW23+I ) ) / 1029 $ ( -DWORK( IW29+I ) ) 1030 ELSE 1031 TEMP = MIN( TEMP, ( TAU - DWORK( IW23+I ) ) / 1032 $ ( -DWORK( IW29+I ) ) ) 1033 END IF 1034 END IF 1035 680 CONTINUE 1036 END IF 1037 IF( L.GT.0 ) STSIZE = MIN( STSIZE, TEMP ) 1038 STSIZE = C4*STSIZE 1039 IF( STSIZE.GE.TOL4 ) THEN 1040C 1041C Compute x_new. 1042C 1043 DO 700 I = 1, MT 1044 DWORK( IW20+I ) = X( I ) - STSIZE*DWORK( IW27+I ) 1045 700 CONTINUE 1046C 1047C Store xD. 1048C 1049 CALL DCOPY( M-1, DWORK( IW20+1 ), 1, DWORK( IW22+1 ), 1 ) 1050 IF( MR.GT.0 ) THEN 1051C 1052C Store xG. 1053C 1054 CALL DCOPY( MR, DWORK( IW20+M ), 1, DWORK( IW23+1 ), 1055 $ 1 ) 1056 END IF 1057C 1058C Compute A(:) = A0 + AA*x_new. 1059C 1060 DO 710 I = 1, MT 1061 ZWORK( IZ10+I ) = DCMPLX( DWORK( IW20+I ) ) 1062 710 CONTINUE 1063 CALL ZCOPY( N*N, ZWORK( IZ5+1 ), 1, ZWORK( IZ7+1 ), 1 ) 1064 CALL ZGEMV( 'N', N*N, MT, CONE, ZWORK( IZ6+1 ), N*N, 1065 $ ZWORK( IZ10+1 ), 1, CONE, ZWORK( IZ7+1 ), 1 ) 1066C 1067C Compute B = B0d + BBd*xD. 1068C 1069 CALL DCOPY( N, DWORK( IW7+1 ), 1, DWORK( IW24+1 ), 1 ) 1070 CALL DGEMV( 'N', N, M-1, ONE, DWORK( IW8+1 ), N, 1071 $ DWORK( IW22+1 ), 1, ONE, DWORK( IW24+1 ), 1 ) 1072C 1073C Compute lambda*diag(B) - A. 1074C 1075 DO 730 J = 1, N 1076 DO 720 I = 1, N 1077 IF( I.EQ.J ) THEN 1078 ZWORK( IZ15+I+(I-1)*N ) = DCMPLX( DLAMBD* 1079 $ DWORK( IW24+I ) ) - ZWORK( IZ7+I+(I-1)*N ) 1080 ELSE 1081 ZWORK( IZ15+I+(J-1)*N ) = 1082 $ -ZWORK( IZ7+I+(J-1)*N ) 1083 END IF 1084 720 CONTINUE 1085 730 CONTINUE 1086C 1087C Compute eig( lambda*diag(B)-A ). 1088C 1089 CALL ZGEES( 'N', 'N', SELECT, N, ZWORK( IZ15+1 ), N, 1090 $ SDIM, ZWORK( IZ16+1 ), ZWORK, N, 1091 $ ZWORK( IZWRK+1 ), LZWORK-IZWRK, 1092 $ DWORK( IWRK+1 ), BWORK, INFO2 ) 1093 IF( INFO2.GT.0 ) THEN 1094 INFO = 6 1095 RETURN 1096 END IF 1097 LZA = INT( ZWORK( IZWRK+1 ) ) 1098 LZAMAX = MAX( LZA, LZAMAX ) 1099 EMIN = DREAL( ZWORK( IZ16+1 ) ) 1100 IF( N.GT.1 ) THEN 1101 DO 740 I = 2, N 1102 IF( DREAL( ZWORK( IZ16+I ) ).LT.EMIN ) 1103 $ EMIN = DREAL( ZWORK( IZ16+I ) ) 1104 740 CONTINUE 1105 END IF 1106 DO 750 I = 1, N 1107 DWORK( IW30+I ) = DREAL( ZWORK( IZ16+I ) ) 1108 750 CONTINUE 1109 DO 760 I = 1, M-1 1110 DWORK( IW30+N+I ) = DWORK( IW22+I ) - BETA 1111 DWORK( IW30+N+M-1+I ) = ALPHA - DWORK( IW22+I ) 1112 760 CONTINUE 1113 IF( MR.GT.0 ) THEN 1114 DO 770 I = 1, MR 1115 DWORK( IW30+N+2*(M-1)+I ) = DWORK( IW23+I ) + TAU 1116 DWORK( IW30+N+2*(M-1)+MR+I ) = TAU - 1117 $ DWORK( IW23+I ) 1118 770 CONTINUE 1119 END IF 1120 PROD = ONE 1121 DO 780 I = 1, N+2*MT 1122 PROD = PROD*DWORK( IW30+I ) 1123 780 CONTINUE 1124 IF( EMIN.LE.ZERO .OR. ( -LOG( PROD ) ).GE.PHI ) THEN 1125 STSIZE = STSIZE / TEN 1126 ELSE 1127 CALL DCOPY( MT, DWORK( IW20+1 ), 1, X, 1 ) 1128 END IF 1129 END IF 1130 IF( STSIZE.LT.TOL4 ) GO TO 810 1131 800 CONTINUE 1132C 1133 810 CONTINUE 1134C 1135C Store xD. 1136C 1137 CALL DCOPY( M-1, X, 1, DWORK( IW22+1 ), 1 ) 1138 IF( MR.GT.0 ) THEN 1139C 1140C Store xG. 1141C 1142 CALL DCOPY( MR, X( M ), 1, DWORK( IW23+1 ), 1 ) 1143 END IF 1144C 1145C Compute A(:) = A0 + AA*x. 1146C 1147 DO 820 I = 1, MT 1148 ZWORK( IZ10+I ) = DCMPLX( X( I ) ) 1149 820 CONTINUE 1150 CALL ZCOPY( N*N, ZWORK( IZ5+1 ), 1, ZWORK( IZ7+1 ), 1 ) 1151 CALL ZGEMV( 'N', N*N, MT, CONE, ZWORK( IZ6+1 ), N*N, 1152 $ ZWORK( IZ10+1 ), 1, CONE, ZWORK( IZ7+1 ), 1 ) 1153C 1154C Compute diag( B ) = B0d + BBd*xD. 1155C 1156 CALL DCOPY( N, DWORK( IW7+1 ), 1, DWORK( IW24+1 ), 1 ) 1157 CALL DGEMV( 'N', N, M-1, ONE, DWORK( IW8+1 ), N, 1158 $ DWORK( IW22+1 ), 1, ONE, DWORK( IW24+1 ), 1 ) 1159C 1160C Compute F. 1161C 1162 DO 840 J = 1, N 1163 DO 830 I = 1, N 1164 IF( I.EQ.J ) THEN 1165 ZWORK( IZ15+I+(I-1)*N ) = DCMPLX( DLAMBD* 1166 $ DWORK( IW24+I ) ) - ZWORK( IZ7+I+(I-1)*N ) 1167 ELSE 1168 ZWORK( IZ15+I+(J-1)*N ) = -ZWORK( IZ7+I+(J-1)*N ) 1169 END IF 1170 830 CONTINUE 1171 840 CONTINUE 1172 CALL ZLACPY( 'Full', N, N, ZWORK( IZ15+1 ), N, 1173 $ ZWORK( IZ17+1 ), N ) 1174C 1175C Compute det( F ). 1176C 1177 CALL ZGEES( 'N', 'N', SELECT, N, ZWORK( IZ15+1 ), N, SDIM, 1178 $ ZWORK( IZ16+1 ), ZWORK, N, ZWORK( IZWRK+1 ), 1179 $ LZWORK-IZWRK, DWORK( IWRK+1 ), BWORK, INFO2 ) 1180 IF( INFO2.GT.0 ) THEN 1181 INFO = 6 1182 RETURN 1183 END IF 1184 LZA = INT( ZWORK( IZWRK+1 ) ) 1185 LZAMAX = MAX( LZA, LZAMAX ) 1186 DETF = CONE 1187 DO 850 I = 1, N 1188 DETF = DETF*ZWORK( IZ16+I ) 1189 850 CONTINUE 1190C 1191C Compute Finv. 1192C 1193 CALL ZGETRF( N, N, ZWORK( IZ17+1 ), N, IWORK, INFO2 ) 1194 IF( INFO2.GT.0 ) THEN 1195 INFO = 5 1196 RETURN 1197 END IF 1198 CALL ZGETRI( N, ZWORK( IZ17+1 ), N, IWORK, ZWORK( IZWRK+1 ), 1199 $ LDWORK-IWRK, INFO2 ) 1200 LZA = INT( ZWORK( IZWRK+1 ) ) 1201 LZAMAX = MAX( LZA, LZAMAX ) 1202C 1203C Compute the barrier function. 1204C 1205 DO 860 I = 1, M-1 1206 DWORK( IW25+I ) = DWORK( IW22+I ) - BETA 1207 DWORK( IW25+M-1+I ) = ALPHA - DWORK( IW22+I ) 1208 860 CONTINUE 1209 IF( MR.GT.0 ) THEN 1210 DO 870 I = 1, MR 1211 DWORK( IW25+2*(M-1)+I ) = DWORK( IW23+I ) + TAU 1212 DWORK( IW25+2*(M-1)+MR+I ) = TAU - DWORK( IW23+I ) 1213 870 CONTINUE 1214 END IF 1215 PROD = ONE 1216 DO 880 I = 1, 2*MT 1217 PROD = PROD*DWORK( IW25+I ) 1218 880 CONTINUE 1219 TEMP = DREAL( DETF ) 1220 IF( TEMP.LT.EPS ) TEMP = EPS 1221 PHI = -LOG( TEMP ) - LOG( PROD ) 1222C 1223C Compute the gradient of the barrier function. 1224C 1225 DO 900 J = 1, MT 1226 DO 890 I = 1, N*N 1227 ZWORK( IZ18+I+(J-1)*N*N ) = DCMPLX( DLAMBD* 1228 $ DWORK( IW9+I+(J-1)*N*N ) ) - ZWORK( IZ6+I+(J-1)*N*N ) 1229 890 CONTINUE 1230 900 CONTINUE 1231 CALL ZGEMV( 'C', N*N, MT, CONE, ZWORK( IZ18+1 ), N*N, 1232 $ ZWORK( IZ17+1 ), 1, CZERO, ZWORK( IZ19+1 ), 1 ) 1233 DO 910 I = 1, M-1 1234 DWORK( IW26+I ) = ONE / ( DWORK( IW22+I ) - BETA ) - 1235 $ ONE / ( ALPHA - DWORK( IW22+I ) ) 1236 910 CONTINUE 1237 IF( MR.GT.0 ) THEN 1238 DO 920 I = 1, MR 1239 DWORK( IW26+M-1+I ) = ONE / ( DWORK( IW23+I ) + TAU ) 1240 $ -ONE / ( TAU - DWORK( IW23+I ) ) 1241 920 CONTINUE 1242 END IF 1243 DO 925 I = 1, MT 1244 DWORK( IW26+I ) = -DREAL( ZWORK( IZ19+I ) ) - 1245 $ DWORK( IW26+I ) 1246 925 CONTINUE 1247C 1248C Compute the Hessian of the barrier function. 1249C 1250 CALL ZGEMM( 'N', 'N', N, N*MT, N, CONE, ZWORK( IZ17+1 ), N, 1251 $ ZWORK( IZ18+1 ), N, CZERO, ZWORK( IZ20+1 ), N ) 1252 1253 CALL DLASET( 'Full', MT, MT, ZERO, ZERO, DWORK( IW11+1 ), 1254 $ MT ) 1255 DO 960 K = 1, MT 1256 CALL ZCOPY( N*N, ZWORK( IZ20+1+(K-1)*N*N ), 1, 1257 $ ZWORK( IZ22+1 ), 1 ) 1258 DO 940 J = 1, N 1259 DO 930 I = 1, N 1260 ZWORK( IZ23+I+(J-1)*N ) = 1261 $ DCONJG( ZWORK( IZ22+J+(I-1)*N ) ) 1262 930 CONTINUE 1263 940 CONTINUE 1264 CALL ZGEMV( 'C', N*N, K, CONE, ZWORK( IZ20+1 ), N*N, 1265 $ ZWORK( IZ23+1 ), 1, CZERO, ZWORK( IZ24+1 ), 1266 $ 1 ) 1267 DO 950 J = 1, K 1268 DWORK( IW11+K+(J-1)*MT ) = 1269 $ DREAL( DCONJG( ZWORK( IZ24+J ) ) ) 1270 950 CONTINUE 1271 960 CONTINUE 1272 DO 970 I = 1, M-1 1273 DWORK( IW10+I ) = ONE / ( DWORK( IW22+I ) - BETA )**2 + 1274 $ ONE / ( ALPHA - DWORK( IW22+I ) )**2 1275 970 CONTINUE 1276 IF( MR.GT.0 ) THEN 1277 DO 980 I = 1, MR 1278 DWORK( IW10+M-1+I ) = 1279 $ ONE / ( DWORK( IW23+I ) + TAU )**2 + 1280 $ ONE / ( TAU - DWORK( IW23+I ) )**2 1281 980 CONTINUE 1282 END IF 1283 DO 990 I = 1, MT 1284 DWORK( IW11+I+(I-1)*MT ) = DWORK( IW11+I+(I-1)*MT ) + 1285 $ DWORK( IW10+I ) 1286 990 CONTINUE 1287 DO 1100 J = 1, MT 1288 DO 1000 I = 1, J 1289 IF( I.NE.J ) THEN 1290 T1 = DWORK( IW11+I+(J-1)*MT ) 1291 T2 = DWORK( IW11+J+(I-1)*MT ) 1292 DWORK( IW11+I+(J-1)*MT ) = T1 + T2 1293 DWORK( IW11+J+(I-1)*MT ) = T1 + T2 1294 END IF 1295 1000 CONTINUE 1296 1100 CONTINUE 1297C 1298C Compute norm( H ). 1299C 1300 1110 HNORM = DLANGE( 'F', MT, MT, DWORK( IW11+1 ), MT, DWORK ) 1301C 1302C Compute rcond( H ). 1303C 1304 CALL DLACPY( 'Full', MT, MT, DWORK( IW11+1 ), MT, 1305 $ DWORK( IW31+1 ), MT ) 1306 HNORM1 = DLANGE( '1', MT, MT, DWORK( IW31+1 ), MT, DWORK ) 1307 CALL DSYTRF( 'U', MT, DWORK( IW31+1 ), MT, IWORK, 1308 $ DWORK( IWRK+1 ), LDWORK-IWRK, INFO2 ) 1309 IF( INFO2.GT.0 ) THEN 1310 INFO = 5 1311 RETURN 1312 END IF 1313 LWA = INT( DWORK( IWRK+1 ) ) 1314 LWAMAX = MAX( LWA, LWAMAX ) 1315 CALL DSYCON( 'U', MT, DWORK( IW31+1 ), MT, IWORK, HNORM1, 1316 $ RCOND, DWORK( IWRK+1 ), IWORK( MT+1 ), INFO2 ) 1317 IF( RCOND.LT.TOL3 ) THEN 1318 DO 1120 I = 1, MT 1319 DWORK( IW11+I+(I-1)*MT ) = DWORK( IW11+I+(I-1)*MT ) + 1320 $ HNORM*REGPAR 1321 1120 CONTINUE 1322 GO TO 1110 1323 END IF 1324C 1325C Compute the tangent line to path of center. 1326C 1327 CALL DCOPY( MT, DWORK( IW26+1 ), 1, DWORK( IW27+1 ), 1 ) 1328 CALL DSYTRS( 'U', MT, 1, DWORK( IW31+1 ), MT, IWORK, 1329 $ DWORK( IW27+1 ), MT, INFO2 ) 1330C 1331C Check if x-h satisfies the Goldstein test. 1332C 1333 GTEST = .FALSE. 1334 DO 1130 I = 1, MT 1335 DWORK( IW20+I ) = X( I ) - DWORK( IW27+I ) 1336 1130 CONTINUE 1337C 1338C Store xD. 1339C 1340 CALL DCOPY( M-1, DWORK( IW20+1 ), 1, DWORK( IW22+1 ), 1 ) 1341 IF( MR.GT.0 ) THEN 1342C 1343C Store xG. 1344C 1345 CALL DCOPY( MR, DWORK( IW20+M ), 1, DWORK( IW23+1 ), 1 ) 1346 END IF 1347C 1348C Compute A(:) = A0 + AA*x_new. 1349C 1350 DO 1140 I = 1, MT 1351 ZWORK( IZ10+I ) = DCMPLX( DWORK( IW20+I ) ) 1352 1140 CONTINUE 1353 CALL ZCOPY( N*N, ZWORK( IZ5+1 ), 1, ZWORK( IZ7+1 ), 1 ) 1354 CALL ZGEMV( 'N', N*N, MT, CONE, ZWORK( IZ6+1 ), N*N, 1355 $ ZWORK( IZ10+1 ), 1, CONE, ZWORK( IZ7+1 ), 1 ) 1356C 1357C Compute diag( B ) = B0d + BBd*xD. 1358C 1359 CALL DCOPY( N, DWORK( IW7+1 ), 1, DWORK( IW24+1 ), 1 ) 1360 CALL DGEMV( 'N', N, M-1, ONE, DWORK( IW8+1 ), N, 1361 $ DWORK( IW22+1 ), 1, ONE, DWORK( IW24+1 ), 1 ) 1362C 1363C Compute lambda*diag(B) - A. 1364C 1365 DO 1160 J = 1, N 1366 DO 1150 I = 1, N 1367 IF( I.EQ.J ) THEN 1368 ZWORK( IZ15+I+(I-1)*N ) = DCMPLX( DLAMBD* 1369 $ DWORK( IW24+I ) ) - ZWORK( IZ7+I+(I-1)*N ) 1370 ELSE 1371 ZWORK( IZ15+I+(J-1)*N ) = -ZWORK( IZ7+I+(J-1)*N ) 1372 END IF 1373 1150 CONTINUE 1374 1160 CONTINUE 1375C 1376C Compute eig( lambda*diag(B)-A ). 1377C 1378 CALL ZGEES( 'N', 'N', SELECT, N, ZWORK( IZ15+1 ), N, SDIM, 1379 $ ZWORK( IZ16+1 ), ZWORK, N, ZWORK( IZWRK+1 ), 1380 $ LZWORK-IZWRK, DWORK( IWRK+1 ), BWORK, INFO2 ) 1381 IF( INFO2.GT.0 ) THEN 1382 INFO = 6 1383 RETURN 1384 END IF 1385 LZA = INT( ZWORK( IZWRK+1 ) ) 1386 LZAMAX = MAX( LZA, LZAMAX ) 1387 DO 1190 I = 1, N 1388 DWORK( IW30+I ) = DREAL( ZWORK( IZ16+I ) ) 1389 1190 CONTINUE 1390 DO 1200 I = 1, M-1 1391 DWORK( IW30+N+I ) = DWORK( IW22+I ) - BETA 1392 DWORK( IW30+N+M-1+I ) = ALPHA - DWORK( IW22+I ) 1393 1200 CONTINUE 1394 IF( MR.GT.0 ) THEN 1395 DO 1210 I = 1, MR 1396 DWORK( IW30+N+2*(M-1)+I ) = DWORK( IW23+I ) + TAU 1397 DWORK( IW30+N+2*(M-1)+MR+I ) = TAU - DWORK( IW23+I ) 1398 1210 CONTINUE 1399 END IF 1400 EMIN = DWORK( IW30+1 ) 1401 DO 1220 I = 1, N+2*MT 1402 IF( DWORK( IW30+I ).LT.EMIN ) EMIN = DWORK( IW30+I ) 1403 1220 CONTINUE 1404 IF( EMIN.LE.ZERO ) THEN 1405 GTEST = .FALSE. 1406 ELSE 1407 PP = DDOT( MT, DWORK( IW26+1 ), 1, DWORK( IW27+1 ), 1 ) 1408 PROD = ONE 1409 DO 1230 I = 1, N+2*MT 1410 PROD = PROD*DWORK( IW30+I ) 1411 1230 CONTINUE 1412 T1 = -LOG( PROD ) 1413 T2 = PHI - C2*PP 1414 T3 = PHI - C4*PP 1415 IF( T1.GE.T3 .AND. T1.LT.T2 ) GTEST = .TRUE. 1416 END IF 1417C 1418C Use x-h if Goldstein test is satisfied. Otherwise use 1419C Nesterov-Nemirovsky's stepsize length. 1420C 1421 PP = DDOT( MT, DWORK( IW26+1 ), 1, DWORK( IW27+1 ), 1 ) 1422 DELTA = SQRT( PP ) 1423 IF( GTEST .OR. DELTA.LE.C3 ) THEN 1424 DO 1240 I = 1, MT 1425 X( I ) = X( I ) - DWORK( IW27+I ) 1426 1240 CONTINUE 1427 ELSE 1428 DO 1250 I = 1, MT 1429 X( I ) = X( I ) - DWORK( IW27+I ) / ( ONE + DELTA ) 1430 1250 CONTINUE 1431 END IF 1432C 1433C Analytic center is found if delta is sufficiently small. 1434C 1435 IF( DELTA.LT.TOL5 ) GO TO 1260 1436 GO TO 810 1437C 1438C Set yf. 1439C 1440 1260 DWORK( IW14+1 ) = DLAMBD 1441 CALL DCOPY( MT, X, 1, DWORK( IW14+2 ), 1 ) 1442C 1443C Set yw. 1444C 1445 CALL DCOPY( MT+1, DWORK( IW14+1 ), 1, DWORK( IW15+1 ), 1 ) 1446C 1447C Compute Fb. 1448C 1449 DO 1280 J = 1, N 1450 DO 1270 I = 1, N 1451 ZWORK( IZ21+I+(J-1)*N ) = DCMPLX( DWORK( IW24+I ) )* 1452 $ DCONJG( ZWORK( IZ17+J+(I-1)*N ) ) 1453 1270 CONTINUE 1454 1280 CONTINUE 1455 CALL ZGEMV( 'C', N*N, MT, CONE, ZWORK( IZ20+1 ), N*N, 1456 $ ZWORK( IZ21+1 ), 1, CZERO, ZWORK( IZ24+1 ), 1 ) 1457 DO 1300 I = 1, MT 1458 DWORK( IW32+I ) = DREAL( ZWORK( IZ24+I ) ) 1459 1300 CONTINUE 1460C 1461C Compute h1. 1462C 1463 CALL DLACPY( 'Full', MT, MT, DWORK( IW11+1 ), MT, 1464 $ DWORK( IW31+1 ), MT ) 1465 CALL DSYSV( 'U', MT, 1, DWORK( IW31+1 ), MT, IWORK, 1466 $ DWORK( IW32+1 ), MT, DWORK( IWRK+1 ), 1467 $ LDWORK-IWRK, INFO2 ) 1468 IF( INFO2.GT.0 ) THEN 1469 INFO = 5 1470 RETURN 1471 END IF 1472 LWA = INT( DWORK( IWRK+1 ) ) 1473 LWAMAX = MAX( LWA, LWAMAX ) 1474C 1475C Compute hn. 1476C 1477 HN = DLANGE( 'F', MT, 1, DWORK( IW32+1 ), MT, DWORK ) 1478C 1479C Compute y. 1480C 1481 DWORK( IW13+1 ) = DLAMBD - C / HN 1482 DO 1310 I = 1, MT 1483 DWORK( IW13+1+I ) = X( I ) + C*DWORK( IW32+I ) / HN 1484 1310 CONTINUE 1485C 1486C Store xD. 1487C 1488 CALL DCOPY( M-1, DWORK( IW13+2 ), 1, DWORK( IW22+1 ), 1 ) 1489 IF( MR.GT.0 ) THEN 1490C 1491C Store xG. 1492C 1493 CALL DCOPY( MR, DWORK( IW13+M+1 ), 1, DWORK( IW23+1 ), 1 ) 1494 END IF 1495C 1496C Compute A(:) = A0 + AA*y(2:mt+1). 1497C 1498 DO 1320 I = 1, MT 1499 ZWORK( IZ10+I ) = DCMPLX( DWORK( IW13+1+I ) ) 1500 1320 CONTINUE 1501 CALL ZCOPY( N*N, ZWORK( IZ5+1 ), 1, ZWORK( IZ7+1 ), 1 ) 1502 CALL ZGEMV( 'N', N*N, MT, CONE, ZWORK( IZ6+1 ), N*N, 1503 $ ZWORK( IZ10+1 ), 1, CONE, ZWORK( IZ7+1 ), 1 ) 1504C 1505C Compute B = B0d + BBd*xD. 1506C 1507 CALL DCOPY( N, DWORK( IW7+1 ), 1, DWORK( IW24+1 ), 1 ) 1508 CALL DGEMV( 'N', N, M-1, ONE, DWORK( IW8+1 ), N, 1509 $ DWORK( IW22+1 ), 1, ONE, DWORK( IW24+1 ), 1 ) 1510C 1511C Compute y(1)*diag(B) - A. 1512C 1513 DO 1340 J = 1, N 1514 DO 1330 I = 1, N 1515 IF( I.EQ.J ) THEN 1516 ZWORK( IZ15+I+(I-1)*N ) = DCMPLX( DWORK( IW13+1 )* 1517 $ DWORK( IW24+I ) ) - ZWORK( IZ7+I+(I-1)*N ) 1518 ELSE 1519 ZWORK( IZ15+I+(J-1)*N ) = -ZWORK( IZ7+I+(J-1)*N ) 1520 END IF 1521 1330 CONTINUE 1522 1340 CONTINUE 1523C 1524C Compute eig( y(1)*diag(B)-A ). 1525C 1526 CALL ZGEES( 'N', 'N', SELECT, N, ZWORK( IZ15+1 ), N, SDIM, 1527 $ ZWORK( IZ16+1 ), ZWORK, N, ZWORK( IZWRK+1 ), 1528 $ LZWORK-IZWRK, DWORK( IWRK+1 ), BWORK, INFO2 ) 1529 IF( INFO2.GT.0 ) THEN 1530 INFO = 6 1531 RETURN 1532 END IF 1533 LZA = INT( ZWORK( IZWRK+1 ) ) 1534 LZAMAX = MAX( LZA, LZAMAX ) 1535 EMIN = DREAL( ZWORK( IZ16+1 ) ) 1536 IF( N.GT.1 ) THEN 1537 DO 1350 I = 2, N 1538 IF( DREAL( ZWORK( IZ16+I ) ).LT.EMIN ) 1539 $ EMIN = DREAL( ZWORK( IZ16+I ) ) 1540 1350 CONTINUE 1541 END IF 1542 POS = .TRUE. 1543 DO 1360 I = 1, M-1 1544 DWORK( IW25+I ) = DWORK( IW22+I ) - BETA 1545 DWORK( IW25+M-1+I ) = ALPHA - DWORK( IW22+I ) 1546 1360 CONTINUE 1547 IF( MR.GT.0 ) THEN 1548 DO 1370 I = 1, MR 1549 DWORK( IW25+2*(M-1)+I ) = DWORK( IW23+I ) + TAU 1550 DWORK( IW25+2*(M-1)+MR+I ) = TAU - DWORK( IW23+I ) 1551 1370 CONTINUE 1552 END IF 1553 TEMP = DWORK( IW25+1 ) 1554 DO 1380 I = 2, 2*MT 1555 IF( DWORK( IW25+I ).LT.TEMP ) TEMP = DWORK( IW25+I ) 1556 1380 CONTINUE 1557 IF( TEMP.LE.ZERO .OR. EMIN.LE.ZERO ) POS = .FALSE. 1558 1390 IF( POS ) THEN 1559C 1560C Set y2 = y. 1561C 1562 CALL DCOPY( MT+1, DWORK( IW13+1 ), 1, DWORK( IW17+1 ), 1 ) 1563C 1564C Compute y = y + 1.5*( y - yw ). 1565C 1566 DO 1400 I = 1, MT+1 1567 DWORK( IW13+I ) = DWORK( IW13+I ) + 1568 $ C5*( DWORK( IW13+I ) - DWORK( IW15+I ) ) 1569 1400 CONTINUE 1570C 1571C Store xD. 1572C 1573 CALL DCOPY( M-1, DWORK( IW13+2 ), 1, DWORK( IW22+1 ), 1 ) 1574 IF( MR.GT.0 ) THEN 1575C 1576C Store xG. 1577C 1578 CALL DCOPY( MR, DWORK( IW13+M+1 ), 1, 1579 $ DWORK( IW23+1 ), 1 ) 1580 END IF 1581C 1582C Compute A(:) = A0 + AA*y(2:mt+1). 1583C 1584 DO 1420 I = 1, MT 1585 ZWORK( IZ10+I ) = DCMPLX( DWORK( IW13+1+I ) ) 1586 1420 CONTINUE 1587 CALL ZCOPY( N*N, ZWORK( IZ5+1 ), 1, ZWORK( IZ7+1 ), 1 ) 1588 CALL ZGEMV( 'N', N*N, MT, CONE, ZWORK( IZ6+1 ), N*N, 1589 $ ZWORK( IZ10+1 ), 1, CONE, ZWORK( IZ7+1 ), 1 ) 1590C 1591C Compute diag( B ) = B0d + BBd*xD. 1592C 1593 CALL DCOPY( N, DWORK( IW7+1 ), 1, DWORK( IW24+1 ), 1 ) 1594 CALL DGEMV( 'N', N, M-1, ONE, DWORK( IW8+1 ), N, 1595 $ DWORK( IW22+1 ), 1, ONE, DWORK( IW24+1 ), 1 ) 1596C 1597C Set yw = y2. 1598C 1599 CALL DCOPY( MT+1, DWORK( IW17+1 ), 1, DWORK( IW15+1 ), 1 ) 1600C 1601C Compute y(1)*diag(B) - A. 1602C 1603 DO 1440 J = 1, N 1604 DO 1430 I = 1, N 1605 IF( I.EQ.J ) THEN 1606 ZWORK( IZ15+I+(I-1)*N ) = DCMPLX( DWORK( IW13+1 )* 1607 $ DWORK( IW24+I ) ) - ZWORK( IZ7+I+(I-1)*N ) 1608 ELSE 1609 ZWORK( IZ15+I+(J-1)*N ) = -ZWORK( IZ7+I+(J-1)*N ) 1610 END IF 1611 1430 CONTINUE 1612 1440 CONTINUE 1613C 1614C Compute eig( y(1)*diag(B)-A ). 1615C 1616 CALL ZGEES( 'N', 'N', SELECT, N, ZWORK( IZ15+1 ), N, SDIM, 1617 $ ZWORK( IZ16+1 ), ZWORK, N, ZWORK( IZWRK+1 ), 1618 $ LZWORK-IZWRK, DWORK( IWRK+1 ), BWORK, INFO2 ) 1619 IF( INFO2.GT.0 ) THEN 1620 INFO = 6 1621 RETURN 1622 END IF 1623 LZA = INT( ZWORK( IZWRK+1 ) ) 1624 LZAMAX = MAX( LZA, LZAMAX ) 1625 EMIN = DREAL( ZWORK( IZ16+1 ) ) 1626 IF( N.GT.1 ) THEN 1627 DO 1450 I = 2, N 1628 IF( DREAL( ZWORK( IZ16+I ) ).LT.EMIN ) 1629 $ EMIN = DREAL( ZWORK( IZ16+I ) ) 1630 1450 CONTINUE 1631 END IF 1632 POS = .TRUE. 1633 DO 1460 I = 1, M-1 1634 DWORK( IW25+I ) = DWORK( IW22+I ) - BETA 1635 DWORK( IW25+M-1+I ) = ALPHA - DWORK( IW22+I ) 1636 1460 CONTINUE 1637 IF( MR.GT.0 ) THEN 1638 DO 1470 I = 1, MR 1639 DWORK( IW25+2*(M-1)+I ) = DWORK( IW23+I ) + TAU 1640 DWORK( IW25+2*(M-1)+MR+I ) = TAU - DWORK( IW23+I ) 1641 1470 CONTINUE 1642 END IF 1643 TEMP = DWORK( IW25+1 ) 1644 DO 1480 I = 2, 2*MT 1645 IF( DWORK( IW25+I ).LT.TEMP ) TEMP = DWORK( IW25+I ) 1646 1480 CONTINUE 1647 IF( TEMP.LE.ZERO .OR. EMIN.LE.ZERO ) POS = .FALSE. 1648 GO TO 1390 1649 END IF 1650 1490 CONTINUE 1651C 1652C Set y1 = ( y + yw ) / 2. 1653C 1654 DO 1500 I = 1, MT+1 1655 DWORK( IW16+I ) = ( DWORK( IW13+I ) + DWORK( IW15+I ) ) 1656 $ / TWO 1657 1500 CONTINUE 1658C 1659C Store xD. 1660C 1661 CALL DCOPY( M-1, DWORK( IW16+2 ), 1, DWORK( IW22+1 ), 1 ) 1662 IF( MR.GT.0 ) THEN 1663C 1664C Store xG. 1665C 1666 CALL DCOPY( MR, DWORK( IW16+M+1 ), 1, DWORK( IW23+1 ), 1 ) 1667 END IF 1668C 1669C Compute A(:) = A0 + AA*y1(2:mt+1). 1670C 1671 DO 1510 I = 1, MT 1672 ZWORK( IZ10+I ) = DCMPLX( DWORK( IW16+1+I ) ) 1673 1510 CONTINUE 1674 CALL ZCOPY( N*N, ZWORK( IZ5+1 ), 1, ZWORK( IZ7+1 ), 1 ) 1675 CALL ZGEMV( 'N', N*N, MT, CONE, ZWORK( IZ6+1 ), N*N, 1676 $ ZWORK( IZ10+1 ), 1, CONE, ZWORK( IZ7+1 ), 1 ) 1677C 1678C Compute diag( B ) = B0d + BBd*xD. 1679C 1680 CALL DCOPY( N, DWORK( IW7+1 ), 1, DWORK( IW24+1 ), 1 ) 1681 CALL DGEMV( 'N', N, M-1, ONE, DWORK( IW8+1 ), N, 1682 $ DWORK( IW22+1 ), 1, ONE, DWORK( IW24+1 ), 1 ) 1683C 1684C Compute y1(1)*diag(B) - A. 1685C 1686 DO 1530 J = 1, N 1687 DO 1520 I = 1, N 1688 IF( I.EQ.J ) THEN 1689 ZWORK( IZ15+I+(I-1)*N ) = DCMPLX( DWORK( IW16+1 )* 1690 $ DWORK( IW24+I ) ) - ZWORK( IZ7+I+(I-1)*N ) 1691 ELSE 1692 ZWORK( IZ15+I+(J-1)*N ) = -ZWORK( IZ7+I+(J-1)*N ) 1693 END IF 1694 1520 CONTINUE 1695 1530 CONTINUE 1696C 1697C Compute eig( y1(1)*diag(B)-A ). 1698C 1699 CALL ZGEES( 'N', 'N', SELECT, N, ZWORK( IZ15+1 ), N, SDIM, 1700 $ ZWORK( IZ16+1 ), ZWORK, N, ZWORK( IZWRK+1 ), 1701 $ LZWORK-IZWRK, DWORK( IWRK+1 ), BWORK, INFO2 ) 1702 IF( INFO2.GT.0 ) THEN 1703 INFO = 6 1704 RETURN 1705 END IF 1706 LZA = INT( ZWORK( IZWRK+1 ) ) 1707 LZAMAX = MAX( LZA, LZAMAX ) 1708 EMIN = DREAL( ZWORK( IZ16+1 ) ) 1709 IF( N.GT.1 ) THEN 1710 DO 1540 I = 2, N 1711 IF( DREAL( ZWORK( IZ16+I ) ).LT.EMIN ) 1712 $ EMIN = DREAL( ZWORK( IZ16+I ) ) 1713 1540 CONTINUE 1714 END IF 1715 POS = .TRUE. 1716 DO 1550 I = 1, M-1 1717 DWORK( IW25+I ) = DWORK( IW22+I ) - BETA 1718 DWORK( IW25+M-1+I ) = ALPHA - DWORK( IW22+I ) 1719 1550 CONTINUE 1720 IF( MR.GT.0 ) THEN 1721 DO 1560 I = 1, MR 1722 DWORK( IW25+2*(M-1)+I ) = DWORK( IW23+I ) + TAU 1723 DWORK( IW25+2*(M-1)+MR+I ) = TAU - DWORK( IW23+I ) 1724 1560 CONTINUE 1725 END IF 1726 TEMP = DWORK( IW25+1 ) 1727 DO 1570 I = 2, 2*MT 1728 IF( DWORK( IW25+I ).LT.TEMP ) TEMP = DWORK( IW25+I ) 1729 1570 CONTINUE 1730 IF( TEMP.LE.ZERO .OR. EMIN.LE.ZERO ) POS = .FALSE. 1731 IF( POS ) THEN 1732C 1733C Set yw = y1. 1734C 1735 CALL DCOPY( MT+1, DWORK( IW16+1 ), 1, DWORK( IW15+1 ), 1 ) 1736 ELSE 1737C 1738C Set y = y1. 1739C 1740 CALL DCOPY( MT+1, DWORK( IW16+1 ), 1, DWORK( IW13+1 ), 1 ) 1741 END IF 1742 DO 1580 I = 1, MT+1 1743 DWORK( IW33+I ) = DWORK( IW13+I ) - DWORK( IW15+I ) 1744 1580 CONTINUE 1745 YNORM1 = DLANGE( 'F', MT+1, 1, DWORK( IW33+1 ), MT+1, DWORK ) 1746 DO 1590 I = 1, MT+1 1747 DWORK( IW33+I ) = DWORK( IW13+I ) - DWORK( IW14+I ) 1748 1590 CONTINUE 1749 YNORM2 = DLANGE( 'F', MT+1, 1, DWORK( IW33+1 ), MT+1, DWORK ) 1750 IF( YNORM1.LT.YNORM2*THETA ) GO TO 1600 1751 GO TO 1490 1752C 1753C Compute c. 1754C 1755 1600 DO 1610 I = 1, MT+1 1756 DWORK( IW33+I ) = DWORK( IW15+I ) - DWORK( IW14+I ) 1757 1610 CONTINUE 1758 C = DLANGE( 'F', MT+1, 1, DWORK( IW33+1 ), MT+1, DWORK ) 1759C 1760C Set x = yw(2:mt+1). 1761C 1762 CALL DCOPY( MT, DWORK( IW15+2 ), 1, X, 1 ) 1763 GO TO 390 1764C 1765C *** Last line of AB13MD *** 1766 END 1767