LAPACK 3.3.0

cgghrd.f

Go to the documentation of this file.
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
 All Files Functions