1 SUBROUTINE AG08BY( FIRST, N, M, P, SVLMAX, ABCD, LDABCD, E, LDE, 2 $ NR, PR, NINFZ, DINFZ, NKRONL, INFZ, KRONL, 3 $ TOL, IWORK, DWORK, LDWORK, INFO ) 4C 5C WARNING 6C 7C This file alters the SLICOT routine AG08BY. It is intended 8C for use from the Octave Control Systems Package and modifies 9C the original SLICOT implementation. This file itself is *NOT* 10C part of SLICOT and the authors from NICONET e.V. are *NOT* 11C responsible for it. See file DESCRIPTION for the current 12C maintainer of the Octave control package. 13C 14C In altered implementation the deprecated LAPACK routine DLATZM 15C is replaced by ODLTZM. 16C 17C 18C SLICOT RELEASE 5.0. 19C 20C Copyright (c) 2002-2010 NICONET e.V. 21C 22C This program is free software: you can redistribute it and/or 23C modify it under the terms of the GNU General Public License as 24C published by the Free Software Foundation, either version 2 of 25C the License, or (at your option) any later version. 26C 27C This program is distributed in the hope that it will be useful, 28C but WITHOUT ANY WARRANTY; without even the implied warranty of 29C MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 30C GNU General Public License for more details. 31C 32C You should have received a copy of the GNU General Public License 33C along with this program. If not, see 34C <http://www.gnu.org/licenses/>. 35C 36C PURPOSE 37C 38C To extract from the (N+P)-by-(M+N) descriptor system pencil 39C 40C S(lambda) = ( B A - lambda*E ) 41C ( D C ) 42C 43C with E nonsingular and upper triangular a 44C (NR+PR)-by-(M+NR) "reduced" descriptor system pencil 45C 46C ( Br Ar-lambda*Er ) 47C Sr(lambda) = ( ) 48C ( Dr Cr ) 49C 50C having the same finite Smith zeros as the pencil 51C S(lambda) but with Dr, a PR-by-M full row rank 52C left upper trapezoidal matrix, and Er, an NR-by-NR 53C upper triangular nonsingular matrix. 54C 55C ARGUMENTS 56C 57C Mode Parameters 58C 59C FIRST LOGICAL 60C Specifies if AG08BY is called first time or it is called 61C for an already reduced system, with D full column rank 62C with the last M rows in upper triangular form: 63C FIRST = .TRUE., first time called; 64C FIRST = .FALSE., not first time called. 65C 66C Input/Output Parameters 67C 68C N (input) INTEGER 69C The number of rows of matrix B, the number of columns of 70C matrix C and the order of square matrices A and E. 71C N >= 0. 72C 73C M (input) INTEGER 74C The number of columns of matrices B and D. M >= 0. 75C M <= P if FIRST = .FALSE. . 76C 77C P (input) INTEGER 78C The number of rows of matrices C and D. P >= 0. 79C 80C SVLMAX (input) DOUBLE PRECISION 81C During each reduction step, the rank-revealing QR 82C factorization of a matrix stops when the estimated minimum 83C singular value is smaller than TOL * MAX(SVLMAX,EMSV), 84C where EMSV is the estimated maximum singular value. 85C SVLMAX >= 0. 86C 87C ABCD (input/output) DOUBLE PRECISION array, dimension 88C (LDABCD,M+N) 89C On entry, the leading (N+P)-by-(M+N) part of this array 90C must contain the compound matrix 91C ( B A ) , 92C ( D C ) 93C where A is an N-by-N matrix, B is an N-by-M matrix, 94C C is a P-by-N matrix and D is a P-by-M matrix. 95C If FIRST = .FALSE., then D must be a full column 96C rank matrix with the last M rows in upper triangular form. 97C On exit, the leading (NR+PR)-by-(M+NR) part of ABCD 98C contains the reduced compound matrix 99C ( Br Ar ) , 100C ( Dr Cr ) 101C where Ar is an NR-by-NR matrix, Br is an NR-by-M matrix, 102C Cr is a PR-by-NR matrix, Dr is a PR-by-M full row rank 103C left upper trapezoidal matrix with the first PR columns 104C in upper triangular form. 105C 106C LDABCD INTEGER 107C The leading dimension of array ABCD. 108C LDABCD >= MAX(1,N+P). 109C 110C E (input/output) DOUBLE PRECISION array, dimension (LDE,N) 111C On entry, the leading N-by-N part of this array must 112C contain the upper triangular nonsingular matrix E. 113C On exit, the leading NR-by-NR part contains the reduced 114C upper triangular nonsingular matrix Er. 115C 116C LDE INTEGER 117C The leading dimension of array E. LDE >= MAX(1,N). 118C 119C NR (output) INTEGER 120C The order of the reduced matrices Ar and Er; also the 121C number of rows of the reduced matrix Br and the number 122C of columns of the reduced matrix Cr. 123C If Dr is invertible, NR is also the number of finite 124C Smith zeros. 125C 126C PR (output) INTEGER 127C The rank of the resulting matrix Dr; also the number of 128C rows of reduced matrices Cr and Dr. 129C 130C NINFZ (output) INTEGER 131C Number of infinite zeros. NINFZ = 0 if FIRST = .FALSE. . 132C 133C DINFZ (output) INTEGER 134C The maximal multiplicity of infinite zeros. 135C DINFZ = 0 if FIRST = .FALSE. . 136C 137C NKRONL (output) INTEGER 138C The maximal dimension of left elementary Kronecker blocks. 139C 140C INFZ (output) INTEGER array, dimension (N) 141C INFZ(i) contains the number of infinite zeros of 142C degree i, where i = 1,2,...,DINFZ. 143C INFZ is not referenced if FIRST = .FALSE. . 144C 145C KRONL (output) INTEGER array, dimension (N+1) 146C KRONL(i) contains the number of left elementary Kronecker 147C blocks of dimension i-by-(i-1), where i = 1,2,...,NKRONL. 148C 149C Tolerances 150C 151C TOL DOUBLE PRECISION 152C A tolerance used in rank decisions to determine the 153C effective rank, which is defined as the order of the 154C largest leading (or trailing) triangular submatrix in the 155C QR (or RQ) factorization with column (or row) pivoting 156C whose estimated condition number is less than 1/TOL. 157C If the user sets TOL <= 0, then an implicitly computed, 158C default tolerance TOLDEF = (N+P)*(N+M)*EPS, is used 159C instead, where EPS is the machine precision 160C (see LAPACK Library routine DLAMCH). 161C NOTE that when SVLMAX > 0, the estimated ranks could be 162C less than those defined above (see SVLMAX). TOL <= 1. 163C 164C Workspace 165C 166C IWORK INTEGER array, dimension (M) 167C If FIRST = .FALSE., IWORK is not referenced. 168C 169C DWORK DOUBLE PRECISION array, dimension (LDWORK) 170C On exit, if INFO = 0, DWORK(1) returns the optimal value 171C of LDWORK. 172C 173C LDWORK INTEGER 174C The length of the array DWORK. 175C LDWORK >= 1, if P = 0; otherwise 176C LDWORK >= MAX( 1, N+M-1, MIN(P,M) + MAX(3*M-1,N), 5*P ), 177C if FIRST = .TRUE.; 178C LDWORK >= MAX( 1, N+M-1, 5*P ), if FIRST = .FALSE. . 179C The second term is not needed if M = 0. 180C For optimum performance LDWORK should be larger. 181C 182C If LDWORK = -1, then a workspace query is assumed; 183C the routine only calculates the optimal size of the 184C DWORK array, returns this value as the first entry of 185C the DWORK array, and no error message related to LDWORK 186C is issued by XERBLA. 187C 188C Error Indicator 189C 190C INFO INTEGER 191C = 0: successful exit; 192C < 0: if INFO = -i, the i-th argument had an illegal 193C value. 194C 195C METHOD 196C 197C The subroutine is based on the reduction algorithm of [1]. 198C 199C REFERENCES 200C 201C [1] P. Misra, P. Van Dooren and A. Varga. 202C Computation of structural invariants of generalized 203C state-space systems. 204C Automatica, 30, pp. 1921-1936, 1994. 205C 206C NUMERICAL ASPECTS 207C 208C The algorithm is numerically backward stable and requires 209C 0( (P+N)*(M+N)*N ) floating point operations. 210C 211C FURTHER COMMENTS 212C 213C The number of infinite zeros is computed as 214C 215C DINFZ 216C NINFZ = Sum (INFZ(i)*i) . 217C i=1 218C Note that each infinite zero of multiplicity k corresponds to 219C an infinite eigenvalue of multiplicity k+1. 220C The multiplicities of the infinite eigenvalues can be determined 221C from PR, DINFZ and INFZ(i), i = 1, ..., DINFZ, as follows: 222C 223C DINFZ 224C - there are PR - Sum (INFZ(i)) simple infinite eigenvalues; 225C i=1 226C 227C - there are INFZ(i) infinite eigenvalues with multiplicity i+1, 228C for i = 1, ..., DINFZ. 229C 230C The left Kronecker indices are: 231C 232C [ 0 0 ... 0 | 1 1 ... 1 | .... | NKRONL ... NKRONL ] 233C |<- KRONL(1) ->|<- KRONL(2) ->| |<- KRONL(NKRONL) ->| 234C 235C CONTRIBUTOR 236C 237C A. Varga, German Aerospace Center, DLR Oberpfaffenhofen. 238C May 1999. Based on the RASP routine SRISEP. 239C 240C REVISIONS 241C 242C V. Sima, Research Institute for Informatics, Bucharest, Sep. 1999, 243C Jan. 2009, Apr. 2009. 244C A. Varga, DLR Oberpfaffenhofen, March 2002. 245C V. Sima, Jan. 2010, following Bujanovic and Drmac's suggestion. 246C 247C KEYWORDS 248C 249C Generalized eigenvalue problem, Kronecker indices, multivariable 250C system, orthogonal transformation, structural invariant. 251C 252C ****************************************************************** 253C 254C .. Parameters .. 255 INTEGER IMAX, IMIN 256 PARAMETER ( IMAX = 1, IMIN = 2 ) 257 DOUBLE PRECISION ONE, ZERO 258 PARAMETER ( ONE = 1.0D0, ZERO = 0.0D0 ) 259C .. Scalar Arguments .. 260 INTEGER DINFZ, INFO, LDABCD, LDE, LDWORK, M, N, NINFZ, 261 $ NKRONL, NR, P, PR 262 DOUBLE PRECISION SVLMAX, TOL 263 LOGICAL FIRST 264C .. Array Arguments .. 265 INTEGER INFZ( * ), IWORK(*), KRONL( * ) 266 DOUBLE PRECISION ABCD( LDABCD, * ), DWORK( * ), E( LDE, * ) 267C .. Local Scalars .. 268 LOGICAL LQUERY 269 INTEGER I, ICOL, ILAST, IRC, IROW, ISMAX, ISMIN, ITAU, 270 $ J, JLAST, JWORK1, JWORK2, K, MN, MN1, MNR, 271 $ MNTAU, MP1, MPM, MUI, MUIM1, N1, NB, NBLCKS, 272 $ PN, RANK, RO, RO1, SIGMA, TAUI, WRKOPT 273 DOUBLE PRECISION C, C1, C2, RCOND, S, S1, S2, SMAX, SMAXPR, 274 $ SMIN, SMINPR, T, TOLZ, TT 275C .. Local Arrays .. 276 DOUBLE PRECISION DUM(1), SVAL(3) 277C .. External Functions .. 278 INTEGER IDAMAX, ILAENV 279 DOUBLE PRECISION DLAMCH, DNRM2 280 EXTERNAL DLAMCH, DNRM2, IDAMAX, ILAENV 281C .. External Subroutines .. 282 EXTERNAL DCOPY, DLAIC1, DLAPMT, DLARFG, DLARTG, DLASET, 283 $ ODLTZM, DORMQR, DROT, DSWAP, MB03OY, XERBLA 284C .. Intrinsic Functions .. 285 INTRINSIC ABS, DBLE, INT, MAX, MIN, SQRT 286C .. Executable Statements .. 287C 288C Test the input parameters. 289C 290 LQUERY = ( LDWORK.EQ.-1 ) 291 INFO = 0 292 PN = P + N 293 MN = M + N 294 MPM = MIN( P, M ) 295 IF( N.LT.0 ) THEN 296 INFO = -2 297 ELSE IF( M.LT.0 .OR. ( .NOT.FIRST .AND. M.GT.P ) ) THEN 298 INFO = -3 299 ELSE IF( P.LT.0 ) THEN 300 INFO = -4 301 ELSE IF( SVLMAX.LT.ZERO ) THEN 302 INFO = -5 303 ELSE IF( LDABCD.LT.MAX( 1, PN ) ) THEN 304 INFO = -7 305 ELSE IF( LDE.LT.MAX( 1, N ) ) THEN 306 INFO = -9 307 ELSE IF( TOL.GT.ONE ) THEN 308 INFO = -17 309 ELSE 310 WRKOPT = MAX( 1, 5*P ) 311 IF( P.GT.0 ) THEN 312 IF( M.GT.0 ) THEN 313 WRKOPT = MAX( WRKOPT, MN-1 ) 314 IF( FIRST ) THEN 315 WRKOPT = MAX( WRKOPT, MPM + MAX( 3*M-1, N ) ) 316 IF( LQUERY ) THEN 317 NB = MIN( 64, ILAENV( 1, 'DORMQR', 'LT', P, N, 318 $ MPM, -1 ) ) 319 WRKOPT = MAX( WRKOPT, MPM + MAX( 1, N )*NB ) 320 END IF 321 END IF 322 END IF 323 END IF 324 IF( LDWORK.LT.WRKOPT .AND. .NOT.LQUERY ) THEN 325 INFO = -20 326 END IF 327 END IF 328 IF( INFO.NE.0 ) THEN 329 CALL XERBLA( 'AG08BY', -INFO ) 330 RETURN 331 ELSE IF( LQUERY ) THEN 332 DWORK(1) = WRKOPT 333 RETURN 334 END IF 335C 336C Initialize output variables. 337C 338 PR = P 339 NR = N 340 DINFZ = 0 341 NINFZ = 0 342 NKRONL = 0 343C 344C Quick return if possible. 345C 346 IF( P.EQ.0 ) THEN 347 DWORK(1) = ONE 348 RETURN 349 END IF 350 IF( N.EQ.0 .AND. M.EQ.0 ) THEN 351 PR = 0 352 NKRONL = 1 353 KRONL(1) = P 354 DWORK(1) = ONE 355 RETURN 356 END IF 357C 358 TOLZ = SQRT( DLAMCH( 'Epsilon' ) ) 359 RCOND = TOL 360 IF( RCOND.LE.ZERO ) THEN 361C 362C Use the default tolerance in rank determination. 363C 364 RCOND = DBLE( PN*MN )*DLAMCH( 'EPSILON' ) 365 END IF 366C 367C The D matrix is (RO+SIGMA)-by-M, where RO = P - SIGMA and 368C SIGMA = 0 for FIRST = .TRUE. and SIGMA = M for FIRST = .FALSE.. 369C The leading (RO+SIGMA)-by-SIGMA submatrix of D has full column 370C rank, with the trailing SIGMA-by-SIGMA submatrix upper triangular. 371C 372 IF( FIRST ) THEN 373 SIGMA = 0 374 ELSE 375 SIGMA = M 376 END IF 377 RO = P - SIGMA 378 MP1 = M + 1 379 MUI = 0 380 DUM(1) = ZERO 381C 382 ITAU = 1 383 JWORK1 = ITAU + MPM 384 ISMIN = 2*P + 1 385 ISMAX = ISMIN + P 386 JWORK2 = ISMAX + P 387 NBLCKS = 0 388 WRKOPT = 1 389C 390 10 IF( PR.EQ.0 ) GO TO 90 391C 392C (NR+1,ICOL+1) points to the current position of matrix D. 393C 394 RO1 = RO 395 MNR = M + NR 396 IF( M.GT.0 ) THEN 397C 398C Compress rows of D; first exploit the trapezoidal shape of the 399C (RO+SIGMA)-by-SIGMA matrix in the first SIGMA columns of D; 400C compress the first SIGMA columns without column pivoting: 401C 402C ( x x x x x ) ( x x x x x ) 403C ( x x x x x ) ( 0 x x x x ) 404C ( x x x x x ) - > ( 0 0 x x x ) 405C ( 0 x x x x ) ( 0 0 0 x x ) 406C ( 0 0 x x x ) ( 0 0 0 x x ) 407C 408C where SIGMA = 3 and RO = 2. 409C Workspace: need maximum M+N-1. 410C 411 IROW = NR 412 DO 20 ICOL = 1, SIGMA 413 IROW = IROW + 1 414 CALL DLARFG( RO+1, ABCD(IROW,ICOL), ABCD(IROW+1,ICOL), 1, 415 $ T ) 416 CALL ODLTZM( 'L', RO+1, MNR-ICOL, ABCD(IROW+1,ICOL), 1, T, 417 $ ABCD(IROW,ICOL+1), ABCD(IROW+1,ICOL+1), 418 $ LDABCD, DWORK ) 419 CALL DCOPY( PR-ICOL, DUM, 0, ABCD(IROW+1,ICOL), 1 ) 420 20 CONTINUE 421 WRKOPT = MAX( WRKOPT, MN - 1 ) 422C 423 IF( FIRST ) THEN 424C 425C Continue with Householder with column pivoting. 426C 427C ( x x x x x ) ( x x x x x ) 428C ( 0 x x x x ) ( 0 x x x x ) 429C ( 0 0 x x x ) - > ( 0 0 x x x ) 430C ( 0 0 0 x x ) ( 0 0 0 x x ) 431C ( 0 0 0 x x ) ( 0 0 0 0 0 ) 432C 433C Real workspace: need maximum min(P,M)+3*M-1; 434C Integer workspace: need maximum M. 435C 436 IROW = MIN( NR+SIGMA+1, PN ) 437 ICOL = MIN( SIGMA+1, M ) 438 CALL MB03OY( RO1, M-SIGMA, ABCD(IROW,ICOL), LDABCD, 439 $ RCOND, SVLMAX, RANK, SVAL, IWORK, DWORK(ITAU), 440 $ DWORK(JWORK1), INFO ) 441 WRKOPT = MAX( WRKOPT, JWORK1 + 3*M - 2 ) 442C 443C Apply the column permutations to B and part of D. 444C 445 CALL DLAPMT( .TRUE., NR+SIGMA, M-SIGMA, ABCD(1,ICOL), 446 $ LDABCD, IWORK ) 447C 448 IF( RANK.GT.0 ) THEN 449C 450C Apply the Householder transformations to the submatrix C. 451C Workspace: need maximum min(P,M) + N; 452C prefer maximum min(P,M) + N*NB. 453C 454 CALL DORMQR( 'Left', 'Transpose', RO1, NR, RANK, 455 $ ABCD(IROW,ICOL), LDABCD, DWORK(ITAU), 456 $ ABCD(IROW,MP1), LDABCD, DWORK(JWORK1), 457 $ LDWORK-JWORK1+1, INFO ) 458 WRKOPT = MAX( WRKOPT, JWORK1 + INT( DWORK(JWORK1) ) - 1 ) 459 CALL DLASET( 'Lower', RO1-1, MIN( RO1-1, RANK ), ZERO, 460 $ ZERO, ABCD(MIN( IROW+1, PN ),ICOL), LDABCD ) 461 RO1 = RO1 - RANK 462 END IF 463 END IF 464C 465C Terminate if Dr has maximal row rank. 466C 467 IF( RO1.EQ.0 ) GO TO 90 468C 469 END IF 470C 471C Update SIGMA. 472C 473 SIGMA = PR - RO1 474C 475 NBLCKS = NBLCKS + 1 476 TAUI = RO1 477C 478C Compress the columns of current C to separate a TAUI-by-MUI 479C full column rank block. 480C 481 IF( NR.EQ.0 ) THEN 482C 483C Finish for zero state dimension. 484C 485 PR = SIGMA 486 RANK = 0 487 ELSE 488C 489C Perform RQ-decomposition with row pivoting on the current C 490C while keeping E upper triangular. 491C The current C is the TAUI-by-NR matrix delimited by rows 492C IRC+1 to IRC+TAUI and columns M+1 to M+NR of ABCD. 493C The rank of current C is computed in MUI. 494C Workspace: need maximum 5*P. 495C 496 IRC = NR + SIGMA 497 N1 = NR 498 IF( TAUI.GT.1 ) THEN 499C 500C Compute norms. 501C 502 DO 30 I = 1, TAUI 503 DWORK(I) = DNRM2( NR, ABCD(IRC+I,MP1), LDABCD ) 504 DWORK(P+I) = DWORK(I) 505 30 CONTINUE 506 END IF 507C 508 RANK = 0 509 MNTAU = MIN( TAUI, NR ) 510C 511C ICOL and IROW will point to the current pivot position in C. 512C 513 ILAST = NR + PR 514 JLAST = M + NR 515 IROW = ILAST 516 ICOL = JLAST 517 I = TAUI 518 40 IF( RANK.LT.MNTAU ) THEN 519 MN1 = M + N1 520C 521C Pivot if necessary. 522C 523 IF( I.NE.1 ) THEN 524 J = IDAMAX( I, DWORK, 1 ) 525 IF( J.NE.I ) THEN 526 DWORK(J) = DWORK(I) 527 DWORK(P+J) = DWORK(P+I) 528 CALL DSWAP( N1, ABCD(IROW,MP1), LDABCD, 529 $ ABCD(IRC+J,MP1), LDABCD ) 530 END IF 531 END IF 532C 533C Zero elements left to ABCD(IROW,ICOL). 534C 535 DO 50 K = 1, N1-1 536 J = M + K 537C 538C Rotate columns J, J+1 to zero ABCD(IROW,J). 539C 540 T = ABCD(IROW,J+1) 541 CALL DLARTG( T, ABCD(IROW,J), C, S, ABCD(IROW,J+1) ) 542 ABCD(IROW,J) = ZERO 543 CALL DROT( IROW-1, ABCD(1,J+1), 1, ABCD(1,J), 1, C, S ) 544 CALL DROT( K+1, E(1,K+1), 1, E(1,K), 1, C, S ) 545C 546C Rotate rows K, K+1 to zero E(K+1,K). 547C 548 T = E(K,K) 549 CALL DLARTG( T, E(K+1,K), C, S, E(K,K) ) 550 E(K+1,K) = ZERO 551 CALL DROT( N1-K, E(K,K+1), LDE, E(K+1,K+1), LDE, C, S ) 552 CALL DROT( MN1, ABCD(K,1), LDABCD, ABCD(K+1,1), LDABCD, 553 $ C, S ) 554 50 CONTINUE 555C 556 IF( RANK.EQ.0 ) THEN 557C 558C Initialize; exit if matrix is zero (RANK = 0). 559C 560 SMAX = ABS( ABCD(ILAST,JLAST) ) 561 IF ( SMAX.EQ.ZERO ) GO TO 80 562 SMIN = SMAX 563 SMAXPR = SMAX 564 SMINPR = SMIN 565 C1 = ONE 566 C2 = ONE 567 ELSE 568C 569C One step of incremental condition estimation. 570C 571 CALL DCOPY( RANK, ABCD(IROW,ICOL+1), LDABCD, 572 $ DWORK(JWORK2), 1 ) 573 CALL DLAIC1( IMIN, RANK, DWORK( ISMIN ), SMIN, 574 $ DWORK(JWORK2), ABCD(IROW,ICOL), SMINPR, S1, 575 $ C1 ) 576 CALL DLAIC1( IMAX, RANK, DWORK( ISMAX ), SMAX, 577 $ DWORK(JWORK2), ABCD(IROW,ICOL), SMAXPR, S2, 578 $ C2 ) 579 WRKOPT = MAX( WRKOPT, 5*P ) 580 END IF 581C 582C Check the rank; finish the loop if rank loss occurs. 583C 584 IF( SVLMAX*RCOND.LE.SMAXPR ) THEN 585 IF( SVLMAX*RCOND.LE.SMINPR ) THEN 586 IF( SMAXPR*RCOND.LE.SMINPR ) THEN 587C 588C Finish the loop if last row. 589C 590 IF( N1.EQ.0 ) THEN 591 RANK = RANK + 1 592 GO TO 80 593 END IF 594C 595 IF( N1.GT.1 ) THEN 596C 597C Update norms. 598C 599 IF( I-1.GT.1 ) THEN 600 DO 60 J = 1, I - 1 601 IF( DWORK(J).NE.ZERO ) THEN 602 T = ABS( ABCD(IRC+J,ICOL) ) / DWORK(J) 603 T = MAX( ( ONE + T )*( ONE - T ), ZERO) 604 TT = T*( DWORK(J)/DWORK(P+J) )**2 605 IF( TT.GT.TOLZ ) THEN 606 DWORK(J) = DWORK(J)*SQRT( T ) 607 ELSE 608 DWORK(J) = DNRM2( N1-1, 609 $ ABCD(IRC+J,MP1), LDABCD ) 610 DWORK(P+J) = DWORK(J) 611 END IF 612 END IF 613 60 CONTINUE 614 END IF 615 END IF 616C 617 DO 70 J = 1, RANK 618 DWORK( ISMIN+J-1 ) = S1*DWORK( ISMIN+J-1 ) 619 DWORK( ISMAX+J-1 ) = S2*DWORK( ISMAX+J-1 ) 620 70 CONTINUE 621C 622 DWORK( ISMIN+RANK ) = C1 623 DWORK( ISMAX+RANK ) = C2 624 SMIN = SMINPR 625 SMAX = SMAXPR 626 RANK = RANK + 1 627 ICOL = ICOL - 1 628 IROW = IROW - 1 629 N1 = N1 - 1 630 I = I - 1 631 GO TO 40 632 END IF 633 END IF 634 END IF 635 END IF 636 END IF 637C 638 80 CONTINUE 639 MUI = RANK 640 NR = NR - MUI 641 PR = SIGMA + MUI 642C 643C Set number of left Kronecker blocks of order (i-1)-by-i. 644C 645 KRONL(NBLCKS) = TAUI - MUI 646C 647C Set number of infinite divisors of order i-1. 648C 649 IF( FIRST .AND. NBLCKS.GT.1 ) 650 $ INFZ(NBLCKS-1) = MUIM1 - TAUI 651 MUIM1 = MUI 652 RO = MUI 653C 654C Continue reduction if rank of current C is positive. 655C 656 IF( MUI.GT.0 ) 657 $ GO TO 10 658C 659C Determine the maximal degree of infinite zeros and 660C the number of infinite zeros. 661C 662 90 CONTINUE 663 IF( FIRST ) THEN 664 IF( MUI.EQ.0 ) THEN 665 DINFZ = MAX( 0, NBLCKS - 1 ) 666 ELSE 667 DINFZ = NBLCKS 668 INFZ(NBLCKS) = MUI 669 END IF 670 K = DINFZ 671 DO 100 I = K, 1, -1 672 IF( INFZ(I).NE.0 ) GO TO 110 673 DINFZ = DINFZ - 1 674 100 CONTINUE 675 110 CONTINUE 676 DO 120 I = 1, DINFZ 677 NINFZ = NINFZ + INFZ(I)*I 678 120 CONTINUE 679 END IF 680C 681C Determine the maximal order of left elementary Kronecker blocks. 682C 683 NKRONL = NBLCKS 684 DO 130 I = NBLCKS, 1, -1 685 IF( KRONL(I).NE.0 ) GO TO 140 686 NKRONL = NKRONL - 1 687 130 CONTINUE 688 140 CONTINUE 689C 690 DWORK(1) = WRKOPT 691 RETURN 692C *** Last line of AG08BY *** 693 END 694