Return to zgghrd.f CVS log

Up to [local] / rpl / lapack / lapack

Patches pour OS/2

    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