SUBROUTINE FSTAIR (A, E, Q, Z, M, N, ISTAIR, RANKE, TOL, * NBLCKS, IMUK, INUK, IMUK0, INUK0, * MNEI, WRK, IWRK,IERR) C PURPOSE: C C Given a pencil sE-A where matrix E is in column echelon form the C subroutine FSTAIR computes according to the wishes of the user a C unitary transformed pencil Q(sE-A)Z which is more or less similar C to the generalized Schur form of the pencil sE-A. C The subroutine yields also part of the Kronecker structure of C the given pencil. C C CONTRIBUTOR: Th.G.J. Beelen (Philips Glass Eindhoven). C Copyright SLICOT C C REVISIONS: 1988, February 02. C C*********************************************************************** C C Philips Glass Eindhoven C 5600 MD Eindhoven, Netherlands C C*********************************************************************** C FSTAIR - SLICOT Library Routine Document C C 1 PURPOSE: C C Given a pencil sE-A where matrix E is in column echelon form the C subroutine FSTAIR computes according to the wishes of the user a C unitary transformed pencil Q(sE-A)Z which is more or less similar C to the generalized Schur form of the pencil sE-A. The computed form C yields part of the Kronecker structure of the given pencil. C C 2 SPECIFICATION: C C SUBROUTINE FSTAIR(A, LDA, E, Q, LDQ, Z, LDZ, M, N, ISTAIR, RANKE, C NBLCKS, IMUK, INUK, IMUK0, INUK0, MNEI, C WRK, IWRK, TOL, MODE, IERR) C INTEGER LDA, LDQ, LDZ, M, N, RANKE, NBLCKS, MODE, IERR C INTEGER ISTAIR(M), IMUK(N), INUK(M+1), IMUK0(N), INUK0(M+1), C INTEGER MNEI(4), IWRK(N) C DOUBLE PRECISION TOL C DOUBLE PRECISION WRK(N) C DOUBLE PRECISION A(LDA,N), E(LDA,N), Q(LDQ,M), Z(LDZ,N) C C 3 ARGUMENT LIST: C C 3.1 ARGUMENTS IN C C A - DOUBLE PRECISION array of DIMENSION (LDA,N). C The leading M x N part of this array contains the M x N C matrix A that has to be row compressed. C NOTE: this array is overwritten. C C LDA - INTEGER C LDA is the leading dimension of the arrays A and E. C (LDA >= M) C C E - DOUBLE PRECISION array of DIMENSION (LDA,N). C The leading M x N part of this array contains the M x N C matrix E which will be transformed equivalent to matrix C A. C On entry, matrix E has to be in column echelon form. C This may be accomplished by subroutine EREDUC. C NOTE: this array is overwritten. C C Q - DOUBLE PRECISION array of DIMENSION (LDQ,M). C The leading M x M part of this array contains an M x M C unitary row transformation matrix corresponding to the C row transformations of the matrices A and E which are C needed to transform an arbitrary pencil to a pencil C where E is in column echelon form. C NOTE: this array is overwritten. C C LDQ - INTEGER C LDQ is the leading dimension of the array Q. C (LDQ >= M) C C Z - DOUBLE PRECISION array of DIMENSION (LDZ,N). C The leading N x N part of this array contains an N x N C unitary column transformation matrix corresponding to C the column transformations of the matrices A and E C which are needed to transform an arbitrary pencil to C a pencil where E is in column echelon form. C NOTE: this array is overwritten. C C LDZ - INTEGER C LDZ is the leading dimension of the array Z. C (LDZ >= N) C C M - INTEGER C M is the number of rows of the matrices A, E and Q. C C N - INTEGER C N is the number of columns of the matrices A, E and Z. C C ISTAIR - INTEGER array of DIMENSION (M). C ISTAIR contains the information on the column echelon C form of the input matrix E. This may be accomplished C by subroutine EREDUC. C ISTAIR(i) = + j if the boundary element E(i,j) is a C corner point. C - j if the boundary element E(i,j) is not C a corner point. C (i=1,...,M) C NOTE: this array is destroyed. C C RANKE - INTEGER C RANKE is the rank of the input matrix E being in column C echelon form. C C 3.2 ARGUMENTS OUT C C A - DOUBLE PRECISION array of DIMENSION (LDA,N). C The leading M x N part of this array contains the M x N C matrix A that has been row compressed while keeping E C in column echelon form. C C E - DOUBLE PRECISION array of DIMENSION (LDA,N). C The leading M x N part of this array contains the M x N C matrix E that has been transformed equivalent to matrix C A. C C Q - DOUBLE PRECISION array of DIMENSION (LDQ,M). C The leading M x M part of this array contains the M x M C unitary matrix Q which is the product of the input C matrix Q and the row transformation matrix which has C transformed the rows of the matrices A and E. C C Z - DOUBLE PRECISION array of DIMENSION (LDZ,N). C The leading N x N part of this array contains the N x N C unitary matrix Z which is the product of the input C matrix Z and the column transformation matrix which has C transformed the columns of the matrices A and E. C C NBLCKS - INTEGER C NBLCKS is the number of submatrices having C full row rank >= 0 detected in matrix A. C C IMUK - INTEGER array of DIMENSION (N). C Array IMUK contains the column dimensions mu(k) C (k=1,...,NBLCKS) of the submatrices having full column C rank in the pencil sE(x)-A(x) C where x = eps,inf if MODE = 1 or 2 C eps MODE = 3 . C C INUK - INTEGER array of DIMENSION (M+1). C Array INUK contains the row dimensions nu(k) C (k=1,...,NBLCKS) of the submatrices having full row C rank in the pencil sE(x)-A(x) C where x = eps,inf if MODE = 1 or 2 C eps MODE = 3 . C C IMUK0 - INTEGER array of DIMENSION (N). C Array IMUK0 contains the column dimensions mu(k) C (k=1,...,NBLCKS) of the submatrices having full column C rank in the pencil sE(eps,inf)-A(eps,inf). C C INUK0 - INTEGER array of DIMENSION (M+1). C Array INUK0 contains the row dimensions nu(k) C (k=1,...,NBLCKS) of the submatrices having full row C rank in the pencil sE(eps,inf)-A(eps,inf). C C MNEI - INTEGER array of DIMENSION (4). C If MODE = 3 then C MNEI(1) = row dimension of sE(eps)-A(eps) C 2 = column dimension of sE(eps)-A(eps) C 3 = row dimension of sE(inf)-A(inf) C 4 = column dimension of sE(inf)-A(inf) C If MODE = 1 or 2 then the array MNEI is empty. C C 3.3 WORK SPACE C C WRK - DOUBLE PRECISION array of DIMENSION (N). C C IWRK - INTEGER array of DIMENSION (N). C C 3.4 TOLERANCES C C TOL - DOUBLE PRECISION C TOL is the tolerance used when considering matrix C elements to be zero. TOL should be set to C TOL = RE * max( ||A|| , ||E|| ) + AE , where C ||.|| is the Frobenius norm. AE and RE are the absolute C and relative accuracy. C A recommanded choice is AE = EPS and RE = 100*EPS, C where EPS is the machine precision. C C 3.5 MODE PARAMETERS C C MODE - INTEGER C According to the value of MODE the subroutine FSTAIR C computes a generalized Schur form of the pencil sE-A, C where the structure of the generalized Schur form is C specified more the higher the value of MODE is. C (See also 6 DESCRIPTION). C C 3.6 ERROR INDICATORS C C IERR - INTEGER C On return, IERR contains 0 unless the subroutine C has failed. C C 4 ERROR INDICATORS and WARNINGS: C C IERR = -1: the value of MODE is not 1, 2 or 3. C IERR = 0: succesfull completion. C IERR = 1: the algorithm has failed. C C 5 AUXILARY ROUTINES and COMMON BLOCKS: C C BAE, SQUAEK, TRIRED from SLICOT. C C 6 DESCRIPTION: C C On entry, matrix E is assumed to be in column echelon form. C Depending on the value of the parameter MODE, submatrices of A C and E will be reduced to a specific form. The higher the value of C MODE, the more the submatrices are transformed. C C Step 1 of the algorithm. C If MODE = 1 then subroutine FSTAIR transforms the matrices A and C E to the following generalized Schur form by unitary transformations C Q1 and Z1, using subroutine BAE. (See also Algorithm 3.2.1 in [1]). C C | sE(eps,inf)-A(eps,inf) | X | C Q1(sE-A)Z1 = |------------------------|------------| C | O | sE(r)-A(r) | C C Here the pencil sE(eps,inf)-A(eps,inf) is in staircase form. C This pencil contains all Kronecker column indices and infinite C elementary divisors of the pencil sE-A. C The pencil sE(r)-A(r) contains all Kronecker row indices and C elementary divisors of sE-A. C NOTE: X is a pencil. C C Step 2 of the algorithm. C If MODE = 2 then furthermore the submatrices having full row and C column rank in the pencil sE(eps,inf)-A(eps,inf) are triangularized C by applying unitary transformations Q2 and Z2 to Q1*(sE-A)*Z1. This C is done by subroutine TRIRED. (see also Algorithm 3.3.1 [1]). C C Step 3 of the algorithm. C If MODE = 3 then moreover the pencils sE(eps)-A(eps) and C sE(inf)-A(inf) are separated in sE(eps,inf)-A(eps,inf) by applying C unitary transformations Q3 and Z3 to Q2*Q1*(sE-A)*Z1*Z2. This is C done by subroutine SQUAEK. (See also Algorithm 3.3.3 in [1]). C We then obtain C C | sE(eps)-A(eps) | X | X | C |----------------|----------------|------------| C | O | sE(inf)-A(inf) | X | C Q(sE-A)Z = |=================================|============| (1) C | | | C | O | sE(r)-A(r) | C C where Q = Q3*Q2*Q1 and Z = Z1*Z2*Z3. C The accumulated row and column transformations are multiplied on C the left and right respectively with the contents of the arrays Q C and Z on entry and the results are stored in Q and Z, respectively. C NOTE: the pencil sE(r)-A(r) is not reduced furthermore. C C Now let sE-A be an arbitrary pencil. This pencil has to be C transformed into a pencil with E in column echelon form before C calling FSTAIR. This may be accomplished by the subroutine EREDUC. C Let the therefore needed unitary row and column transformations be C Q0 and Z0, respectively. C Let, on entry, the arrays Q and Z contain Q0 and Z0, and let ISTAIR C contain the appropiate information. Then, on return with MODE = 3, C the contents of the arrays Q and Z are Q3*Q2*Q1*Q0 and Z0*Z1*Z2*Z3 C which are the transformation matrices that transform the arbitrary C pencil sE-A into the form (1). C C 7 REFERENCES: C C [1] Th.G.J. Beelen, New Algorithms for Computing the Kronecker C structure of a Pencil with Applications to Systems and Control C Theory, Ph.D.Thesis, Eindhoven University of Technology, C The Netherlands, 1987. C C 8 NUMERICAL ASPECTS: C C It is shown in [1] that the algorithm is numerically backward C stable. The operations count is proportional to (max(M,N))**3 . C C 9 FURTHER REMARKS: C C - The difference mu(k)-nu(k) = # Kronecker blocks of size kx(k+1). C The number of these blocks is given by NBLCKS. C - The difference nu(k)-mu(k+1) = # infinite elementary divisors of C degree k (here mu(NBLCKS+1) := 0). C - MNEI(3) = MNEI(4) since pencil sE(inf)-A(inf) is regular. C - In the pencil sE(r)-A(r) the pencils sE(f)-A(f) and sE(eta)-A(eta) C can be separated by pertransposing the pencil sE(r)-A(r) and C using the last part of subroutine FSTAIR. The result has got to be C pertransposed again. (For more details see section 3.3.1 in [1]). C C*********************************************************************** C C .. Scalar arguments .. C INTEGER LDA, LDQ, LDZ, M, N, RANKE, NBLCKS, MODE, IERR DOUBLE PRECISION TOL C C .. Array arguments .. C INTEGER ISTAIR(M), IMUK(N), INUK(M+1), IMUK0(N), INUK0(M+1), * MNEI(4), IWRK(N) DOUBLE PRECISION A(M,N), E(M,N), Q(M,M), Z(N,N), * WRK(N) C C EXTERNAL SUBROUTINES/FUNCTIONS: C C BAE, SQUAEK, TRIRED from SLICOT. C C Local variables. C INTEGER MEI, NEI, IFIRA, IFICA, NRA, NCA, JK, RANKA, * ISMUK, ISNUK, I, K C LDA=M LDE=M LDQ=M LDZ=N MODE=3 IERR = 0 C C A(k) is the submatrix in A that will be row compressed. C C ISMUK= sum(i=1,..,k) MU(i), ISNUK= sum(i=1,...,k) NU(i), C IFIRA, IFICA: first row and first column index of A(k) in A. C NRA, NCA: number of rows and columns in A(k). C IFIRA = 1 IFICA = 1 NRA = M NCA = N - RANKE ISNUK = 0 ISMUK = 0 C C NBLCKS = # blocks detected in A with full row rank >= 0. C NBLCKS = 0 K = 0 C C Initialization of the arrays INUK and IMUK. C DO 10 I = 1, M + 1 INUK(I) = -1 10 CONTINUE C C Note: it is necessary that array INUK has dimension M+1 since it C is possible that M = 1 and NBLCKS = 2. C Example sE-A = (0 0 s -1). C DO 20 I = 1, N IMUK(I) = -1 20 CONTINUE C C Compress the rows of A while keeping E in column echelon form. C C REPEAT C 30 K = K + 1 CALL BAE(A, LDA, E, Q, LDQ, Z, LDZ, M, N, ISTAIR, IFIRA, * IFICA, NCA, RANKA, WRK, IWRK, TOL) IMUK(K) = NCA ISMUK = ISMUK + NCA INUK(K) = RANKA ISNUK = ISNUK + RANKA NBLCKS = NBLCKS + 1 C C If rank of A(k) = nrb then A has full row rank ; C JK = first column index (in A) after right most column of C matrix A(k+1). C (in case A(k+1) is empty, then JK = N+1). C IFIRA = 1 + ISNUK IFICA = 1 + ISMUK IF (IFIRA .GT. M) THEN JK = N + 1 ELSE JK = IABS(ISTAIR(IFIRA)) END IF NRA = M - ISNUK NCA = JK - 1 - ISMUK C C If NCA > 0 then there can be done some more row compression C of matrix A while keeping matrix E in column echelon form. C IF (NCA .GT. 0) GOTO 30 C UNTIL NCA <= 0 C C Matrix E(k+1) has full column rank since NCA = 0. C Reduce A and E by ignoring all rows and columns corresponding C to E(k+1). C Ignoring these columns in E changes the ranks of the C submatrices E(i), (i=1,...,k-1). C C MEI and NEI are the dimensions of the pencil C sE(eps,inf)-A(eps,inf), i.e., the pencil that contains only C Kronecker column indices and infinity elementary divisors. C MEI = ISNUK NEI = ISMUK C C Save dimensions of the submatrices having full row or column rank C in pencil sE(eps,inf)-A(eps,inf), i.e., copy the arrays C IMUK and INUK to IMUK0 and INUK0, respectively. C DO 40 I = 1, M + 1 INUK0(I) = INUK(I) 40 CONTINUE C DO 50 I = 1, N IMUK0(I) = IMUK(I) 50 CONTINUE C IF (MODE .EQ. 1) RETURN C C Triangularization of the submatrices in A and E. C CALL TRIRED(A, LDA, E, Q, LDQ, Z, LDZ, M, N, NBLCKS, INUK, IMUK, * IERR) C IF (IERR .NE. 0) then c write(6,*) 'error: fstair failed!' return endif C IF (MODE .EQ. 2) RETURN C C Reduction to square submatrices E(k)'s in E. C CALL SQUAEK(A, LDA, E, Q, LDQ, Z, LDZ, M, N, NBLCKS, INUK, IMUK, * MNEI) C RETURN C *** Last line of FSTAIR ********************************************* END SUBROUTINE SQUAEK(A, LDA, E, Q, LDQ, Z, LDZ, M, N, NBLCKS, * INUK, IMUK, MNEI) C C PURPOSE: C C On entry, it is assumed that the M by N matrices A and E have C been obtained after applying the Algorithms 3.2.1 and 3.3.1 to C the pencil s*E - A as described in [1], i.e., C C | s*E(eps,inf)-A(eps,inf) | X | C Q(s*E - A)Z = |-------------------------|-------------| C | 0 | s*E(r)-A(r) | C C Here the pencil s*E(eps,inf)-A(eps,inf) is in staircase form. C This pencil contains all Kronecker column indices and infinite C elementary divisors of the pencil s*E - A. C The pencil s*E(r)-A(r) contains all Kronecker row indices and C finite elementary divisors of s*E - A. C Furthermore, the submatrices having full row and column rank in C the pencil s*E(eps,inf)-A(eps,inf) are assumed to be triangu- C larized. C Subroutine SQUAEK separates the pencils s*E(eps)-A(eps) and C s*E(inf)-A(inf) in s*E(eps,inf)-A(eps,inf) using Algorithm 3.3.3 C in [1]. The result then is C C Q(s*E - A)Z = C C | s*E(eps)-A(eps) | X | X | C |-----------------|-----------------|-------------| C | 0 | s*E(inf)-A(inf) | X | C |===================================|=============| C | | | C | 0 | s*E(r)-A(r) | C C Note that the pencil s*E(r)-A(r) is not reduced furthermore. C REMARK: This routine is intended to be called only from the C SLICOT routine FSTAIR. C C PARAMETERS: C C A - DOUBLE PRECISION array of dimension (LDA,N). C On entry, it contains the matrix AA to be reduced. C On return, it contains the transformed matrix AA. C E - DOUBLE PRECISION array of dimension (LDA,N). C On entry, it contains the matrix EE to be reduced. C On return, it contains the transformed matrix EE. C Q - DOUBLE PRECISION array of dimension (LDQ,M). C On entry, Q contains the row transformations corresponding to C to the input matrices A and E. C On return, it contains the product of the input matrix Q and C the row transformation matrix that has transformed the rows C of the matrices A and E. C Z - DOUBLE PRECISION array of dimension (LDZ,N). C On entry, Z contains the column transformations corresponding C to the input matrices A and E. C On return, it contains the product of the input matrix Z and C the column transformation matrix that has transformed the C columns of the matrices A and E. C M - INTEGER. C Number of rows of A and E. 1 <= M <= LDA. C N - INTEGER. C Number of columns of A and E. N >= 1. C NBLCKS - INTEGER. C Number of submatrices having full row rank >=0 in A(eps,inf). C INUK - INTEGER array of dimension (NBLCKS). C On entry, INUK contains the row dimensions nu(k), C (k=1,..,NBLCKS) of the submatrices having full row rank in the C pencil s*E(eps,inf)-A(eps,inf). C On return, INUK contains the row dimensions nu(k), C (k=1,..,NBLCKS) of the submatrices having full row rank in the C pencil s*E(eps)-A(eps). C IMUK - INTEGER array of dimension (NBLCKS). C On entry, IMUK contains the column dimensions mu(k), C (k=1,..,NBLCKS) of the submatrices having full column rank in C the pencil s*E(eps,inf)-A(eps,inf). C On return, IMUK contains the column dimensions mnu(k), C (k=1,..,NBLCKS) of the submatrices having full column rank in C the pencil s*E(eps)-A(eps). C MNEI - INTEGER array of dimension (4). C MNEI(1) = MEPS = row dimension of s*E(eps)-A(eps), C 2 = NEPS = column dimension of s*E(eps)-A(eps), C 3 = MINF = row dimension of s*E(inf)-A(inf), C 4 = NINF = column dimension of s*E(inf)-A(inf). C C REFERENCES: C C [1] Th.G.J. Beelen, New Algorithms for Computing the Kronecker C structure of a Pencil with Applications to Systems and C Control Theory, Ph.D.Thesis, Eindhoven University of C Technology, The Netherlands, 1987. C C CONTRIBUTOR: Th.G.J. Beelen (Philips Glas Eindhoven) C C REVISIONS: 1988, February 02. C C Specification of the parameters. C C .. Scalar arguments .. C INTEGER LDA, LDQ, LDZ, M, N, NBLCKS C C .. Array arguments .. C DOUBLE PRECISION A(LDA,N), E(LDA,N), Q(LDQ,M), Z(LDZ,N) INTEGER INUK(NBLCKS), IMUK(NBLCKS), MNEI(4) C C EXTERNAL SUBROUTINES: C C DGIV, DROTI from SLICOT. C C Local variables. C DOUBLE PRECISION SC, SS INTEGER SK1P1, TK1P1, TP1, TP INTEGER ISMUK, ISNUK, MUKP1, MUK, NUK INTEGER IP, J, K, MUP, MUP1, NUP, NELM INTEGER MEPS, NEPS, MINF, NINF INTEGER RA, CA, RE, CE, RJE, CJE, CJA C C Initialisation. C ISMUK = 0 ISNUK = 0 DO 10 K = 1, NBLCKS ISMUK = ISMUK + IMUK(K) ISNUK = ISNUK + INUK(K) 10 CONTINUE C C MEPS, NEPS are the dimensions of the pencil s*E(eps)-A(eps). C MEPS = Sum(k=1,...,nblcks) NU(k), C NEPS = Sum(k=1,...,nblcks) MU(k). C MINF, NINF are the dimensions of the pencil s*E(inf)-A(inf). C MEPS = ISNUK NEPS = ISMUK MINF = 0 NINF = 0 C C MUKP1 = mu(k+1). N.B. It is assumed that mu(NBLCKS + 1) = 0. C MUKP1 = 0 C DO 60 K = NBLCKS, 1, -1 NUK = INUK(K) MUK = IMUK(K) C C Reduce submatrix E(k,k+1) to square matrix. C NOTE that always NU(k) >= MU(k+1) >= 0. C C WHILE ( NU(k) > MU(k+1) ) DO 20 IF (NUK .GT. MUKP1) THEN C C sk1p1 = sum(i=k+1,...,p-1) NU(i) C tk1p1 = sum(i=k+1,...,p-1) MU(i) C ismuk = sum(i=1,...,k) MU(i) C tp1 = sum(i=1,...,p-1) MU(i) = ismuk + tk1p1. C SK1P1 = 0 TK1P1 = 0 DO 50 IP = K + 1, NBLCKS C C Annihilate the elements originally present in the last C row of E(k,p+1) and A(k,p). C Start annihilating the first MU(p) - MU(p+1) elements by C applying column Givens rotations plus interchanging C elements. C Use original bottom diagonal element of A(k,k) as pivot. C Start position pivot in A = (ra,ca). C TP1 = ISMUK + TK1P1 RA = ISNUK + SK1P1 CA = TP1 C MUP = IMUK(IP) MUP1 = INUK(IP) NUP = MUP1 C DO 30 J = 1, (MUP - NUP) C C CJA = current column index of pivot in A. C CJA = CA + J - 1 CALL DGIV(A(RA,CJA), A(RA,CJA+1), SC, SS) C C Apply transformations to A- and E-matrix. C Interchange columns simultaneously. C Update column transformation matrix Z. C NELM = RA CALL DROTI(NELM, A(1,CJA), 1, A(1,CJA+1), 1, SC, SS) A(RA,CJA) = 0.0D0 CALL DROTI(NELM, E(1,CJA), 1, E(1,CJA+1), 1, SC, SS) CALL DROTI(N, Z(1,CJA), 1, Z(1,CJA+1), 1, SC, SS) 30 CONTINUE C C Annihilate the remaining elements originally present in C the last row of E(k,p+1) and A(k,p) by alternatingly C applying row and column rotations plus interchanging C elements. C Use diagonal elements of E(p,p+1) and original bottom C diagonal element of A(k,k) as pivots, respectively. C (re,ce) and (ra,ca) are the starting positions of the C pivots in E and A. C RE = RA + 1 TP = TP1 + MUP CE = 1 + TP CA = TP - MUP1 C DO 40 J = 1, MUP1 C C (RJE,CJE) = current position pivot in E. C RJE = RE + J - 1 CJE = CE + J - 1 CJA = CA + J - 1 C C Determine the row transformations. C Apply these transformations to E- and A-matrix . C Interchange the rows simultaneously. C Update row transformation matrix Q. C CALL DGIV(E(RJE,CJE), E(RJE-1,CJE), SC, SS) NELM = N - CJE + 1 CALL DROTI(NELM, E(RJE,CJE), LDA, E(RJE-1,CJE), LDA, * SC, SS) E(RJE,CJE) = 0.0D0 NELM = N - CJA + 1 CALL DROTI(NELM, A(RJE,CJA), LDA, A(RJE-1,CJA), LDA, * SC, SS) CALL DROTI(M, Q(RJE,1), LDQ, Q(RJE-1,1), LDQ, SC, SS) C C Determine the column transformations. C Apply these transformations to A- and E-matrix. C Interchange the columns simultaneously. C Update column transformation matrix Z. C CALL DGIV(A(RJE,CJA), A(RJE,CJA+1), SC, SS) NELM = RJE CALL DROTI(NELM, A(1,CJA), 1, A(1,CJA+1), 1, SC, SS) A(RJE,CJA) = 0.0D0 CALL DROTI(NELM, E(1,CJA), 1, E(1,CJA+1), 1, SC, SS) CALL DROTI(N, Z(1,CJA), 1, Z(1,CJA+1), 1, SC, SS) 40 CONTINUE C SK1P1 = SK1P1 + NUP TK1P1 = TK1P1 + MUP C 50 CONTINUE C C Reduce A=A(eps,inf) and E=E(eps,inf) by ignoring their last C row and right most column. The row and column ignored C belong to the pencil s*E(inf)-A(inf). C Redefine blocks in new A and E. C MUK = MUK - 1 NUK = NUK - 1 IMUK(K) = MUK INUK(K) = NUK ISMUK = ISMUK - 1 ISNUK = ISNUK - 1 MEPS = MEPS - 1 NEPS = NEPS - 1 MINF = MINF + 1 NINF = NINF + 1 C GOTO 20 END IF C END WHILE 20 C C Now submatrix E(k,k+1) is square. C C Consider next submatrix (k:=k-1). C ISNUK = ISNUK - NUK ISMUK = ISMUK - MUK MUKP1 = MUK 60 CONTINUE C C If mu(NBLCKS) = 0, then the last submatrix counted in NBLCKS is C a 0 by 0 (empty) matrix. This "matrix" must be removed. C IF (IMUK(NBLCKS) .EQ. 0) NBLCKS = NBLCKS - 1 C C Store dimensions of the pencils s*E(eps)-A(eps) and C s*E(inf)-A(inf) in array MNEI. C MNEI(1) = MEPS MNEI(2) = NEPS MNEI(3) = MINF MNEI(4) = NINF C RETURN C *** Last line of SQUAEK ********************************************* END ** END OF SQUAEKTEXT SUBROUTINE TRIAAK(A, LDA, E, Z, LDZ, N, NRA, NCA, IFIRA, IFICA) C C PURPOSE: C C Subroutine TRIAAK reduces a submatrix A(k) of A to upper triangu- C lar form by column Givens rotations only. C Here A(k) = A(IFIRA:ma,IFICA:na) where ma = IFIRA - 1 + NRA, C na = IFICA - 1 + NCA. C Matrix A(k) is assumed to have full row rank on entry. Hence, no C pivoting is done during the reduction process. See Algorithm 2.3.1 C and Remark 2.3.4 in [1]. C The constructed column transformations are also applied to matrix C E(k) = E(1:IFIRA-1,IFICA:na). C Note that in E columns are transformed with the same column C indices as in A, but with row indices different from those in A. C REMARK: This routine is intended to be called only from the C SLICOT auxiliary routine TRIRED. C C PARAMETERS: C C A - DOUBLE PRECISION array of dimension (LDA,N). C On entry, it contains the submatrix A(k) of full row rank C to be reduced to upper triangular form. C On return, it contains the transformed matrix A(k). C E - DOUBLE PRECISION array of dimension (LDA,N). C On entry, it contains the submatrix E(k). C On return, it contains the transformed matrix E(k). C Z - DOUBLE PRECISION array of dimension (LDZ,N). C On entry, Z contains the column transformations corresponding C to the input matrices A and E. C On return, it contains the product of the input matrix Z and C the column transformation matrix that has transformed the C columns of the matrices A and E. C N - INTEGER. C Number of columns of A and E. (N >= 1). C NRA - INTEGER. C Number of rows in A(k) to be transformed (1 <= NRA <= LDA). C NCA - INTEGER. C Number of columns in A(k) to be transformed (1 <= NCA <= N). C IFIRA - INTEGER. C Number of first row in A(k) to be transformed. C IFICA - INTEGER. C Number of first column in A(k) to be transformed. C C REFERENCES: C C [1] Th.G.J. Beelen, New Algorithms for Computing the Kronecker C structure of a Pencil with Applications to Systems and C Control Theory, Ph.D.Thesis, Eindhoven University of C Technology, The Netherlands, 1987. C C CONTRIBUTOR: Th.G.J. Beelen (Philips Glas Eindhoven) C C REVISIONS: 1988, January 29. C C Specification of the parameters. C C .. Scalar arguments .. C INTEGER LDA, LDZ, N, NRA, NCA, IFIRA, IFICA C C .. Array arguments .. C DOUBLE PRECISION A(LDA,N), E(LDA,N), Z(LDZ,N) C C EXTERNAL SUBROUTINES: C C DROT from BLAS C DGIV from SLICOT. C C Local variables. C DOUBLE PRECISION SC, SS INTEGER I, II, J, JJ, JJPVT, IFICA1, IFIRA1, MNI, NELM C IFIRA1 = IFIRA - 1 IFICA1 = IFICA - 1 C DO 20 I = NRA, 1, -1 II = IFIRA1 + I MNI = NCA - NRA + I JJPVT = IFICA1 + MNI NELM = IFIRA1 + I DO 10 J = MNI - 1, 1, -1 JJ = IFICA1 + J C C Determine the Givens transformation on columns jj and jjpvt. C Apply the transformation to these columns to annihilate C A(ii,jj) (from rows 1 up to ifira1+i). C Apply the transformation also to the E-matrix C (from rows 1 up to ifira1). C Update column transformation matrix Z. C CALL DGIV(A(II,JJPVT), A(II,JJ), SC, SS) CALL DROT(NELM, A(1,JJPVT), 1, A(1,JJ), 1, SC, SS) A(II,JJ) = 0.0D0 CALL DROT(IFIRA1, E(1,JJPVT), 1, E(1,JJ), 1, SC, SS) CALL DROT(N, Z(1,JJPVT), 1, Z(1,JJ), 1, SC, SS) 10 CONTINUE 20 CONTINUE C RETURN C *** Last line of TRIAAK ********************************************* END ** END OF TRIAAKTEXT *UPTODATE TRIAEKTEXT SUBROUTINE TRIAEK(A, LDA, E, Q, LDQ, M, N, NRE, NCE, IFIRE, * IFICE, IFICA) C C PURPOSE: C C Subroutine TRIAEK reduces a submatrix E(k) of E to upper triangu- C lar form by row Givens rotations only. C Here E(k) = E(IFIRE:me,IFICE:ne), where me = IFIRE - 1 + NRE, C ne = IFICE - 1 + NCE. C Matrix E(k) is assumed to have full column rank on entry. Hence, C no pivoting is done during the reduction process. See Algorithm C 2.3.1 and Remark 2.3.4 in [1]. C The constructed row transformations are also applied to matrix C A(k) = A(IFIRE:me,IFICA:N). C Note that in A(k) rows are transformed with the same row indices C as in E but with column indices different from those in E. C REMARK: This routine is intended to be called only from the C SLICOT auxiliary routine TRIRED. C C PARAMETERS: C C A - DOUBLE PRECISION array of dimension (LDA,N). C On entry, it contains the submatrix A(k). C On return, it contains the transformed matrix A(k). C E - DOUBLE PRECISION array of dimension (LDA,N). C On entry, it contains the submatrix E(k) of full column C rank to be reduced to upper triangular form. C On return, it contains the transformed matrix E(k). C Q - DOUBLE PRECISION array of dimension (LDQ,M). C On entry, Q contains the row transformations corresponding C to the input matrices A and E. C On return, it contains the product of the input matrix Q and C the row transformation matrix that has transformed the rows C of the matrices A and E. C M - INTEGER. C Number of rows of A and E. (1 <= M <= LDA). C N - INTEGER. C Number of columns of A and E. (N >= 1). C NRE - INTEGER C Number of rows in E to be transformed (1 <= NRE <= M). C NCE - INTEGER. C Number of columns in E to be transformed (1 <= NCE <= N). C IFIRE - INTEGER. C Index of first row in E to be transformed. C IFICE - INTEGER. C Index of first column in E to be transformed. C IFICA - INTEGER. C Index of first column in A to be transformed. C C REFERENCES: C C [1] Th.G.J. Beelen, New Algorithms for Computing the Kronecker C structure of a Pencil with Applications to Systems and C Control Theory, Ph.D.Thesis, Eindhoven University of C Technology, The Netherlands, 1987. C C CONTRIBUTOR: Th.G.J. Beelen (Philips Glas Eindhoven) C C REVISIONS: 1988, January 29. C C Specification of the parameters. C C .. Scalar arguments .. C INTEGER LDA, LDQ, M, N, NRE, NCE, IFIRE, IFICE, IFICA C C .. Array arguments .. C DOUBLE PRECISION A(LDA,N), E(LDA,N), Q(LDQ,M) C C EXTERNAL SUBROUTINES: C C DROT from BLAS C DGIV from SLICOT. C C Local variables. C DOUBLE PRECISION SC, SS INTEGER I, II, IIPVT, J, JJ, IFICE1, IFIRE1, NELM C IFIRE1 = IFIRE - 1 IFICE1 = IFICE - 1 C DO 20 J = 1, NCE JJ = IFICE1 + J IIPVT = IFIRE1 + J DO 10 I = J + 1, NRE II = IFIRE1 + I C C Determine the Givens transformation on rows ii and iipvt. C Apply the transformation to these rows (in whole E-matrix) C to annihilate E(ii,jj) (from columns jj up to n). C Apply the transformations also to the A-matrix C (from columns ifica up to n). C Update the row transformation matrix Q. C CALL DGIV(E(IIPVT,JJ), E(II,JJ), SC, SS) NELM = N - JJ + 1 CALL DROT(NELM, E(IIPVT,JJ), LDA, E(II,JJ), LDA, SC, SS) E(II,JJ) = 0.0D0 NELM = N - IFICA + 1 CALL DROT(NELM, A(IIPVT,IFICA), LDA, A(II,IFICA), LDA, * SC, SS) CALL DROT(M, Q(IIPVT,1), LDQ, Q(II,1), LDQ, SC, SS) 10 CONTINUE 20 CONTINUE C RETURN C *** Last line of TRIAEK ********************************************* END ** END OF TRIAEKTEXT *UPTODATE TRIREDTEXT SUBROUTINE TRIRED(A, LDA, E, Q, LDQ, Z, LDZ, M, N, NBLCKS, * INUK, IMUK, IERR) C C PURPOSE: C C On entry, it is assumed that the M by N matrices A and E have C been transformed to generalized Schur form by unitary trans- C formations (see Algorithm 3.2.1 in [1]), i.e., C C | s*E(eps,inf)-A(eps,inf) | X | C s*E - A = |-------------------------|-------------| . C | 0 | s*E(r)-A(r) | C C Here the pencil s*E(eps,inf)-A(eps,inf) is in staircase form. C This pencil contains all Kronecker column indices and infinite C elementary divisors of the pencil s*E - A. C The pencil s*E(r)-A(r) contains all Kronecker row indices and C finite elementary divisors of s*E - A. C Subroutine TRIRED performs the triangularization of the sub- C matrices having full row and column rank in the pencil C s*E(eps,inf)-A(eps,inf) using Algorithm 3.3.1 in [1]. C REMARK: This routine is intended to be called only from the C SLICOT routine FSTAIR. C C PARAMETERS: C C A - DOUBLE PRECISION array of dimension (LDA,N). C On entry, it contains the matrix A to be reduced. C On return, it contains the transformed matrix A. C E - DOUBLE PRECISION array of dimension (LDA,N). C On entry, it contains the matrix E to be reduced. C On return, it contains the transformed matrix E. C Q - DOUBLE PRECISION array of dimension (LDQ,M). C On entry, Q contains the row transformations corresponding C to the input matrices A and E. C On return, it contains the product of the input matrix Q and C the row transformation matrix that has transformed the rows C of the matrices A and E. C Z - DOUBLE PRECISION array of dimension (LDZ,N). C On entry, Z contains the column transformations corresponding C to the input matrices A and E. C On return, it contains the product of the input matrix Z and C the column transformation matrix that has transformed the C columns of the matrices A and E. C M - INTEGER. C Number of rows in A and E, 1 <= M <= LDA. C N - INTEGER. C Number of columns in A and E, N >= 1. C NBLCKS - INTEGER. C Number of submatrices having full row rank >=0 in A(eps,inf). C INUK - INTEGER array of dimension (NBLCKS). C Array containing the row dimensions nu(k) (k=1,..,NBLCKS) C of the submatrices having full row rank in the pencil C s*E(eps,inf)-A(eps,inf). C IMUK - INTEGER array of dimension (NBLCKS). C Array containing the column dimensions mu(k) (k=1,..,NBLCKS) C of the submatrices having full column rank in the pencil. C IERR - INTEGER. C IERR = 0, successful completion, C 1, incorrect dimensions of a full row rank submatrix, C 2, incorrect dimensions of a full column rank submatrix C C REFERENCES: C C [1] Th.G.J. Beelen, New Algorithms for Computing the Kronecker C structure of a Pencil with Applications to Systems and C Control Theory, Ph.D.Thesis, Eindhoven University of C Technology, The Netherlands, 1987. C C CONTRIBUTOR: Th.G.J. Beelen (Philips Glas Eindhoven) C C REVISIONS: 1988, January 29. C C Specification of the parameters. C C .. Scalar arguments .. C INTEGER LDA, LDQ, LDZ, M, N, NBLCKS, IERR C C .. Array arguments .. C DOUBLE PRECISION A(LDA,N), E(LDA,N), Q(LDQ,M), Z(LDZ,N) INTEGER INUK(NBLCKS), IMUK(NBLCKS) C C EXTERNAL SUBROUTINES: C C TRIAAK, TRIAEK from SLICOT. C C Local variables. C INTEGER ISMUK, ISNUK1, IFIRA, IFICA, IFIRE, IFICE INTEGER I, K, MUK, MUKP1, NUK C C ISMUK = sum(i=1,...,k) MU(i), C ISNUK1 = sum(i=1,...,k-1) NU(i). C ISMUK = 0 ISNUK1 = 0 DO 10 I = 1, NBLCKS ISMUK = ISMUK + IMUK(I) ISNUK1 = ISNUK1 + INUK(I) 10 CONTINUE C C NOTE: ISNUK1 has not yet the correct value. C IERR = 0 MUKP1 = 0 DO 20 K = NBLCKS, 1, -1 MUK = IMUK(K) NUK = INUK(K) ISNUK1 = ISNUK1 - NUK C C Determine left upper absolute coordinates of E(k) in E-matrix. C IFIRE = 1 + ISNUK1 IFICE = 1 + ISMUK C C Determine left upper absolute coordinates of A(k) in A-matrix. C IFIRA = IFIRE IFICA = IFICE - MUK C C Reduce E(k) to upper triangular form using Givens C transformations on rows only. Apply the same transformations C to the rows of A(k). C IF (MUKP1 .GT. NUK) THEN IERR = 1 RETURN END IF C CALL TRIAEK(A, LDA, E, Q, LDQ, M, N, NUK, MUKP1, IFIRE, IFICE, * IFICA) C C Reduce A(k) to upper triangular form using Givens C transformations on columns only. Apply the same transformations C to the columns in the E-matrix. C IF (NUK .GT. MUK) THEN IERR = 2 RETURN END IF C CALL TRIAAK(A, LDA, E, Z, LDZ, N, NUK, MUK, IFIRA, IFICA) C ISMUK = ISMUK - MUK MUKP1 = MUK 20 CONTINUE C RETURN C *** Last line of TRIRED ********************************************* END SUBROUTINE BAE(A, LDA, E, Q, LDQ, Z, LDZ, M, N, ISTAIR, IFIRA, * IFICA, NCA, RANK, WRK, IWRK, TOL) C C LIBRARY INDEX: C C C C PURPOSE: C C Let A and E be M x N matrices with E in column echelon form. C Let AA and EE be the following submatrices of A and E: C AA := A(IFIRA : M ; IFICA : N) C EE := E(IFIRA : M ; IFICA : N). C Let Aj and Ej be the following submatrices of AA and EE: C Aj := A(IFIRA : M ; IFICA : IFICA + NCA -1) and C Ej := E(IFIRA : M ; IFICA + NCA : N). C C The subroutine BAE transforms (AA,EE) such that Aj is row C compressed while keeping matrix Ej in column echelon form C (which may be different from the form on entry). C In fact BAE performs the j-th step of Algorithm 3.2.1 in [1]. C Furthermore, BAE determines the rank RANK of the submatrix Ej. C This is equal to the number of corner points in submatrix Ej. C REMARK: This routine is intended to be called only from the C SLICOT routine FSTAIR. C C PARAMETERS: C C A - DOUBLE PRECISION array of DIMENSION (LDA,N). C On entry, A(IFIRA : M ; IFICA : IFICA + NCA -1) contains the C matrix AA. C On return, it contains the matrix AA that has been row com- C pressed while keeping EE in column echelon form. C LDA - INTEGER. C LDA is the leading dimension of the arrays A and E. LDA >= M. C E - DOUBLE PRECISION array of DIMENSION (LDA,N). C On entry, E(IFIRA : M ; IFICA + NCA : N) contains the matrix C EE which is in column echelon form. C On return, it contains the transformed matrix EE which is kept C in column echelon form. C Q - DOUBLE PRECISION array of DIMENSION (LDQ,M). C On entry, the array Q contains the row transformations C corresponding to the input matrices A and E. C On return, it contains the M x M unitary matrix Q which is the C product of the input matrix Q and the row transformation C matrix that has transformed the rows of the matrices A and E. C LDQ - INTEGER. C LDQ is the leading dimension of the matrix Q. LDQ >= M. C Z - DOUBLE PRECISION array of DIMENSION (LDZ,N). C On entry, the array Z contains the column transformations C corresponding to the input matrices A and E. C On return, it contains the N x N unitary matrix Z which is the C product of the input matrix Z and the column transformation C matrix that has transformed the columns of A and E. C LDZ - INTEGER. C LDZ is the leading dimension of the matrix Z. LDZ >= N. C M - INTEGER. C M is the number of rows of the matrices A, E and Q. M >= 1. C N - INTEGER. C N is the number of columns of the matrices A, E and Z. N >= 1. C ISTAIR - INTEGER array of DIMENSION (M). C On entry, ISTAIR contains information on the column echelon C form of the input matrix E as follows: C ISTAIR(i) = + j: the boundary element E(i,j) is a corner point C - j: the boundary element E(i,j) is not a corner C point. C (i=1,...,M) C On return, ISTAIR contains the same information for the trans- C formed matrix E. C IFIRA - INTEGER. C IFIRA is the first row index of the submatrix Aj and Ej in C matrix A and E, respectively. C IFICA - INTEGER. C IFICA and IFICA + NCA are the first column index of the C submatrices Aj and Ej in the matrices A and E, respectively. C NCA - INTEGER. C NCA is the number of columns of the submatrix Aj in A. C RANK - INTEGER. C Numerical rank of the submatrix Ej in E (based on TOL). C WRK - DOUBLE PRECISION array of DIMENSION (N). C A real work space array. C IWRK - INTEGER array of DIMENSION (N). C An integer work space array. C TOL - DOUBLE PECISION. C TOL is the tolerance used when considering matrix elements to C be zero. TOL should be set to RE * max(||A||,||E||) + AE, C where ||.|| is the Frobenius norm. AE and RE are the absolute C and relative accuracy respectively. C A recommanded choice is AE = EPS and RE = 100*EPS, where EPS C is the machine precision. C C REFERENCES: C C [1] Th.G.J. Beelen, New Algorithms for Computing the Kronecker C structure of a Pencil with Applications to Systems and C Control Theory, Ph.D.Thesis, Eindhoven University of C Technology, The Netherlands, 1987. C C CONTRIBUTOR: Th.G.J. Beelen (Philips Glass Eindhoven). C C REVISIONS: 1988, January 29. C C Specification of the parameters. C C .. Scalar arguments .. C INTEGER LDA, LDQ, LDZ, M, N, IFIRA, IFICA, NCA, RANK DOUBLE PRECISION TOL C C .. Array arguments .. C INTEGER ISTAIR(M), IWRK(N) DOUBLE PRECISION A(LDA,N), E(LDA,N), Q(LDQ,M), Z(LDZ,N), WRK(N) C C EXTERNAL SUBROUTINES/FUNCTIONS: C C IDAMAX, DROT, DSWAP from BLAS. C DGIV from SLICOT. C C Local variables. C INTEGER I, II, IMX, IP, IR, IST1, IST2, ISTPVT, ITYPE, * IFIRA1, IFICA1, JPVT, JC1, JC2, NROWS, * K, K1, KK, L, LSAV, LL, MK1, MXRANK, NELM, MJ, NJ DOUBLE PRECISION BMXNRM, BMX, SC, SS, EIJPVT LOGICAL LZERO C C Initialisation. C C NJ = number of columns in submatrix Aj, C MJ = number of rows in submatrices Aj and Ej. C NJ = NCA MJ = M + 1 - IFIRA IFIRA1 = IFIRA - 1 IFICA1 = IFICA - 1 DO 10 I = 1, NJ IWRK(I) = I 10 CONTINUE K = 1 LZERO = .FALSE. RANK = MIN0(NJ,MJ) MXRANK = RANK C C WHILE (K <= MXRANK) and (LZERO = FALSE) DO 20 IF ((K .LE. MXRANK) .AND. (.NOT.LZERO)) THEN C C Determine column in Aj with largest max-norm. C BMXNRM = 0.0D0 LSAV = K DO 30 L = K, NJ C C IMX is relative index in column L of Aj where max el. is C found. C NOTE: the first el. in column L is in row K of matrix Aj. C KK = IFIRA1 + K LL = IFICA1 + L IMX = IDAMAX(MJ - K + 1, A(KK,LL), 1) BMX = DABS(A(IMX + KK - 1, LL)) IF (BMX .GT. BMXNRM) THEN BMXNRM = BMX LSAV = L END IF 30 CONTINUE C IF (BMXNRM .LT. TOL) THEN C C Set submatrix of Aj to zero. C DO 50 L = K, NJ LL = IFICA1 + L DO 40 I = K, MJ II = IFIRA1 + I A(II,LL) = 0.0D0 40 CONTINUE 50 CONTINUE LZERO = .TRUE. RANK = K - 1 ELSE C C Check whether columns have to be interchanged. C IF (LSAV .NE. K) THEN C C Interchange the columns in A which correspond to the C columns lsav and k in Aj. Store the permutation in IWRK. C CALL DSWAP(M, A(1,IFICA1 + K), 1, A(1,IFICA1 + LSAV), 1) IP = IWRK(LSAV) IWRK(LSAV) = IWRK(K) IWRK(K) = IP END IF C K1 = K + 1 MK1 = NJ - K + 1 + (N - NCA - IFICA1) KK = IFICA1 + K C DO 90 IR = K1, MJ C I = MJ - IR + K1 C C II = absolute row number in A corresponding to row i in C Aj. C II = IFIRA1 + I C C Construct Givens transformation to annihilate Aj(i,k). C Apply the row transformation to whole matrix A. C (NOT only to Aj). C Update row transformation matrix Q. C CALL DGIV(A(II - 1,KK), A(II,KK), SC, SS) CALL DROT(MK1, A(II - 1,KK), LDA, A(II,KK), LDA, SC, SS) A(II,KK) = 0.0D0 CALL DROT(M, Q(II - 1,1), LDQ, Q(II,1), LDQ, SC, SS) C C Determine boundary type of matrix E at rows II-1 and II. C IST1 = ISTAIR(II - 1) IST2 = ISTAIR(II) IF ((IST1 * IST2) .GT. 0) THEN IF (IST1 .GT. 0) THEN C C boundary form = (* x) C (0 *) C ITYPE = 1 ELSE C C boundary form = (x x) C (x x) C ITYPE = 3 END IF ELSE IF (IST1 .LT. 0) THEN C C boundary form = (x x) C (* x) C ITYPE=2 ELSE C C boundary form = (* x) C (0 x) C ITYPE = 4 END IF END IF C C Apply row transformation also to matrix E. C C JC1 = absolute number of the column in E in which stair C element of row i-1 of Ej is present. C JC2 = absolute number of the column in E in which stair C element of row i of Ej is present. C C NOTE: JC1 < JC2 if ITYPE = 1. C JC1 = JC2 if ITYPE = 2, 3 or 4. C JC1 = IABS(IST1) JC2 = IABS(IST2) JPVT = MIN0(JC1,JC2) NELM = N - JPVT + 1 C CALL DROT(NELM, E(II-1,JPVT), LDA, E(II,JPVT), LDA, * SC, SS) EIJPVT = E(II,JPVT) C GOTO (80, 60, 90, 70), ITYPE C 60 IF (DABS(EIJPVT) .LT. TOL) THEN C (x x) (* x) C Boundary form has been changed from (* x) to (0 x) C ISTPVT = ISTAIR(II) ISTAIR(II - 1) = ISTPVT ISTAIR(II) = -(ISTPVT + 1) E(II, JPVT) = 0.0D0 END IF GOTO 90 C 70 IF (DABS(EIJPVT) .GE. TOL) THEN C C (* x) (x x) C Boundary form has been changed from (0 x) to (* x) C ISTPVT = ISTAIR(II - 1) ISTAIR(II - 1) = -ISTPVT ISTAIR(II) = ISTPVT END IF GOTO 90 C C Construct column Givens transformation to annihilate C E(ii,jpvt). C Apply column Givens transformation to matrix E. C (NOT only to Ej). C 80 CALL DGIV(E(II,JPVT + 1), E(II,JPVT), SC, SS) CALL DROT(II, E(1,JPVT + 1), 1, E(1,JPVT), 1, SC, SS) E(II,JPVT) = 0.0D0 C C Apply this transformation also to matrix A. C (NOT only to Aj). C Update column transformation matrix Z. C CALL DROT(M, A(1,JPVT + 1), 1, A(1,JPVT), 1, SC, SS) CALL DROT(N, Z(1,JPVT + 1), 1, Z(1,JPVT), 1, SC, SS) C 90 CONTINUE C K = K + 1 END IF GOTO 20 END IF C END WHILE 20 C C Permute columns of Aj to original order. C NROWS = IFIRA1 + RANK DO 120 I = 1, NROWS DO 100 K = 1, NJ KK = IFICA1 + K WRK(IWRK(K)) = A(I,KK) 100 CONTINUE DO 110 K = 1, NJ KK = IFICA1 + K A(I,KK) = WRK(K) 110 CONTINUE 120 CONTINUE C RETURN C *** Last line of BAE ************************************************ END ** END OF BAETEXT *UPTODATE DGIVTEXT SUBROUTINE DGIV(DA, DB, DC, DS) C C LIBRARY INDEX: C C 2.1.4 Decompositions and transformations. C C PURPOSE: C C This routine constructs the Givens transformation C C ( DC DS ) C G = ( ), DC**2 + DS**2 = 1.0D0 , C (-DS DC ) C T T C such that the vector (DA,DB) is transformed into (R,0), i.e., C C ( DC DS ) ( DA ) ( R ) C ( ) ( ) = ( ) C (-DS DC ) ( DB ) ( 0 ) . C C This routine is a modification of the BLAS routine DROTG C (Algorithm 539) in order to leave the arguments DA and DB C unchanged. The value or R is not returned. C C CONTRIBUTOR: P. Van Dooren (PRLB). C C REVISIONS: 1987, November 24. C C Specification of parameters. C C .. Scalar Arguments .. C DOUBLE PRECISION DA, DB, DC, DS C C Local variables. C DOUBLE PRECISION R, U, V C IF (DABS(DA) .GT. DABS(DB)) THEN U = DA + DA V = DB/U R = DSQRT(0.25D0 + V**2) * U DC = DA/R DS = V * (DC + DC) ELSE IF (DB .NE. 0.0D0) THEN U = DB + DB V = DA/U R = DSQRT(0.25D0 + V**2) * U DS = DB/R DC = V * (DS + DS) ELSE DC = 1.0D0 DS = 0.0D0 END IF END IF RETURN C *** Last line of DGIV *********************************************** END ** END OF DGIVTEXT *UPTODATE DROTITEXT SUBROUTINE DROTI (N, X, INCX, Y, INCY, C, S) C C LIBRARY INDEX: C C 2.1.4 Decompositions and transfromations. C C PURPOSE: C C The subroutine DROTI performs the Givens transformation, defined C by C (cos) and S (sin), and interchanges the vectors involved, C i.e., C C |X(i)| | 0 1 | | C S | |X(i)| C | | := | | x | | x | |, i = 1,...N. C |Y(i)| | 1 0 | |-S C | |Y(i)| C C REMARK. This routine is a modification of DROT from BLAS. C C CONTRIBUTOR: Th.G.J. Beelen (Philips Glass Eindhoven) C C REVISIONS: 1988, February 02. C C Specification of the parameters. C C .. Scalar argumants .. C INTEGER INCX, INCY, N DOUBLE PRECISION C, S C C .. Array arguments .. C DOUBLE PRECISION X(*), Y(*) C C Local variables. C DOUBLE PRECISION DTEMP INTEGER I, IX, IY C IF (N .LE. 0) RETURN IF ((INCX.NE.1) .OR. (INCY.NE.1)) THEN C C Code for unequal increments or equal increments not equal to 1. C IX = 1 IY = 1 IF (INCX.LT.0) IX = (-N+1) * INCX + 1 IF (INCY.LT.0) IY = (-N+1) * INCY + 1 DO 10 I = 1, N DTEMP = C * Y(IY) - S * X(IX) Y(IY) = C * X(IX) + S * Y(IY) X(IX) = DTEMP IX = IX + INCX IY = IY + INCY 10 CONTINUE ELSE C C Code for both increments equal to 1. C DO 20 I = 1, N DTEMP = C * Y(I) - S * X(I) Y(I) = C * X(I) + S * Y(I) X(I) = DTEMP 20 CONTINUE END IF RETURN C *** Last line if DROTI ********************************************** END