Annotation of rpl/lapack/lapack/zgghrd.f, revision 1.1
1.1 ! bertrand 1: SUBROUTINE ZGGHRD( COMPQ, COMPZ, N, ILO, IHI, A, LDA, B, LDB, Q,
! 2: $ LDQ, Z, LDZ, INFO )
! 3: *
! 4: * -- LAPACK routine (version 3.2) --
! 5: * -- LAPACK is a software package provided by Univ. of Tennessee, --
! 6: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
! 7: * November 2006
! 8: *
! 9: * .. Scalar Arguments ..
! 10: CHARACTER COMPQ, COMPZ
! 11: INTEGER IHI, ILO, INFO, LDA, LDB, LDQ, LDZ, N
! 12: * ..
! 13: * .. Array Arguments ..
! 14: COMPLEX*16 A( LDA, * ), B( LDB, * ), Q( LDQ, * ),
! 15: $ Z( LDZ, * )
! 16: * ..
! 17: *
! 18: * Purpose
! 19: * =======
! 20: *
! 21: * ZGGHRD reduces a pair of complex matrices (A,B) to generalized upper
! 22: * Hessenberg form using unitary transformations, where A is a
! 23: * general matrix and B is upper triangular. The form of the
! 24: * generalized eigenvalue problem is
! 25: * A*x = lambda*B*x,
! 26: * and B is typically made upper triangular by computing its QR
! 27: * factorization and moving the unitary matrix Q to the left side
! 28: * of the equation.
! 29: *
! 30: * This subroutine simultaneously reduces A to a Hessenberg matrix H:
! 31: * Q**H*A*Z = H
! 32: * and transforms B to another upper triangular matrix T:
! 33: * Q**H*B*Z = T
! 34: * in order to reduce the problem to its standard form
! 35: * H*y = lambda*T*y
! 36: * where y = Z**H*x.
! 37: *
! 38: * The unitary matrices Q and Z are determined as products of Givens
! 39: * rotations. They may either be formed explicitly, or they may be
! 40: * postmultiplied into input matrices Q1 and Z1, so that
! 41: * Q1 * A * Z1**H = (Q1*Q) * H * (Z1*Z)**H
! 42: * Q1 * B * Z1**H = (Q1*Q) * T * (Z1*Z)**H
! 43: * If Q1 is the unitary matrix from the QR factorization of B in the
! 44: * original equation A*x = lambda*B*x, then ZGGHRD reduces the original
! 45: * problem to generalized Hessenberg form.
! 46: *
! 47: * Arguments
! 48: * =========
! 49: *
! 50: * COMPQ (input) CHARACTER*1
! 51: * = 'N': do not compute Q;
! 52: * = 'I': Q is initialized to the unit matrix, and the
! 53: * unitary matrix Q is returned;
! 54: * = 'V': Q must contain a unitary matrix Q1 on entry,
! 55: * and the product Q1*Q is returned.
! 56: *
! 57: * COMPZ (input) CHARACTER*1
! 58: * = 'N': do not compute Q;
! 59: * = 'I': Q is initialized to the unit matrix, and the
! 60: * unitary matrix Q is returned;
! 61: * = 'V': Q must contain a unitary matrix Q1 on entry,
! 62: * and the product Q1*Q is returned.
! 63: *
! 64: * N (input) INTEGER
! 65: * The order of the matrices A and B. N >= 0.
! 66: *
! 67: * ILO (input) INTEGER
! 68: * IHI (input) INTEGER
! 69: * ILO and IHI mark the rows and columns of A which are to be
! 70: * reduced. It is assumed that A is already upper triangular
! 71: * in rows and columns 1:ILO-1 and IHI+1:N. ILO and IHI are
! 72: * normally set by a previous call to ZGGBAL; otherwise they
! 73: * should be set to 1 and N respectively.
! 74: * 1 <= ILO <= IHI <= N, if N > 0; ILO=1 and IHI=0, if N=0.
! 75: *
! 76: * A (input/output) COMPLEX*16 array, dimension (LDA, N)
! 77: * On entry, the N-by-N general matrix to be reduced.
! 78: * On exit, the upper triangle and the first subdiagonal of A
! 79: * are overwritten with the upper Hessenberg matrix H, and the
! 80: * rest is set to zero.
! 81: *
! 82: * LDA (input) INTEGER
! 83: * The leading dimension of the array A. LDA >= max(1,N).
! 84: *
! 85: * B (input/output) COMPLEX*16 array, dimension (LDB, N)
! 86: * On entry, the N-by-N upper triangular matrix B.
! 87: * On exit, the upper triangular matrix T = Q**H B Z. The
! 88: * elements below the diagonal are set to zero.
! 89: *
! 90: * LDB (input) INTEGER
! 91: * The leading dimension of the array B. LDB >= max(1,N).
! 92: *
! 93: * Q (input/output) COMPLEX*16 array, dimension (LDQ, N)
! 94: * On entry, if COMPQ = 'V', the unitary matrix Q1, typically
! 95: * from the QR factorization of B.
! 96: * On exit, if COMPQ='I', the unitary matrix Q, and if
! 97: * COMPQ = 'V', the product Q1*Q.
! 98: * Not referenced if COMPQ='N'.
! 99: *
! 100: * LDQ (input) INTEGER
! 101: * The leading dimension of the array Q.
! 102: * LDQ >= N if COMPQ='V' or 'I'; LDQ >= 1 otherwise.
! 103: *
! 104: * Z (input/output) COMPLEX*16 array, dimension (LDZ, N)
! 105: * On entry, if COMPZ = 'V', the unitary matrix Z1.
! 106: * On exit, if COMPZ='I', the unitary matrix Z, and if
! 107: * COMPZ = 'V', the product Z1*Z.
! 108: * Not referenced if COMPZ='N'.
! 109: *
! 110: * LDZ (input) INTEGER
! 111: * The leading dimension of the array Z.
! 112: * LDZ >= N if COMPZ='V' or 'I'; LDZ >= 1 otherwise.
! 113: *
! 114: * INFO (output) INTEGER
! 115: * = 0: successful exit.
! 116: * < 0: if INFO = -i, the i-th argument had an illegal value.
! 117: *
! 118: * Further Details
! 119: * ===============
! 120: *
! 121: * This routine reduces A to Hessenberg and B to triangular form by
! 122: * an unblocked reduction, as described in _Matrix_Computations_,
! 123: * by Golub and van Loan (Johns Hopkins Press).
! 124: *
! 125: * =====================================================================
! 126: *
! 127: * .. Parameters ..
! 128: COMPLEX*16 CONE, CZERO
! 129: PARAMETER ( CONE = ( 1.0D+0, 0.0D+0 ),
! 130: $ CZERO = ( 0.0D+0, 0.0D+0 ) )
! 131: * ..
! 132: * .. Local Scalars ..
! 133: LOGICAL ILQ, ILZ
! 134: INTEGER ICOMPQ, ICOMPZ, JCOL, JROW
! 135: DOUBLE PRECISION C
! 136: COMPLEX*16 CTEMP, S
! 137: * ..
! 138: * .. External Functions ..
! 139: LOGICAL LSAME
! 140: EXTERNAL LSAME
! 141: * ..
! 142: * .. External Subroutines ..
! 143: EXTERNAL XERBLA, ZLARTG, ZLASET, ZROT
! 144: * ..
! 145: * .. Intrinsic Functions ..
! 146: INTRINSIC DCONJG, MAX
! 147: * ..
! 148: * .. Executable Statements ..
! 149: *
! 150: * Decode COMPQ
! 151: *
! 152: IF( LSAME( COMPQ, 'N' ) ) THEN
! 153: ILQ = .FALSE.
! 154: ICOMPQ = 1
! 155: ELSE IF( LSAME( COMPQ, 'V' ) ) THEN
! 156: ILQ = .TRUE.
! 157: ICOMPQ = 2
! 158: ELSE IF( LSAME( COMPQ, 'I' ) ) THEN
! 159: ILQ = .TRUE.
! 160: ICOMPQ = 3
! 161: ELSE
! 162: ICOMPQ = 0
! 163: END IF
! 164: *
! 165: * Decode COMPZ
! 166: *
! 167: IF( LSAME( COMPZ, 'N' ) ) THEN
! 168: ILZ = .FALSE.
! 169: ICOMPZ = 1
! 170: ELSE IF( LSAME( COMPZ, 'V' ) ) THEN
! 171: ILZ = .TRUE.
! 172: ICOMPZ = 2
! 173: ELSE IF( LSAME( COMPZ, 'I' ) ) THEN
! 174: ILZ = .TRUE.
! 175: ICOMPZ = 3
! 176: ELSE
! 177: ICOMPZ = 0
! 178: END IF
! 179: *
! 180: * Test the input parameters.
! 181: *
! 182: INFO = 0
! 183: IF( ICOMPQ.LE.0 ) THEN
! 184: INFO = -1
! 185: ELSE IF( ICOMPZ.LE.0 ) THEN
! 186: INFO = -2
! 187: ELSE IF( N.LT.0 ) THEN
! 188: INFO = -3
! 189: ELSE IF( ILO.LT.1 ) THEN
! 190: INFO = -4
! 191: ELSE IF( IHI.GT.N .OR. IHI.LT.ILO-1 ) THEN
! 192: INFO = -5
! 193: ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
! 194: INFO = -7
! 195: ELSE IF( LDB.LT.MAX( 1, N ) ) THEN
! 196: INFO = -9
! 197: ELSE IF( ( ILQ .AND. LDQ.LT.N ) .OR. LDQ.LT.1 ) THEN
! 198: INFO = -11
! 199: ELSE IF( ( ILZ .AND. LDZ.LT.N ) .OR. LDZ.LT.1 ) THEN
! 200: INFO = -13
! 201: END IF
! 202: IF( INFO.NE.0 ) THEN
! 203: CALL XERBLA( 'ZGGHRD', -INFO )
! 204: RETURN
! 205: END IF
! 206: *
! 207: * Initialize Q and Z if desired.
! 208: *
! 209: IF( ICOMPQ.EQ.3 )
! 210: $ CALL ZLASET( 'Full', N, N, CZERO, CONE, Q, LDQ )
! 211: IF( ICOMPZ.EQ.3 )
! 212: $ CALL ZLASET( 'Full', N, N, CZERO, CONE, Z, LDZ )
! 213: *
! 214: * Quick return if possible
! 215: *
! 216: IF( N.LE.1 )
! 217: $ RETURN
! 218: *
! 219: * Zero out lower triangle of B
! 220: *
! 221: DO 20 JCOL = 1, N - 1
! 222: DO 10 JROW = JCOL + 1, N
! 223: B( JROW, JCOL ) = CZERO
! 224: 10 CONTINUE
! 225: 20 CONTINUE
! 226: *
! 227: * Reduce A and B
! 228: *
! 229: DO 40 JCOL = ILO, IHI - 2
! 230: *
! 231: DO 30 JROW = IHI, JCOL + 2, -1
! 232: *
! 233: * Step 1: rotate rows JROW-1, JROW to kill A(JROW,JCOL)
! 234: *
! 235: CTEMP = A( JROW-1, JCOL )
! 236: CALL ZLARTG( CTEMP, A( JROW, JCOL ), C, S,
! 237: $ A( JROW-1, JCOL ) )
! 238: A( JROW, JCOL ) = CZERO
! 239: CALL ZROT( N-JCOL, A( JROW-1, JCOL+1 ), LDA,
! 240: $ A( JROW, JCOL+1 ), LDA, C, S )
! 241: CALL ZROT( N+2-JROW, B( JROW-1, JROW-1 ), LDB,
! 242: $ B( JROW, JROW-1 ), LDB, C, S )
! 243: IF( ILQ )
! 244: $ CALL ZROT( N, Q( 1, JROW-1 ), 1, Q( 1, JROW ), 1, C,
! 245: $ DCONJG( S ) )
! 246: *
! 247: * Step 2: rotate columns JROW, JROW-1 to kill B(JROW,JROW-1)
! 248: *
! 249: CTEMP = B( JROW, JROW )
! 250: CALL ZLARTG( CTEMP, B( JROW, JROW-1 ), C, S,
! 251: $ B( JROW, JROW ) )
! 252: B( JROW, JROW-1 ) = CZERO
! 253: CALL ZROT( IHI, A( 1, JROW ), 1, A( 1, JROW-1 ), 1, C, S )
! 254: CALL ZROT( JROW-1, B( 1, JROW ), 1, B( 1, JROW-1 ), 1, C,
! 255: $ S )
! 256: IF( ILZ )
! 257: $ CALL ZROT( N, Z( 1, JROW ), 1, Z( 1, JROW-1 ), 1, C, S )
! 258: 30 CONTINUE
! 259: 40 CONTINUE
! 260: *
! 261: RETURN
! 262: *
! 263: * End of ZGGHRD
! 264: *
! 265: END
CVSweb interface <joel.bertrand@systella.fr>
![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to zgghrd.f CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [local] / rpl / lapack / lapack