1 SUBROUTINE FSTAIR (A, E, Q, Z, M, N, ISTAIR, RANKE, TOL, 2 * NBLCKS, IMUK, INUK, IMUK0, INUK0, 3 * MNEI, WRK, IWRK,IERR) 4C PURPOSE: 5C 6C Given a pencil sE-A where matrix E is in column echelon form the 7C subroutine FSTAIR computes according to the wishes of the user a 8C unitary transformed pencil Q(sE-A)Z which is more or less similar 9C to the generalized Schur form of the pencil sE-A. 10C The subroutine yields also part of the Kronecker structure of 11C the given pencil. 12C 13C CONTRIBUTOR: Th.G.J. Beelen (Philips Glass Eindhoven). 14C Copyright SLICOT 15C 16C REVISIONS: 1988, February 02. 17C 18C*********************************************************************** 19C 20C Philips Glass Eindhoven 21C 5600 MD Eindhoven, Netherlands 22C 23C*********************************************************************** 24C FSTAIR - SLICOT Library Routine Document 25C 26C 1 PURPOSE: 27C 28C Given a pencil sE-A where matrix E is in column echelon form the 29C subroutine FSTAIR computes according to the wishes of the user a 30C unitary transformed pencil Q(sE-A)Z which is more or less similar 31C to the generalized Schur form of the pencil sE-A. The computed form 32C yields part of the Kronecker structure of the given pencil. 33C 34C 2 SPECIFICATION: 35C 36C SUBROUTINE FSTAIR(A, LDA, E, Q, LDQ, Z, LDZ, M, N, ISTAIR, RANKE, 37C NBLCKS, IMUK, INUK, IMUK0, INUK0, MNEI, 38C WRK, IWRK, TOL, MODE, IERR) 39C INTEGER LDA, LDQ, LDZ, M, N, RANKE, NBLCKS, MODE, IERR 40C INTEGER ISTAIR(M), IMUK(N), INUK(M+1), IMUK0(N), INUK0(M+1), 41C INTEGER MNEI(4), IWRK(N) 42C DOUBLE PRECISION TOL 43C DOUBLE PRECISION WRK(N) 44C DOUBLE PRECISION A(LDA,N), E(LDA,N), Q(LDQ,M), Z(LDZ,N) 45C 46C 3 ARGUMENT LIST: 47C 48C 3.1 ARGUMENTS IN 49C 50C A - DOUBLE PRECISION array of DIMENSION (LDA,N). 51C The leading M x N part of this array contains the M x N 52C matrix A that has to be row compressed. 53C NOTE: this array is overwritten. 54C 55C LDA - INTEGER 56C LDA is the leading dimension of the arrays A and E. 57C (LDA >= M) 58C 59C E - DOUBLE PRECISION array of DIMENSION (LDA,N). 60C The leading M x N part of this array contains the M x N 61C matrix E which will be transformed equivalent to matrix 62C A. 63C On entry, matrix E has to be in column echelon form. 64C This may be accomplished by subroutine EREDUC. 65C NOTE: this array is overwritten. 66C 67C Q - DOUBLE PRECISION array of DIMENSION (LDQ,M). 68C The leading M x M part of this array contains an M x M 69C unitary row transformation matrix corresponding to the 70C row transformations of the matrices A and E which are 71C needed to transform an arbitrary pencil to a pencil 72C where E is in column echelon form. 73C NOTE: this array is overwritten. 74C 75C LDQ - INTEGER 76C LDQ is the leading dimension of the array Q. 77C (LDQ >= M) 78C 79C Z - DOUBLE PRECISION array of DIMENSION (LDZ,N). 80C The leading N x N part of this array contains an N x N 81C unitary column transformation matrix corresponding to 82C the column transformations of the matrices A and E 83C which are needed to transform an arbitrary pencil to 84C a pencil where E is in column echelon form. 85C NOTE: this array is overwritten. 86C 87C LDZ - INTEGER 88C LDZ is the leading dimension of the array Z. 89C (LDZ >= N) 90C 91C M - INTEGER 92C M is the number of rows of the matrices A, E and Q. 93C 94C N - INTEGER 95C N is the number of columns of the matrices A, E and Z. 96C 97C ISTAIR - INTEGER array of DIMENSION (M). 98C ISTAIR contains the information on the column echelon 99C form of the input matrix E. This may be accomplished 100C by subroutine EREDUC. 101C ISTAIR(i) = + j if the boundary element E(i,j) is a 102C corner point. 103C - j if the boundary element E(i,j) is not 104C a corner point. 105C (i=1,...,M) 106C NOTE: this array is destroyed. 107C 108C RANKE - INTEGER 109C RANKE is the rank of the input matrix E being in column 110C echelon form. 111C 112C 3.2 ARGUMENTS OUT 113C 114C A - DOUBLE PRECISION array of DIMENSION (LDA,N). 115C The leading M x N part of this array contains the M x N 116C matrix A that has been row compressed while keeping E 117C in column echelon form. 118C 119C E - DOUBLE PRECISION array of DIMENSION (LDA,N). 120C The leading M x N part of this array contains the M x N 121C matrix E that has been transformed equivalent to matrix 122C A. 123C 124C Q - DOUBLE PRECISION array of DIMENSION (LDQ,M). 125C The leading M x M part of this array contains the M x M 126C unitary matrix Q which is the product of the input 127C matrix Q and the row transformation matrix which has 128C transformed the rows of the matrices A and E. 129C 130C Z - DOUBLE PRECISION array of DIMENSION (LDZ,N). 131C The leading N x N part of this array contains the N x N 132C unitary matrix Z which is the product of the input 133C matrix Z and the column transformation matrix which has 134C transformed the columns of the matrices A and E. 135C 136C NBLCKS - INTEGER 137C NBLCKS is the number of submatrices having 138C full row rank >= 0 detected in matrix A. 139C 140C IMUK - INTEGER array of DIMENSION (N). 141C Array IMUK contains the column dimensions mu(k) 142C (k=1,...,NBLCKS) of the submatrices having full column 143C rank in the pencil sE(x)-A(x) 144C where x = eps,inf if MODE = 1 or 2 145C eps MODE = 3 . 146C 147C INUK - INTEGER array of DIMENSION (M+1). 148C Array INUK contains the row dimensions nu(k) 149C (k=1,...,NBLCKS) of the submatrices having full row 150C rank in the pencil sE(x)-A(x) 151C where x = eps,inf if MODE = 1 or 2 152C eps MODE = 3 . 153C 154C IMUK0 - INTEGER array of DIMENSION (N). 155C Array IMUK0 contains the column dimensions mu(k) 156C (k=1,...,NBLCKS) of the submatrices having full column 157C rank in the pencil sE(eps,inf)-A(eps,inf). 158C 159C INUK0 - INTEGER array of DIMENSION (M+1). 160C Array INUK0 contains the row dimensions nu(k) 161C (k=1,...,NBLCKS) of the submatrices having full row 162C rank in the pencil sE(eps,inf)-A(eps,inf). 163C 164C MNEI - INTEGER array of DIMENSION (4). 165C If MODE = 3 then 166C MNEI(1) = row dimension of sE(eps)-A(eps) 167C 2 = column dimension of sE(eps)-A(eps) 168C 3 = row dimension of sE(inf)-A(inf) 169C 4 = column dimension of sE(inf)-A(inf) 170C If MODE = 1 or 2 then the array MNEI is empty. 171C 172C 3.3 WORK SPACE 173C 174C WRK - DOUBLE PRECISION array of DIMENSION (N). 175C 176C IWRK - INTEGER array of DIMENSION (N). 177C 178C 3.4 TOLERANCES 179C 180C TOL - DOUBLE PRECISION 181C TOL is the tolerance used when considering matrix 182C elements to be zero. TOL should be set to 183C TOL = RE * max( ||A|| , ||E|| ) + AE , where 184C ||.|| is the Frobenius norm. AE and RE are the absolute 185C and relative accuracy. 186C A recommanded choice is AE = EPS and RE = 100*EPS, 187C where EPS is the machine precision. 188C 189C 3.5 MODE PARAMETERS 190C 191C MODE - INTEGER 192C According to the value of MODE the subroutine FSTAIR 193C computes a generalized Schur form of the pencil sE-A, 194C where the structure of the generalized Schur form is 195C specified more the higher the value of MODE is. 196C (See also 6 DESCRIPTION). 197C 198C 3.6 ERROR INDICATORS 199C 200C IERR - INTEGER 201C On return, IERR contains 0 unless the subroutine 202C has failed. 203C 204C 4 ERROR INDICATORS and WARNINGS: 205C 206C IERR = -1: the value of MODE is not 1, 2 or 3. 207C IERR = 0: succesfull completion. 208C IERR = 1: the algorithm has failed. 209C 210C 5 AUXILARY ROUTINES and COMMON BLOCKS: 211C 212C BAE, SQUAEK, TRIRED from SLICOT. 213C 214C 6 DESCRIPTION: 215C 216C On entry, matrix E is assumed to be in column echelon form. 217C Depending on the value of the parameter MODE, submatrices of A 218C and E will be reduced to a specific form. The higher the value of 219C MODE, the more the submatrices are transformed. 220C 221C Step 1 of the algorithm. 222C If MODE = 1 then subroutine FSTAIR transforms the matrices A and 223C E to the following generalized Schur form by unitary transformations 224C Q1 and Z1, using subroutine BAE. (See also Algorithm 3.2.1 in [1]). 225C 226C | sE(eps,inf)-A(eps,inf) | X | 227C Q1(sE-A)Z1 = |------------------------|------------| 228C | O | sE(r)-A(r) | 229C 230C Here the pencil sE(eps,inf)-A(eps,inf) is in staircase form. 231C This pencil contains all Kronecker column indices and infinite 232C elementary divisors of the pencil sE-A. 233C The pencil sE(r)-A(r) contains all Kronecker row indices and 234C elementary divisors of sE-A. 235C NOTE: X is a pencil. 236C 237C Step 2 of the algorithm. 238C If MODE = 2 then furthermore the submatrices having full row and 239C column rank in the pencil sE(eps,inf)-A(eps,inf) are triangularized 240C by applying unitary transformations Q2 and Z2 to Q1*(sE-A)*Z1. This 241C is done by subroutine TRIRED. (see also Algorithm 3.3.1 [1]). 242C 243C Step 3 of the algorithm. 244C If MODE = 3 then moreover the pencils sE(eps)-A(eps) and 245C sE(inf)-A(inf) are separated in sE(eps,inf)-A(eps,inf) by applying 246C unitary transformations Q3 and Z3 to Q2*Q1*(sE-A)*Z1*Z2. This is 247C done by subroutine SQUAEK. (See also Algorithm 3.3.3 in [1]). 248C We then obtain 249C 250C | sE(eps)-A(eps) | X | X | 251C |----------------|----------------|------------| 252C | O | sE(inf)-A(inf) | X | 253C Q(sE-A)Z = |=================================|============| (1) 254C | | | 255C | O | sE(r)-A(r) | 256C 257C where Q = Q3*Q2*Q1 and Z = Z1*Z2*Z3. 258C The accumulated row and column transformations are multiplied on 259C the left and right respectively with the contents of the arrays Q 260C and Z on entry and the results are stored in Q and Z, respectively. 261C NOTE: the pencil sE(r)-A(r) is not reduced furthermore. 262C 263C Now let sE-A be an arbitrary pencil. This pencil has to be 264C transformed into a pencil with E in column echelon form before 265C calling FSTAIR. This may be accomplished by the subroutine EREDUC. 266C Let the therefore needed unitary row and column transformations be 267C Q0 and Z0, respectively. 268C Let, on entry, the arrays Q and Z contain Q0 and Z0, and let ISTAIR 269C contain the appropiate information. Then, on return with MODE = 3, 270C the contents of the arrays Q and Z are Q3*Q2*Q1*Q0 and Z0*Z1*Z2*Z3 271C which are the transformation matrices that transform the arbitrary 272C pencil sE-A into the form (1). 273C 274C 7 REFERENCES: 275C 276C [1] Th.G.J. Beelen, New Algorithms for Computing the Kronecker 277C structure of a Pencil with Applications to Systems and Control 278C Theory, Ph.D.Thesis, Eindhoven University of Technology, 279C The Netherlands, 1987. 280C 281C 8 NUMERICAL ASPECTS: 282C 283C It is shown in [1] that the algorithm is numerically backward 284C stable. The operations count is proportional to (max(M,N))**3 . 285C 286C 9 FURTHER REMARKS: 287C 288C - The difference mu(k)-nu(k) = # Kronecker blocks of size kx(k+1). 289C The number of these blocks is given by NBLCKS. 290C - The difference nu(k)-mu(k+1) = # infinite elementary divisors of 291C degree k (here mu(NBLCKS+1) := 0). 292C - MNEI(3) = MNEI(4) since pencil sE(inf)-A(inf) is regular. 293C - In the pencil sE(r)-A(r) the pencils sE(f)-A(f) and sE(eta)-A(eta) 294C can be separated by pertransposing the pencil sE(r)-A(r) and 295C using the last part of subroutine FSTAIR. The result has got to be 296C pertransposed again. (For more details see section 3.3.1 in [1]). 297C 298C*********************************************************************** 299C 300C .. Scalar arguments .. 301C 302 INTEGER LDA, LDQ, LDZ, M, N, RANKE, NBLCKS, MODE, IERR 303 DOUBLE PRECISION TOL 304C 305C .. Array arguments .. 306C 307 INTEGER ISTAIR(M), IMUK(N), INUK(M+1), IMUK0(N), INUK0(M+1), 308 * MNEI(4), IWRK(N) 309 DOUBLE PRECISION A(M,N), E(M,N), Q(M,M), Z(N,N), 310 * WRK(N) 311C 312C EXTERNAL SUBROUTINES/FUNCTIONS: 313C 314C BAE, SQUAEK, TRIRED from SLICOT. 315C 316C Local variables. 317C 318 INTEGER MEI, NEI, IFIRA, IFICA, NRA, NCA, JK, RANKA, 319 * ISMUK, ISNUK, I, K 320C 321 LDA=M 322 LDE=M 323 LDQ=M 324 LDZ=N 325 MODE=3 326 IERR = 0 327C 328C A(k) is the submatrix in A that will be row compressed. 329C 330C ISMUK= sum(i=1,..,k) MU(i), ISNUK= sum(i=1,...,k) NU(i), 331C IFIRA, IFICA: first row and first column index of A(k) in A. 332C NRA, NCA: number of rows and columns in A(k). 333C 334 IFIRA = 1 335 IFICA = 1 336 NRA = M 337 NCA = N - RANKE 338 ISNUK = 0 339 ISMUK = 0 340C 341C NBLCKS = # blocks detected in A with full row rank >= 0. 342C 343 NBLCKS = 0 344 K = 0 345C 346C Initialization of the arrays INUK and IMUK. 347C 348 DO 10 I = 1, M + 1 349 INUK(I) = -1 350 10 CONTINUE 351C 352C Note: it is necessary that array INUK has dimension M+1 since it 353C is possible that M = 1 and NBLCKS = 2. 354C Example sE-A = (0 0 s -1). 355C 356 DO 20 I = 1, N 357 IMUK(I) = -1 358 20 CONTINUE 359C 360C Compress the rows of A while keeping E in column echelon form. 361C 362C REPEAT 363C 364 30 K = K + 1 365 CALL BAE(A, LDA, E, Q, LDQ, Z, LDZ, M, N, ISTAIR, IFIRA, 366 * IFICA, NCA, RANKA, WRK, IWRK, TOL) 367 IMUK(K) = NCA 368 ISMUK = ISMUK + NCA 369 370 INUK(K) = RANKA 371 ISNUK = ISNUK + RANKA 372 NBLCKS = NBLCKS + 1 373C 374C If rank of A(k) = nrb then A has full row rank ; 375C JK = first column index (in A) after right most column of 376C matrix A(k+1). 377C (in case A(k+1) is empty, then JK = N+1). 378C 379 IFIRA = 1 + ISNUK 380 IFICA = 1 + ISMUK 381 IF (IFIRA .GT. M) THEN 382 JK = N + 1 383 ELSE 384 JK = IABS(ISTAIR(IFIRA)) 385 END IF 386 NRA = M - ISNUK 387 NCA = JK - 1 - ISMUK 388C 389C If NCA > 0 then there can be done some more row compression 390C of matrix A while keeping matrix E in column echelon form. 391C 392 IF (NCA .GT. 0) GOTO 30 393C UNTIL NCA <= 0 394C 395C Matrix E(k+1) has full column rank since NCA = 0. 396C Reduce A and E by ignoring all rows and columns corresponding 397C to E(k+1). 398C Ignoring these columns in E changes the ranks of the 399C submatrices E(i), (i=1,...,k-1). 400C 401C MEI and NEI are the dimensions of the pencil 402C sE(eps,inf)-A(eps,inf), i.e., the pencil that contains only 403C Kronecker column indices and infinity elementary divisors. 404C 405 MEI = ISNUK 406 NEI = ISMUK 407C 408C Save dimensions of the submatrices having full row or column rank 409C in pencil sE(eps,inf)-A(eps,inf), i.e., copy the arrays 410C IMUK and INUK to IMUK0 and INUK0, respectively. 411C 412 DO 40 I = 1, M + 1 413 INUK0(I) = INUK(I) 414 40 CONTINUE 415C 416 DO 50 I = 1, N 417 IMUK0(I) = IMUK(I) 418 50 CONTINUE 419C 420 IF (MODE .EQ. 1) RETURN 421C 422C Triangularization of the submatrices in A and E. 423C 424 CALL TRIRED(A, LDA, E, Q, LDQ, Z, LDZ, M, N, NBLCKS, INUK, IMUK, 425 * IERR) 426C 427 IF (IERR .NE. 0) then 428c write(6,*) 'error: fstair failed!' 429 return 430 endif 431C 432 IF (MODE .EQ. 2) RETURN 433C 434C Reduction to square submatrices E(k)'s in E. 435C 436 CALL SQUAEK(A, LDA, E, Q, LDQ, Z, LDZ, M, N, NBLCKS, INUK, IMUK, 437 * MNEI) 438C 439 RETURN 440C *** Last line of FSTAIR ********************************************* 441 END 442 SUBROUTINE SQUAEK(A, LDA, E, Q, LDQ, Z, LDZ, M, N, NBLCKS, 443 * INUK, IMUK, MNEI) 444C 445C PURPOSE: 446C 447C On entry, it is assumed that the M by N matrices A and E have 448C been obtained after applying the Algorithms 3.2.1 and 3.3.1 to 449C the pencil s*E - A as described in [1], i.e., 450C 451C | s*E(eps,inf)-A(eps,inf) | X | 452C Q(s*E - A)Z = |-------------------------|-------------| 453C | 0 | s*E(r)-A(r) | 454C 455C Here the pencil s*E(eps,inf)-A(eps,inf) is in staircase form. 456C This pencil contains all Kronecker column indices and infinite 457C elementary divisors of the pencil s*E - A. 458C The pencil s*E(r)-A(r) contains all Kronecker row indices and 459C finite elementary divisors of s*E - A. 460C Furthermore, the submatrices having full row and column rank in 461C the pencil s*E(eps,inf)-A(eps,inf) are assumed to be triangu- 462C larized. 463C Subroutine SQUAEK separates the pencils s*E(eps)-A(eps) and 464C s*E(inf)-A(inf) in s*E(eps,inf)-A(eps,inf) using Algorithm 3.3.3 465C in [1]. The result then is 466C 467C Q(s*E - A)Z = 468C 469C | s*E(eps)-A(eps) | X | X | 470C |-----------------|-----------------|-------------| 471C | 0 | s*E(inf)-A(inf) | X | 472C |===================================|=============| 473C | | | 474C | 0 | s*E(r)-A(r) | 475C 476C Note that the pencil s*E(r)-A(r) is not reduced furthermore. 477C REMARK: This routine is intended to be called only from the 478C SLICOT routine FSTAIR. 479C 480C PARAMETERS: 481C 482C A - DOUBLE PRECISION array of dimension (LDA,N). 483C On entry, it contains the matrix AA to be reduced. 484C On return, it contains the transformed matrix AA. 485C E - DOUBLE PRECISION array of dimension (LDA,N). 486C On entry, it contains the matrix EE to be reduced. 487C On return, it contains the transformed matrix EE. 488C Q - DOUBLE PRECISION array of dimension (LDQ,M). 489C On entry, Q contains the row transformations corresponding to 490C to the input matrices A and E. 491C On return, it contains the product of the input matrix Q and 492C the row transformation matrix that has transformed the rows 493C of the matrices A and E. 494C Z - DOUBLE PRECISION array of dimension (LDZ,N). 495C On entry, Z contains the column transformations corresponding 496C to the input matrices A and E. 497C On return, it contains the product of the input matrix Z and 498C the column transformation matrix that has transformed the 499C columns of the matrices A and E. 500C M - INTEGER. 501C Number of rows of A and E. 1 <= M <= LDA. 502C N - INTEGER. 503C Number of columns of A and E. N >= 1. 504C NBLCKS - INTEGER. 505C Number of submatrices having full row rank >=0 in A(eps,inf). 506C INUK - INTEGER array of dimension (NBLCKS). 507C On entry, INUK contains the row dimensions nu(k), 508C (k=1,..,NBLCKS) of the submatrices having full row rank in the 509C pencil s*E(eps,inf)-A(eps,inf). 510C On return, INUK contains the row dimensions nu(k), 511C (k=1,..,NBLCKS) of the submatrices having full row rank in the 512C pencil s*E(eps)-A(eps). 513C IMUK - INTEGER array of dimension (NBLCKS). 514C On entry, IMUK contains the column dimensions mu(k), 515C (k=1,..,NBLCKS) of the submatrices having full column rank in 516C the pencil s*E(eps,inf)-A(eps,inf). 517C On return, IMUK contains the column dimensions mnu(k), 518C (k=1,..,NBLCKS) of the submatrices having full column rank in 519C the pencil s*E(eps)-A(eps). 520C MNEI - INTEGER array of dimension (4). 521C MNEI(1) = MEPS = row dimension of s*E(eps)-A(eps), 522C 2 = NEPS = column dimension of s*E(eps)-A(eps), 523C 3 = MINF = row dimension of s*E(inf)-A(inf), 524C 4 = NINF = column dimension of s*E(inf)-A(inf). 525C 526C REFERENCES: 527C 528C [1] Th.G.J. Beelen, New Algorithms for Computing the Kronecker 529C structure of a Pencil with Applications to Systems and 530C Control Theory, Ph.D.Thesis, Eindhoven University of 531C Technology, The Netherlands, 1987. 532C 533C CONTRIBUTOR: Th.G.J. Beelen (Philips Glas Eindhoven) 534C 535C REVISIONS: 1988, February 02. 536C 537C Specification of the parameters. 538C 539C .. Scalar arguments .. 540C 541 INTEGER LDA, LDQ, LDZ, M, N, NBLCKS 542C 543C .. Array arguments .. 544C 545 DOUBLE PRECISION A(LDA,N), E(LDA,N), Q(LDQ,M), Z(LDZ,N) 546 INTEGER INUK(NBLCKS), IMUK(NBLCKS), MNEI(4) 547C 548C EXTERNAL SUBROUTINES: 549C 550C DGIV, DROTI from SLICOT. 551C 552C Local variables. 553C 554 DOUBLE PRECISION SC, SS 555 INTEGER SK1P1, TK1P1, TP1, TP 556 INTEGER ISMUK, ISNUK, MUKP1, MUK, NUK 557 INTEGER IP, J, K, MUP, MUP1, NUP, NELM 558 INTEGER MEPS, NEPS, MINF, NINF 559 INTEGER RA, CA, RE, CE, RJE, CJE, CJA 560C 561C Initialisation. 562C 563 ISMUK = 0 564 ISNUK = 0 565 DO 10 K = 1, NBLCKS 566 ISMUK = ISMUK + IMUK(K) 567 ISNUK = ISNUK + INUK(K) 568 10 CONTINUE 569C 570C MEPS, NEPS are the dimensions of the pencil s*E(eps)-A(eps). 571C MEPS = Sum(k=1,...,nblcks) NU(k), 572C NEPS = Sum(k=1,...,nblcks) MU(k). 573C MINF, NINF are the dimensions of the pencil s*E(inf)-A(inf). 574C 575 MEPS = ISNUK 576 NEPS = ISMUK 577 MINF = 0 578 NINF = 0 579C 580C MUKP1 = mu(k+1). N.B. It is assumed that mu(NBLCKS + 1) = 0. 581C 582 MUKP1 = 0 583C 584 DO 60 K = NBLCKS, 1, -1 585 NUK = INUK(K) 586 MUK = IMUK(K) 587C 588C Reduce submatrix E(k,k+1) to square matrix. 589C NOTE that always NU(k) >= MU(k+1) >= 0. 590C 591C WHILE ( NU(k) > MU(k+1) ) DO 592 20 IF (NUK .GT. MUKP1) THEN 593C 594C sk1p1 = sum(i=k+1,...,p-1) NU(i) 595C tk1p1 = sum(i=k+1,...,p-1) MU(i) 596C ismuk = sum(i=1,...,k) MU(i) 597C tp1 = sum(i=1,...,p-1) MU(i) = ismuk + tk1p1. 598C 599 SK1P1 = 0 600 TK1P1 = 0 601 DO 50 IP = K + 1, NBLCKS 602C 603C Annihilate the elements originally present in the last 604C row of E(k,p+1) and A(k,p). 605C Start annihilating the first MU(p) - MU(p+1) elements by 606C applying column Givens rotations plus interchanging 607C elements. 608C Use original bottom diagonal element of A(k,k) as pivot. 609C Start position pivot in A = (ra,ca). 610C 611 TP1 = ISMUK + TK1P1 612 RA = ISNUK + SK1P1 613 CA = TP1 614C 615 MUP = IMUK(IP) 616 MUP1 = INUK(IP) 617 NUP = MUP1 618C 619 DO 30 J = 1, (MUP - NUP) 620C 621C CJA = current column index of pivot in A. 622C 623 CJA = CA + J - 1 624 CALL DGIV(A(RA,CJA), A(RA,CJA+1), SC, SS) 625C 626C Apply transformations to A- and E-matrix. 627C Interchange columns simultaneously. 628C Update column transformation matrix Z. 629C 630 NELM = RA 631 CALL DROTI(NELM, A(1,CJA), 1, A(1,CJA+1), 1, SC, SS) 632 A(RA,CJA) = 0.0D0 633 CALL DROTI(NELM, E(1,CJA), 1, E(1,CJA+1), 1, SC, SS) 634 CALL DROTI(N, Z(1,CJA), 1, Z(1,CJA+1), 1, SC, SS) 635 30 CONTINUE 636C 637C Annihilate the remaining elements originally present in 638C the last row of E(k,p+1) and A(k,p) by alternatingly 639C applying row and column rotations plus interchanging 640C elements. 641C Use diagonal elements of E(p,p+1) and original bottom 642C diagonal element of A(k,k) as pivots, respectively. 643C (re,ce) and (ra,ca) are the starting positions of the 644C pivots in E and A. 645C 646 RE = RA + 1 647 TP = TP1 + MUP 648 CE = 1 + TP 649 CA = TP - MUP1 650C 651 DO 40 J = 1, MUP1 652C 653C (RJE,CJE) = current position pivot in E. 654C 655 RJE = RE + J - 1 656 CJE = CE + J - 1 657 CJA = CA + J - 1 658C 659C Determine the row transformations. 660C Apply these transformations to E- and A-matrix . 661C Interchange the rows simultaneously. 662C Update row transformation matrix Q. 663C 664 CALL DGIV(E(RJE,CJE), E(RJE-1,CJE), SC, SS) 665 NELM = N - CJE + 1 666 CALL DROTI(NELM, E(RJE,CJE), LDA, E(RJE-1,CJE), LDA, 667 * SC, SS) 668 E(RJE,CJE) = 0.0D0 669 NELM = N - CJA + 1 670 CALL DROTI(NELM, A(RJE,CJA), LDA, A(RJE-1,CJA), LDA, 671 * SC, SS) 672 CALL DROTI(M, Q(RJE,1), LDQ, Q(RJE-1,1), LDQ, SC, SS) 673C 674C Determine the column transformations. 675C Apply these transformations to A- and E-matrix. 676C Interchange the columns simultaneously. 677C Update column transformation matrix Z. 678C 679 CALL DGIV(A(RJE,CJA), A(RJE,CJA+1), SC, SS) 680 NELM = RJE 681 CALL DROTI(NELM, A(1,CJA), 1, A(1,CJA+1), 1, SC, SS) 682 A(RJE,CJA) = 0.0D0 683 CALL DROTI(NELM, E(1,CJA), 1, E(1,CJA+1), 1, SC, SS) 684 CALL DROTI(N, Z(1,CJA), 1, Z(1,CJA+1), 1, SC, SS) 685 40 CONTINUE 686C 687 SK1P1 = SK1P1 + NUP 688 TK1P1 = TK1P1 + MUP 689C 690 50 CONTINUE 691C 692C Reduce A=A(eps,inf) and E=E(eps,inf) by ignoring their last 693C row and right most column. The row and column ignored 694C belong to the pencil s*E(inf)-A(inf). 695C Redefine blocks in new A and E. 696C 697 MUK = MUK - 1 698 NUK = NUK - 1 699 IMUK(K) = MUK 700 INUK(K) = NUK 701 ISMUK = ISMUK - 1 702 ISNUK = ISNUK - 1 703 MEPS = MEPS - 1 704 NEPS = NEPS - 1 705 MINF = MINF + 1 706 NINF = NINF + 1 707C 708 GOTO 20 709 END IF 710C END WHILE 20 711C 712C Now submatrix E(k,k+1) is square. 713C 714C Consider next submatrix (k:=k-1). 715C 716 ISNUK = ISNUK - NUK 717 ISMUK = ISMUK - MUK 718 MUKP1 = MUK 719 60 CONTINUE 720C 721C If mu(NBLCKS) = 0, then the last submatrix counted in NBLCKS is 722C a 0 by 0 (empty) matrix. This "matrix" must be removed. 723C 724 IF (IMUK(NBLCKS) .EQ. 0) NBLCKS = NBLCKS - 1 725C 726C Store dimensions of the pencils s*E(eps)-A(eps) and 727C s*E(inf)-A(inf) in array MNEI. 728C 729 MNEI(1) = MEPS 730 MNEI(2) = NEPS 731 MNEI(3) = MINF 732 MNEI(4) = NINF 733C 734 RETURN 735C *** Last line of SQUAEK ********************************************* 736 END 737** END OF SQUAEKTEXT 738 SUBROUTINE TRIAAK(A, LDA, E, Z, LDZ, N, NRA, NCA, IFIRA, IFICA) 739C 740C PURPOSE: 741C 742C Subroutine TRIAAK reduces a submatrix A(k) of A to upper triangu- 743C lar form by column Givens rotations only. 744C Here A(k) = A(IFIRA:ma,IFICA:na) where ma = IFIRA - 1 + NRA, 745C na = IFICA - 1 + NCA. 746C Matrix A(k) is assumed to have full row rank on entry. Hence, no 747C pivoting is done during the reduction process. See Algorithm 2.3.1 748C and Remark 2.3.4 in [1]. 749C The constructed column transformations are also applied to matrix 750C E(k) = E(1:IFIRA-1,IFICA:na). 751C Note that in E columns are transformed with the same column 752C indices as in A, but with row indices different from those in A. 753C REMARK: This routine is intended to be called only from the 754C SLICOT auxiliary routine TRIRED. 755C 756C PARAMETERS: 757C 758C A - DOUBLE PRECISION array of dimension (LDA,N). 759C On entry, it contains the submatrix A(k) of full row rank 760C to be reduced to upper triangular form. 761C On return, it contains the transformed matrix A(k). 762C E - DOUBLE PRECISION array of dimension (LDA,N). 763C On entry, it contains the submatrix E(k). 764C On return, it contains the transformed matrix E(k). 765C Z - DOUBLE PRECISION array of dimension (LDZ,N). 766C On entry, Z contains the column transformations corresponding 767C to the input matrices A and E. 768C On return, it contains the product of the input matrix Z and 769C the column transformation matrix that has transformed the 770C columns of the matrices A and E. 771C N - INTEGER. 772C Number of columns of A and E. (N >= 1). 773C NRA - INTEGER. 774C Number of rows in A(k) to be transformed (1 <= NRA <= LDA). 775C NCA - INTEGER. 776C Number of columns in A(k) to be transformed (1 <= NCA <= N). 777C IFIRA - INTEGER. 778C Number of first row in A(k) to be transformed. 779C IFICA - INTEGER. 780C Number of first column in A(k) to be transformed. 781C 782C REFERENCES: 783C 784C [1] Th.G.J. Beelen, New Algorithms for Computing the Kronecker 785C structure of a Pencil with Applications to Systems and 786C Control Theory, Ph.D.Thesis, Eindhoven University of 787C Technology, The Netherlands, 1987. 788C 789C CONTRIBUTOR: Th.G.J. Beelen (Philips Glas Eindhoven) 790C 791C REVISIONS: 1988, January 29. 792C 793C Specification of the parameters. 794C 795C .. Scalar arguments .. 796C 797 INTEGER LDA, LDZ, N, NRA, NCA, IFIRA, IFICA 798C 799C .. Array arguments .. 800C 801 DOUBLE PRECISION A(LDA,N), E(LDA,N), Z(LDZ,N) 802C 803C EXTERNAL SUBROUTINES: 804C 805C DROT from BLAS 806C DGIV from SLICOT. 807C 808C Local variables. 809C 810 DOUBLE PRECISION SC, SS 811 INTEGER I, II, J, JJ, JJPVT, IFICA1, IFIRA1, MNI, NELM 812C 813 IFIRA1 = IFIRA - 1 814 IFICA1 = IFICA - 1 815C 816 DO 20 I = NRA, 1, -1 817 II = IFIRA1 + I 818 MNI = NCA - NRA + I 819 JJPVT = IFICA1 + MNI 820 NELM = IFIRA1 + I 821 DO 10 J = MNI - 1, 1, -1 822 JJ = IFICA1 + J 823C 824C Determine the Givens transformation on columns jj and jjpvt. 825C Apply the transformation to these columns to annihilate 826C A(ii,jj) (from rows 1 up to ifira1+i). 827C Apply the transformation also to the E-matrix 828C (from rows 1 up to ifira1). 829C Update column transformation matrix Z. 830C 831 CALL DGIV(A(II,JJPVT), A(II,JJ), SC, SS) 832 CALL DROT(NELM, A(1,JJPVT), 1, A(1,JJ), 1, SC, SS) 833 A(II,JJ) = 0.0D0 834 CALL DROT(IFIRA1, E(1,JJPVT), 1, E(1,JJ), 1, SC, SS) 835 CALL DROT(N, Z(1,JJPVT), 1, Z(1,JJ), 1, SC, SS) 836 10 CONTINUE 837 20 CONTINUE 838C 839 RETURN 840C *** Last line of TRIAAK ********************************************* 841 END 842** END OF TRIAAKTEXT 843*UPTODATE TRIAEKTEXT 844 SUBROUTINE TRIAEK(A, LDA, E, Q, LDQ, M, N, NRE, NCE, IFIRE, 845 * IFICE, IFICA) 846C 847C PURPOSE: 848C 849C Subroutine TRIAEK reduces a submatrix E(k) of E to upper triangu- 850C lar form by row Givens rotations only. 851C Here E(k) = E(IFIRE:me,IFICE:ne), where me = IFIRE - 1 + NRE, 852C ne = IFICE - 1 + NCE. 853C Matrix E(k) is assumed to have full column rank on entry. Hence, 854C no pivoting is done during the reduction process. See Algorithm 855C 2.3.1 and Remark 2.3.4 in [1]. 856C The constructed row transformations are also applied to matrix 857C A(k) = A(IFIRE:me,IFICA:N). 858C Note that in A(k) rows are transformed with the same row indices 859C as in E but with column indices different from those in E. 860C REMARK: This routine is intended to be called only from the 861C SLICOT auxiliary routine TRIRED. 862C 863C PARAMETERS: 864C 865C A - DOUBLE PRECISION array of dimension (LDA,N). 866C On entry, it contains the submatrix A(k). 867C On return, it contains the transformed matrix A(k). 868C E - DOUBLE PRECISION array of dimension (LDA,N). 869C On entry, it contains the submatrix E(k) of full column 870C rank to be reduced to upper triangular form. 871C On return, it contains the transformed matrix E(k). 872C Q - DOUBLE PRECISION array of dimension (LDQ,M). 873C On entry, Q contains the row transformations corresponding 874C to the input matrices A and E. 875C On return, it contains the product of the input matrix Q and 876C the row transformation matrix that has transformed the rows 877C of the matrices A and E. 878C M - INTEGER. 879C Number of rows of A and E. (1 <= M <= LDA). 880C N - INTEGER. 881C Number of columns of A and E. (N >= 1). 882C NRE - INTEGER 883C Number of rows in E to be transformed (1 <= NRE <= M). 884C NCE - INTEGER. 885C Number of columns in E to be transformed (1 <= NCE <= N). 886C IFIRE - INTEGER. 887C Index of first row in E to be transformed. 888C IFICE - INTEGER. 889C Index of first column in E to be transformed. 890C IFICA - INTEGER. 891C Index of first column in A to be transformed. 892C 893C REFERENCES: 894C 895C [1] Th.G.J. Beelen, New Algorithms for Computing the Kronecker 896C structure of a Pencil with Applications to Systems and 897C Control Theory, Ph.D.Thesis, Eindhoven University of 898C Technology, The Netherlands, 1987. 899C 900C CONTRIBUTOR: Th.G.J. Beelen (Philips Glas Eindhoven) 901C 902C REVISIONS: 1988, January 29. 903C 904C Specification of the parameters. 905C 906C .. Scalar arguments .. 907C 908 INTEGER LDA, LDQ, M, N, NRE, NCE, IFIRE, IFICE, IFICA 909C 910C .. Array arguments .. 911C 912 DOUBLE PRECISION A(LDA,N), E(LDA,N), Q(LDQ,M) 913C 914C EXTERNAL SUBROUTINES: 915C 916C DROT from BLAS 917C DGIV from SLICOT. 918C 919C Local variables. 920C 921 DOUBLE PRECISION SC, SS 922 INTEGER I, II, IIPVT, J, JJ, IFICE1, IFIRE1, NELM 923C 924 IFIRE1 = IFIRE - 1 925 IFICE1 = IFICE - 1 926C 927 DO 20 J = 1, NCE 928 JJ = IFICE1 + J 929 IIPVT = IFIRE1 + J 930 DO 10 I = J + 1, NRE 931 II = IFIRE1 + I 932C 933C Determine the Givens transformation on rows ii and iipvt. 934C Apply the transformation to these rows (in whole E-matrix) 935C to annihilate E(ii,jj) (from columns jj up to n). 936C Apply the transformations also to the A-matrix 937C (from columns ifica up to n). 938C Update the row transformation matrix Q. 939C 940 CALL DGIV(E(IIPVT,JJ), E(II,JJ), SC, SS) 941 NELM = N - JJ + 1 942 CALL DROT(NELM, E(IIPVT,JJ), LDA, E(II,JJ), LDA, SC, SS) 943 E(II,JJ) = 0.0D0 944 NELM = N - IFICA + 1 945 CALL DROT(NELM, A(IIPVT,IFICA), LDA, A(II,IFICA), LDA, 946 * SC, SS) 947 CALL DROT(M, Q(IIPVT,1), LDQ, Q(II,1), LDQ, SC, SS) 948 10 CONTINUE 949 20 CONTINUE 950C 951 RETURN 952C *** Last line of TRIAEK ********************************************* 953 END 954** END OF TRIAEKTEXT 955*UPTODATE TRIREDTEXT 956 SUBROUTINE TRIRED(A, LDA, E, Q, LDQ, Z, LDZ, M, N, NBLCKS, 957 * INUK, IMUK, IERR) 958C 959C PURPOSE: 960C 961C On entry, it is assumed that the M by N matrices A and E have 962C been transformed to generalized Schur form by unitary trans- 963C formations (see Algorithm 3.2.1 in [1]), i.e., 964C 965C | s*E(eps,inf)-A(eps,inf) | X | 966C s*E - A = |-------------------------|-------------| . 967C | 0 | s*E(r)-A(r) | 968C 969C Here the pencil s*E(eps,inf)-A(eps,inf) is in staircase form. 970C This pencil contains all Kronecker column indices and infinite 971C elementary divisors of the pencil s*E - A. 972C The pencil s*E(r)-A(r) contains all Kronecker row indices and 973C finite elementary divisors of s*E - A. 974C Subroutine TRIRED performs the triangularization of the sub- 975C matrices having full row and column rank in the pencil 976C s*E(eps,inf)-A(eps,inf) using Algorithm 3.3.1 in [1]. 977C REMARK: This routine is intended to be called only from the 978C SLICOT routine FSTAIR. 979C 980C PARAMETERS: 981C 982C A - DOUBLE PRECISION array of dimension (LDA,N). 983C On entry, it contains the matrix A to be reduced. 984C On return, it contains the transformed matrix A. 985C E - DOUBLE PRECISION array of dimension (LDA,N). 986C On entry, it contains the matrix E to be reduced. 987C On return, it contains the transformed matrix E. 988C Q - DOUBLE PRECISION array of dimension (LDQ,M). 989C On entry, Q contains the row transformations corresponding 990C to the input matrices A and E. 991C On return, it contains the product of the input matrix Q and 992C the row transformation matrix that has transformed the rows 993C of the matrices A and E. 994C Z - DOUBLE PRECISION array of dimension (LDZ,N). 995C On entry, Z contains the column transformations corresponding 996C to the input matrices A and E. 997C On return, it contains the product of the input matrix Z and 998C the column transformation matrix that has transformed the 999C columns of the matrices A and E. 1000C M - INTEGER. 1001C Number of rows in A and E, 1 <= M <= LDA. 1002C N - INTEGER. 1003C Number of columns in A and E, N >= 1. 1004C NBLCKS - INTEGER. 1005C Number of submatrices having full row rank >=0 in A(eps,inf). 1006C INUK - INTEGER array of dimension (NBLCKS). 1007C Array containing the row dimensions nu(k) (k=1,..,NBLCKS) 1008C of the submatrices having full row rank in the pencil 1009C s*E(eps,inf)-A(eps,inf). 1010C IMUK - INTEGER array of dimension (NBLCKS). 1011C Array containing the column dimensions mu(k) (k=1,..,NBLCKS) 1012C of the submatrices having full column rank in the pencil. 1013C IERR - INTEGER. 1014C IERR = 0, successful completion, 1015C 1, incorrect dimensions of a full row rank submatrix, 1016C 2, incorrect dimensions of a full column rank submatrix 1017C 1018C REFERENCES: 1019C 1020C [1] Th.G.J. Beelen, New Algorithms for Computing the Kronecker 1021C structure of a Pencil with Applications to Systems and 1022C Control Theory, Ph.D.Thesis, Eindhoven University of 1023C Technology, The Netherlands, 1987. 1024C 1025C CONTRIBUTOR: Th.G.J. Beelen (Philips Glas Eindhoven) 1026C 1027C REVISIONS: 1988, January 29. 1028C 1029C Specification of the parameters. 1030C 1031C .. Scalar arguments .. 1032C 1033 INTEGER LDA, LDQ, LDZ, M, N, NBLCKS, IERR 1034C 1035C .. Array arguments .. 1036C 1037 DOUBLE PRECISION A(LDA,N), E(LDA,N), Q(LDQ,M), Z(LDZ,N) 1038 INTEGER INUK(NBLCKS), IMUK(NBLCKS) 1039C 1040C EXTERNAL SUBROUTINES: 1041C 1042C TRIAAK, TRIAEK from SLICOT. 1043C 1044C Local variables. 1045C 1046 INTEGER ISMUK, ISNUK1, IFIRA, IFICA, IFIRE, IFICE 1047 INTEGER I, K, MUK, MUKP1, NUK 1048C 1049C ISMUK = sum(i=1,...,k) MU(i), 1050C ISNUK1 = sum(i=1,...,k-1) NU(i). 1051C 1052 ISMUK = 0 1053 ISNUK1 = 0 1054 DO 10 I = 1, NBLCKS 1055 ISMUK = ISMUK + IMUK(I) 1056 ISNUK1 = ISNUK1 + INUK(I) 1057 10 CONTINUE 1058C 1059C NOTE: ISNUK1 has not yet the correct value. 1060C 1061 IERR = 0 1062 MUKP1 = 0 1063 DO 20 K = NBLCKS, 1, -1 1064 MUK = IMUK(K) 1065 NUK = INUK(K) 1066 ISNUK1 = ISNUK1 - NUK 1067C 1068C Determine left upper absolute coordinates of E(k) in E-matrix. 1069C 1070 IFIRE = 1 + ISNUK1 1071 IFICE = 1 + ISMUK 1072C 1073C Determine left upper absolute coordinates of A(k) in A-matrix. 1074C 1075 IFIRA = IFIRE 1076 IFICA = IFICE - MUK 1077C 1078C Reduce E(k) to upper triangular form using Givens 1079C transformations on rows only. Apply the same transformations 1080C to the rows of A(k). 1081C 1082 IF (MUKP1 .GT. NUK) THEN 1083 IERR = 1 1084 RETURN 1085 END IF 1086C 1087 CALL TRIAEK(A, LDA, E, Q, LDQ, M, N, NUK, MUKP1, IFIRE, IFICE, 1088 * IFICA) 1089C 1090C Reduce A(k) to upper triangular form using Givens 1091C transformations on columns only. Apply the same transformations 1092C to the columns in the E-matrix. 1093C 1094 IF (NUK .GT. MUK) THEN 1095 IERR = 2 1096 RETURN 1097 END IF 1098C 1099 CALL TRIAAK(A, LDA, E, Z, LDZ, N, NUK, MUK, IFIRA, IFICA) 1100C 1101 ISMUK = ISMUK - MUK 1102 MUKP1 = MUK 1103 20 CONTINUE 1104C 1105 RETURN 1106C *** Last line of TRIRED ********************************************* 1107 END 1108 SUBROUTINE BAE(A, LDA, E, Q, LDQ, Z, LDZ, M, N, ISTAIR, IFIRA, 1109 * IFICA, NCA, RANK, WRK, IWRK, TOL) 1110C 1111C LIBRARY INDEX: 1112C 1113C 1114C 1115C PURPOSE: 1116C 1117C Let A and E be M x N matrices with E in column echelon form. 1118C Let AA and EE be the following submatrices of A and E: 1119C AA := A(IFIRA : M ; IFICA : N) 1120C EE := E(IFIRA : M ; IFICA : N). 1121C Let Aj and Ej be the following submatrices of AA and EE: 1122C Aj := A(IFIRA : M ; IFICA : IFICA + NCA -1) and 1123C Ej := E(IFIRA : M ; IFICA + NCA : N). 1124C 1125C The subroutine BAE transforms (AA,EE) such that Aj is row 1126C compressed while keeping matrix Ej in column echelon form 1127C (which may be different from the form on entry). 1128C In fact BAE performs the j-th step of Algorithm 3.2.1 in [1]. 1129C Furthermore, BAE determines the rank RANK of the submatrix Ej. 1130C This is equal to the number of corner points in submatrix Ej. 1131C REMARK: This routine is intended to be called only from the 1132C SLICOT routine FSTAIR. 1133C 1134C PARAMETERS: 1135C 1136C A - DOUBLE PRECISION array of DIMENSION (LDA,N). 1137C On entry, A(IFIRA : M ; IFICA : IFICA + NCA -1) contains the 1138C matrix AA. 1139C On return, it contains the matrix AA that has been row com- 1140C pressed while keeping EE in column echelon form. 1141C LDA - INTEGER. 1142C LDA is the leading dimension of the arrays A and E. LDA >= M. 1143C E - DOUBLE PRECISION array of DIMENSION (LDA,N). 1144C On entry, E(IFIRA : M ; IFICA + NCA : N) contains the matrix 1145C EE which is in column echelon form. 1146C On return, it contains the transformed matrix EE which is kept 1147C in column echelon form. 1148C Q - DOUBLE PRECISION array of DIMENSION (LDQ,M). 1149C On entry, the array Q contains the row transformations 1150C corresponding to the input matrices A and E. 1151C On return, it contains the M x M unitary matrix Q which is the 1152C product of the input matrix Q and the row transformation 1153C matrix that has transformed the rows of the matrices A and E. 1154C LDQ - INTEGER. 1155C LDQ is the leading dimension of the matrix Q. LDQ >= M. 1156C Z - DOUBLE PRECISION array of DIMENSION (LDZ,N). 1157C On entry, the array Z contains the column transformations 1158C corresponding to the input matrices A and E. 1159C On return, it contains the N x N unitary matrix Z which is the 1160C product of the input matrix Z and the column transformation 1161C matrix that has transformed the columns of A and E. 1162C LDZ - INTEGER. 1163C LDZ is the leading dimension of the matrix Z. LDZ >= N. 1164C M - INTEGER. 1165C M is the number of rows of the matrices A, E and Q. M >= 1. 1166C N - INTEGER. 1167C N is the number of columns of the matrices A, E and Z. N >= 1. 1168C ISTAIR - INTEGER array of DIMENSION (M). 1169C On entry, ISTAIR contains information on the column echelon 1170C form of the input matrix E as follows: 1171C ISTAIR(i) = + j: the boundary element E(i,j) is a corner point 1172C - j: the boundary element E(i,j) is not a corner 1173C point. 1174C (i=1,...,M) 1175C On return, ISTAIR contains the same information for the trans- 1176C formed matrix E. 1177C IFIRA - INTEGER. 1178C IFIRA is the first row index of the submatrix Aj and Ej in 1179C matrix A and E, respectively. 1180C IFICA - INTEGER. 1181C IFICA and IFICA + NCA are the first column index of the 1182C submatrices Aj and Ej in the matrices A and E, respectively. 1183C NCA - INTEGER. 1184C NCA is the number of columns of the submatrix Aj in A. 1185C RANK - INTEGER. 1186C Numerical rank of the submatrix Ej in E (based on TOL). 1187C WRK - DOUBLE PRECISION array of DIMENSION (N). 1188C A real work space array. 1189C IWRK - INTEGER array of DIMENSION (N). 1190C An integer work space array. 1191C TOL - DOUBLE PECISION. 1192C TOL is the tolerance used when considering matrix elements to 1193C be zero. TOL should be set to RE * max(||A||,||E||) + AE, 1194C where ||.|| is the Frobenius norm. AE and RE are the absolute 1195C and relative accuracy respectively. 1196C A recommanded choice is AE = EPS and RE = 100*EPS, where EPS 1197C is the machine precision. 1198C 1199C REFERENCES: 1200C 1201C [1] Th.G.J. Beelen, New Algorithms for Computing the Kronecker 1202C structure of a Pencil with Applications to Systems and 1203C Control Theory, Ph.D.Thesis, Eindhoven University of 1204C Technology, The Netherlands, 1987. 1205C 1206C CONTRIBUTOR: Th.G.J. Beelen (Philips Glass Eindhoven). 1207C 1208C REVISIONS: 1988, January 29. 1209C 1210C Specification of the parameters. 1211C 1212C .. Scalar arguments .. 1213C 1214 INTEGER LDA, LDQ, LDZ, M, N, IFIRA, IFICA, NCA, RANK 1215 DOUBLE PRECISION TOL 1216C 1217C .. Array arguments .. 1218C 1219 INTEGER ISTAIR(M), IWRK(N) 1220 DOUBLE PRECISION A(LDA,N), E(LDA,N), Q(LDQ,M), Z(LDZ,N), WRK(N) 1221C 1222C EXTERNAL SUBROUTINES/FUNCTIONS: 1223C 1224C IDAMAX, DROT, DSWAP from BLAS. 1225C DGIV from SLICOT. 1226C 1227C Local variables. 1228C 1229 INTEGER I, II, IMX, IP, IR, IST1, IST2, ISTPVT, ITYPE, 1230 * IFIRA1, IFICA1, JPVT, JC1, JC2, NROWS, 1231 * K, K1, KK, L, LSAV, LL, MK1, MXRANK, NELM, MJ, NJ 1232 DOUBLE PRECISION BMXNRM, BMX, SC, SS, EIJPVT 1233 LOGICAL LZERO 1234C 1235C Initialisation. 1236C 1237C NJ = number of columns in submatrix Aj, 1238C MJ = number of rows in submatrices Aj and Ej. 1239C 1240 NJ = NCA 1241 MJ = M + 1 - IFIRA 1242 IFIRA1 = IFIRA - 1 1243 IFICA1 = IFICA - 1 1244 DO 10 I = 1, NJ 1245 IWRK(I) = I 1246 10 CONTINUE 1247 K = 1 1248 LZERO = .FALSE. 1249 RANK = MIN0(NJ,MJ) 1250 MXRANK = RANK 1251C 1252C WHILE (K <= MXRANK) and (LZERO = FALSE) DO 1253 20 IF ((K .LE. MXRANK) .AND. (.NOT.LZERO)) THEN 1254C 1255C Determine column in Aj with largest max-norm. 1256C 1257 BMXNRM = 0.0D0 1258 LSAV = K 1259 DO 30 L = K, NJ 1260C 1261C IMX is relative index in column L of Aj where max el. is 1262C found. 1263C NOTE: the first el. in column L is in row K of matrix Aj. 1264C 1265 KK = IFIRA1 + K 1266 LL = IFICA1 + L 1267 IMX = IDAMAX(MJ - K + 1, A(KK,LL), 1) 1268 BMX = DABS(A(IMX + KK - 1, LL)) 1269 IF (BMX .GT. BMXNRM) THEN 1270 BMXNRM = BMX 1271 LSAV = L 1272 END IF 1273 30 CONTINUE 1274C 1275 IF (BMXNRM .LT. TOL) THEN 1276C 1277C Set submatrix of Aj to zero. 1278C 1279 DO 50 L = K, NJ 1280 LL = IFICA1 + L 1281 DO 40 I = K, MJ 1282 II = IFIRA1 + I 1283 A(II,LL) = 0.0D0 1284 40 CONTINUE 1285 50 CONTINUE 1286 LZERO = .TRUE. 1287 RANK = K - 1 1288 ELSE 1289C 1290C Check whether columns have to be interchanged. 1291C 1292 IF (LSAV .NE. K) THEN 1293C 1294C Interchange the columns in A which correspond to the 1295C columns lsav and k in Aj. Store the permutation in IWRK. 1296C 1297 CALL DSWAP(M, A(1,IFICA1 + K), 1, A(1,IFICA1 + LSAV), 1) 1298 IP = IWRK(LSAV) 1299 IWRK(LSAV) = IWRK(K) 1300 IWRK(K) = IP 1301 END IF 1302C 1303 K1 = K + 1 1304 MK1 = NJ - K + 1 + (N - NCA - IFICA1) 1305 KK = IFICA1 + K 1306C 1307 DO 90 IR = K1, MJ 1308C 1309 I = MJ - IR + K1 1310C 1311C II = absolute row number in A corresponding to row i in 1312C Aj. 1313C 1314 II = IFIRA1 + I 1315C 1316C Construct Givens transformation to annihilate Aj(i,k). 1317C Apply the row transformation to whole matrix A. 1318C (NOT only to Aj). 1319C Update row transformation matrix Q. 1320C 1321 CALL DGIV(A(II - 1,KK), A(II,KK), SC, SS) 1322 CALL DROT(MK1, A(II - 1,KK), LDA, A(II,KK), LDA, SC, SS) 1323 A(II,KK) = 0.0D0 1324 CALL DROT(M, Q(II - 1,1), LDQ, Q(II,1), LDQ, SC, SS) 1325C 1326C Determine boundary type of matrix E at rows II-1 and II. 1327C 1328 IST1 = ISTAIR(II - 1) 1329 IST2 = ISTAIR(II) 1330 IF ((IST1 * IST2) .GT. 0) THEN 1331 IF (IST1 .GT. 0) THEN 1332C 1333C boundary form = (* x) 1334C (0 *) 1335C 1336 ITYPE = 1 1337 ELSE 1338C 1339C boundary form = (x x) 1340C (x x) 1341C 1342 ITYPE = 3 1343 END IF 1344 ELSE 1345 IF (IST1 .LT. 0) THEN 1346C 1347C boundary form = (x x) 1348C (* x) 1349C 1350 ITYPE=2 1351 ELSE 1352C 1353C boundary form = (* x) 1354C (0 x) 1355C 1356 ITYPE = 4 1357 END IF 1358 END IF 1359C 1360C Apply row transformation also to matrix E. 1361C 1362C JC1 = absolute number of the column in E in which stair 1363C element of row i-1 of Ej is present. 1364C JC2 = absolute number of the column in E in which stair 1365C element of row i of Ej is present. 1366C 1367C NOTE: JC1 < JC2 if ITYPE = 1. 1368C JC1 = JC2 if ITYPE = 2, 3 or 4. 1369C 1370 JC1 = IABS(IST1) 1371 JC2 = IABS(IST2) 1372 JPVT = MIN0(JC1,JC2) 1373 NELM = N - JPVT + 1 1374C 1375 CALL DROT(NELM, E(II-1,JPVT), LDA, E(II,JPVT), LDA, 1376 * SC, SS) 1377 EIJPVT = E(II,JPVT) 1378C 1379 GOTO (80, 60, 90, 70), ITYPE 1380C 1381 60 IF (DABS(EIJPVT) .LT. TOL) THEN 1382C (x x) (* x) 1383C Boundary form has been changed from (* x) to (0 x) 1384C 1385 ISTPVT = ISTAIR(II) 1386 ISTAIR(II - 1) = ISTPVT 1387 ISTAIR(II) = -(ISTPVT + 1) 1388 E(II, JPVT) = 0.0D0 1389 END IF 1390 GOTO 90 1391C 1392 70 IF (DABS(EIJPVT) .GE. TOL) THEN 1393C 1394C (* x) (x x) 1395C Boundary form has been changed from (0 x) to (* x) 1396C 1397 ISTPVT = ISTAIR(II - 1) 1398 ISTAIR(II - 1) = -ISTPVT 1399 ISTAIR(II) = ISTPVT 1400 END IF 1401 GOTO 90 1402C 1403C Construct column Givens transformation to annihilate 1404C E(ii,jpvt). 1405C Apply column Givens transformation to matrix E. 1406C (NOT only to Ej). 1407C 1408 80 CALL DGIV(E(II,JPVT + 1), E(II,JPVT), SC, SS) 1409 CALL DROT(II, E(1,JPVT + 1), 1, E(1,JPVT), 1, SC, SS) 1410 E(II,JPVT) = 0.0D0 1411C 1412C Apply this transformation also to matrix A. 1413C (NOT only to Aj). 1414C Update column transformation matrix Z. 1415C 1416 CALL DROT(M, A(1,JPVT + 1), 1, A(1,JPVT), 1, SC, SS) 1417 CALL DROT(N, Z(1,JPVT + 1), 1, Z(1,JPVT), 1, SC, SS) 1418C 1419 90 CONTINUE 1420C 1421 K = K + 1 1422 END IF 1423 GOTO 20 1424 END IF 1425C END WHILE 20 1426C 1427C Permute columns of Aj to original order. 1428C 1429 NROWS = IFIRA1 + RANK 1430 DO 120 I = 1, NROWS 1431 DO 100 K = 1, NJ 1432 KK = IFICA1 + K 1433 WRK(IWRK(K)) = A(I,KK) 1434 100 CONTINUE 1435 DO 110 K = 1, NJ 1436 KK = IFICA1 + K 1437 A(I,KK) = WRK(K) 1438 110 CONTINUE 1439 120 CONTINUE 1440C 1441 RETURN 1442C *** Last line of BAE ************************************************ 1443 END 1444** END OF BAETEXT 1445*UPTODATE DGIVTEXT 1446 SUBROUTINE DGIV(DA, DB, DC, DS) 1447C 1448C LIBRARY INDEX: 1449C 1450C 2.1.4 Decompositions and transformations. 1451C 1452C PURPOSE: 1453C 1454C This routine constructs the Givens transformation 1455C 1456C ( DC DS ) 1457C G = ( ), DC**2 + DS**2 = 1.0D0 , 1458C (-DS DC ) 1459C T T 1460C such that the vector (DA,DB) is transformed into (R,0), i.e., 1461C 1462C ( DC DS ) ( DA ) ( R ) 1463C ( ) ( ) = ( ) 1464C (-DS DC ) ( DB ) ( 0 ) . 1465C 1466C This routine is a modification of the BLAS routine DROTG 1467C (Algorithm 539) in order to leave the arguments DA and DB 1468C unchanged. The value or R is not returned. 1469C 1470C CONTRIBUTOR: P. Van Dooren (PRLB). 1471C 1472C REVISIONS: 1987, November 24. 1473C 1474C Specification of parameters. 1475C 1476C .. Scalar Arguments .. 1477C 1478 DOUBLE PRECISION DA, DB, DC, DS 1479C 1480C Local variables. 1481C 1482 DOUBLE PRECISION R, U, V 1483C 1484 IF (DABS(DA) .GT. DABS(DB)) THEN 1485 U = DA + DA 1486 V = DB/U 1487 R = DSQRT(0.25D0 + V**2) * U 1488 DC = DA/R 1489 DS = V * (DC + DC) 1490 ELSE 1491 IF (DB .NE. 0.0D0) THEN 1492 U = DB + DB 1493 V = DA/U 1494 R = DSQRT(0.25D0 + V**2) * U 1495 DS = DB/R 1496 DC = V * (DS + DS) 1497 ELSE 1498 DC = 1.0D0 1499 DS = 0.0D0 1500 END IF 1501 END IF 1502 RETURN 1503C *** Last line of DGIV *********************************************** 1504 END 1505** END OF DGIVTEXT 1506*UPTODATE DROTITEXT 1507 SUBROUTINE DROTI (N, X, INCX, Y, INCY, C, S) 1508C 1509C LIBRARY INDEX: 1510C 1511C 2.1.4 Decompositions and transfromations. 1512C 1513C PURPOSE: 1514C 1515C The subroutine DROTI performs the Givens transformation, defined 1516C by C (cos) and S (sin), and interchanges the vectors involved, 1517C i.e., 1518C 1519C |X(i)| | 0 1 | | C S | |X(i)| 1520C | | := | | x | | x | |, i = 1,...N. 1521C |Y(i)| | 1 0 | |-S C | |Y(i)| 1522C 1523C REMARK. This routine is a modification of DROT from BLAS. 1524C 1525C CONTRIBUTOR: Th.G.J. Beelen (Philips Glass Eindhoven) 1526C 1527C REVISIONS: 1988, February 02. 1528C 1529C Specification of the parameters. 1530C 1531C .. Scalar argumants .. 1532C 1533 INTEGER INCX, INCY, N 1534 DOUBLE PRECISION C, S 1535C 1536C .. Array arguments .. 1537C 1538 DOUBLE PRECISION X(*), Y(*) 1539C 1540C Local variables. 1541C 1542 DOUBLE PRECISION DTEMP 1543 INTEGER I, IX, IY 1544C 1545 IF (N .LE. 0) RETURN 1546 IF ((INCX.NE.1) .OR. (INCY.NE.1)) THEN 1547C 1548C Code for unequal increments or equal increments not equal to 1. 1549C 1550 IX = 1 1551 IY = 1 1552 IF (INCX.LT.0) IX = (-N+1) * INCX + 1 1553 IF (INCY.LT.0) IY = (-N+1) * INCY + 1 1554 DO 10 I = 1, N 1555 DTEMP = C * Y(IY) - S * X(IX) 1556 Y(IY) = C * X(IX) + S * Y(IY) 1557 X(IX) = DTEMP 1558 IX = IX + INCX 1559 IY = IY + INCY 1560 10 CONTINUE 1561 ELSE 1562C 1563C Code for both increments equal to 1. 1564C 1565 DO 20 I = 1, N 1566 DTEMP = C * Y(I) - S * X(I) 1567 Y(I) = C * X(I) + S * Y(I) 1568 X(I) = DTEMP 1569 20 CONTINUE 1570 END IF 1571 RETURN 1572C *** Last line if DROTI ********************************************** 1573 END 1574