LAPACK 3.3.0

sgghrd.f

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