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