00001 SUBROUTINE CGGHRD( COMPQ, COMPZ, N, ILO, IHI, A, LDA, B, LDB, Q, 00002 $ LDQ, Z, LDZ, INFO ) 00003 * 00004 * -- LAPACK routine (version 3.2) -- 00005 * -- LAPACK is a software package provided by Univ. of Tennessee, -- 00006 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 00007 * November 2006 00008 * 00009 * .. Scalar Arguments .. 00010 CHARACTER COMPQ, COMPZ 00011 INTEGER IHI, ILO, INFO, LDA, LDB, LDQ, LDZ, N 00012 * .. 00013 * .. Array Arguments .. 00014 COMPLEX A( LDA, * ), B( LDB, * ), Q( LDQ, * ), 00015 $ Z( LDZ, * ) 00016 * .. 00017 * 00018 * Purpose 00019 * ======= 00020 * 00021 * CGGHRD reduces a pair of complex matrices (A,B) to generalized upper 00022 * Hessenberg form using unitary transformations, where A is a 00023 * general matrix and B is upper triangular. The form of the generalized 00024 * eigenvalue problem is 00025 * A*x = lambda*B*x, 00026 * and B is typically made upper triangular by computing its QR 00027 * factorization and moving the unitary matrix Q to the left side 00028 * of the equation. 00029 * 00030 * This subroutine simultaneously reduces A to a Hessenberg matrix H: 00031 * Q**H*A*Z = H 00032 * and transforms B to another upper triangular matrix T: 00033 * Q**H*B*Z = T 00034 * in order to reduce the problem to its standard form 00035 * H*y = lambda*T*y 00036 * where y = Z**H*x. 00037 * 00038 * The unitary matrices Q and Z are determined as products of Givens 00039 * rotations. They may either be formed explicitly, or they may be 00040 * postmultiplied into input matrices Q1 and Z1, so that 00041 * Q1 * A * Z1**H = (Q1*Q) * H * (Z1*Z)**H 00042 * Q1 * B * Z1**H = (Q1*Q) * T * (Z1*Z)**H 00043 * If Q1 is the unitary matrix from the QR factorization of B in the 00044 * original equation A*x = lambda*B*x, then CGGHRD reduces the original 00045 * problem to generalized Hessenberg form. 00046 * 00047 * Arguments 00048 * ========= 00049 * 00050 * COMPQ (input) CHARACTER*1 00051 * = 'N': do not compute Q; 00052 * = 'I': Q is initialized to the unit matrix, and the 00053 * unitary matrix Q is returned; 00054 * = 'V': Q must contain a unitary matrix Q1 on entry, 00055 * and the product Q1*Q is returned. 00056 * 00057 * COMPZ (input) CHARACTER*1 00058 * = 'N': do not compute Q; 00059 * = 'I': Q is initialized to the unit matrix, and the 00060 * unitary matrix Q is returned; 00061 * = 'V': Q must contain a unitary matrix Q1 on entry, 00062 * and the product Q1*Q is returned. 00063 * 00064 * N (input) INTEGER 00065 * The order of the matrices A and B. N >= 0. 00066 * 00067 * ILO (input) INTEGER 00068 * IHI (input) INTEGER 00069 * ILO and IHI mark the rows and columns of A which are to be 00070 * reduced. It is assumed that A is already upper triangular 00071 * in rows and columns 1:ILO-1 and IHI+1:N. ILO and IHI are 00072 * normally set by a previous call to CGGBAL; otherwise they 00073 * should be set to 1 and N respectively. 00074 * 1 <= ILO <= IHI <= N, if N > 0; ILO=1 and IHI=0, if N=0. 00075 * 00076 * A (input/output) COMPLEX array, dimension (LDA, N) 00077 * On entry, the N-by-N general matrix to be reduced. 00078 * On exit, the upper triangle and the first subdiagonal of A 00079 * are overwritten with the upper Hessenberg matrix H, and the 00080 * rest is set to zero. 00081 * 00082 * LDA (input) INTEGER 00083 * The leading dimension of the array A. LDA >= max(1,N). 00084 * 00085 * B (input/output) COMPLEX array, dimension (LDB, N) 00086 * On entry, the N-by-N upper triangular matrix B. 00087 * On exit, the upper triangular matrix T = Q**H B Z. The 00088 * elements below the diagonal are set to zero. 00089 * 00090 * LDB (input) INTEGER 00091 * The leading dimension of the array B. LDB >= max(1,N). 00092 * 00093 * Q (input/output) COMPLEX array, dimension (LDQ, N) 00094 * On entry, if COMPQ = 'V', the unitary matrix Q1, typically 00095 * from the QR factorization of B. 00096 * On exit, if COMPQ='I', the unitary matrix Q, and if 00097 * COMPQ = 'V', the product Q1*Q. 00098 * Not referenced if COMPQ='N'. 00099 * 00100 * LDQ (input) INTEGER 00101 * The leading dimension of the array Q. 00102 * LDQ >= N if COMPQ='V' or 'I'; LDQ >= 1 otherwise. 00103 * 00104 * Z (input/output) COMPLEX array, dimension (LDZ, N) 00105 * On entry, if COMPZ = 'V', the unitary matrix Z1. 00106 * On exit, if COMPZ='I', the unitary matrix Z, and if 00107 * COMPZ = 'V', the product Z1*Z. 00108 * Not referenced if COMPZ='N'. 00109 * 00110 * LDZ (input) INTEGER 00111 * The leading dimension of the array Z. 00112 * LDZ >= N if COMPZ='V' or 'I'; LDZ >= 1 otherwise. 00113 * 00114 * INFO (output) INTEGER 00115 * = 0: successful exit. 00116 * < 0: if INFO = -i, the i-th argument had an illegal value. 00117 * 00118 * Further Details 00119 * =============== 00120 * 00121 * This routine reduces A to Hessenberg and B to triangular form by 00122 * an unblocked reduction, as described in _Matrix_Computations_, 00123 * by Golub and van Loan (Johns Hopkins Press). 00124 * 00125 * ===================================================================== 00126 * 00127 * .. Parameters .. 00128 COMPLEX CONE, CZERO 00129 PARAMETER ( CONE = ( 1.0E+0, 0.0E+0 ), 00130 $ CZERO = ( 0.0E+0, 0.0E+0 ) ) 00131 * .. 00132 * .. Local Scalars .. 00133 LOGICAL ILQ, ILZ 00134 INTEGER ICOMPQ, ICOMPZ, JCOL, JROW 00135 REAL C 00136 COMPLEX CTEMP, S 00137 * .. 00138 * .. External Functions .. 00139 LOGICAL LSAME 00140 EXTERNAL LSAME 00141 * .. 00142 * .. External Subroutines .. 00143 EXTERNAL CLARTG, CLASET, CROT, XERBLA 00144 * .. 00145 * .. Intrinsic Functions .. 00146 INTRINSIC CONJG, MAX 00147 * .. 00148 * .. Executable Statements .. 00149 * 00150 * Decode COMPQ 00151 * 00152 IF( LSAME( COMPQ, 'N' ) ) THEN 00153 ILQ = .FALSE. 00154 ICOMPQ = 1 00155 ELSE IF( LSAME( COMPQ, 'V' ) ) THEN 00156 ILQ = .TRUE. 00157 ICOMPQ = 2 00158 ELSE IF( LSAME( COMPQ, 'I' ) ) THEN 00159 ILQ = .TRUE. 00160 ICOMPQ = 3 00161 ELSE 00162 ICOMPQ = 0 00163 END IF 00164 * 00165 * Decode COMPZ 00166 * 00167 IF( LSAME( COMPZ, 'N' ) ) THEN 00168 ILZ = .FALSE. 00169 ICOMPZ = 1 00170 ELSE IF( LSAME( COMPZ, 'V' ) ) THEN 00171 ILZ = .TRUE. 00172 ICOMPZ = 2 00173 ELSE IF( LSAME( COMPZ, 'I' ) ) THEN 00174 ILZ = .TRUE. 00175 ICOMPZ = 3 00176 ELSE 00177 ICOMPZ = 0 00178 END IF 00179 * 00180 * Test the input parameters. 00181 * 00182 INFO = 0 00183 IF( ICOMPQ.LE.0 ) THEN 00184 INFO = -1 00185 ELSE IF( ICOMPZ.LE.0 ) THEN 00186 INFO = -2 00187 ELSE IF( N.LT.0 ) THEN 00188 INFO = -3 00189 ELSE IF( ILO.LT.1 ) THEN 00190 INFO = -4 00191 ELSE IF( IHI.GT.N .OR. IHI.LT.ILO-1 ) THEN 00192 INFO = -5 00193 ELSE IF( LDA.LT.MAX( 1, N ) ) THEN 00194 INFO = -7 00195 ELSE IF( LDB.LT.MAX( 1, N ) ) THEN 00196 INFO = -9 00197 ELSE IF( ( ILQ .AND. LDQ.LT.N ) .OR. LDQ.LT.1 ) THEN 00198 INFO = -11 00199 ELSE IF( ( ILZ .AND. LDZ.LT.N ) .OR. LDZ.LT.1 ) THEN 00200 INFO = -13 00201 END IF 00202 IF( INFO.NE.0 ) THEN 00203 CALL XERBLA( 'CGGHRD', -INFO ) 00204 RETURN 00205 END IF 00206 * 00207 * Initialize Q and Z if desired. 00208 * 00209 IF( ICOMPQ.EQ.3 ) 00210 $ CALL CLASET( 'Full', N, N, CZERO, CONE, Q, LDQ ) 00211 IF( ICOMPZ.EQ.3 ) 00212 $ CALL CLASET( 'Full', N, N, CZERO, CONE, Z, LDZ ) 00213 * 00214 * Quick return if possible 00215 * 00216 IF( N.LE.1 ) 00217 $ RETURN 00218 * 00219 * Zero out lower triangle of B 00220 * 00221 DO 20 JCOL = 1, N - 1 00222 DO 10 JROW = JCOL + 1, N 00223 B( JROW, JCOL ) = CZERO 00224 10 CONTINUE 00225 20 CONTINUE 00226 * 00227 * Reduce A and B 00228 * 00229 DO 40 JCOL = ILO, IHI - 2 00230 * 00231 DO 30 JROW = IHI, JCOL + 2, -1 00232 * 00233 * Step 1: rotate rows JROW-1, JROW to kill A(JROW,JCOL) 00234 * 00235 CTEMP = A( JROW-1, JCOL ) 00236 CALL CLARTG( CTEMP, A( JROW, JCOL ), C, S, 00237 $ A( JROW-1, JCOL ) ) 00238 A( JROW, JCOL ) = CZERO 00239 CALL CROT( N-JCOL, A( JROW-1, JCOL+1 ), LDA, 00240 $ A( JROW, JCOL+1 ), LDA, C, S ) 00241 CALL CROT( N+2-JROW, B( JROW-1, JROW-1 ), LDB, 00242 $ B( JROW, JROW-1 ), LDB, C, S ) 00243 IF( ILQ ) 00244 $ CALL CROT( N, Q( 1, JROW-1 ), 1, Q( 1, JROW ), 1, C, 00245 $ CONJG( S ) ) 00246 * 00247 * Step 2: rotate columns JROW, JROW-1 to kill B(JROW,JROW-1) 00248 * 00249 CTEMP = B( JROW, JROW ) 00250 CALL CLARTG( CTEMP, B( JROW, JROW-1 ), C, S, 00251 $ B( JROW, JROW ) ) 00252 B( JROW, JROW-1 ) = CZERO 00253 CALL CROT( IHI, A( 1, JROW ), 1, A( 1, JROW-1 ), 1, C, S ) 00254 CALL CROT( JROW-1, B( 1, JROW ), 1, B( 1, JROW-1 ), 1, C, 00255 $ S ) 00256 IF( ILZ ) 00257 $ CALL CROT( N, Z( 1, JROW ), 1, Z( 1, JROW-1 ), 1, C, S ) 00258 30 CONTINUE 00259 40 CONTINUE 00260 * 00261 RETURN 00262 * 00263 * End of CGGHRD 00264 * 00265 END
1.7.2