#include "blaswrap.h" #include "f2c.h" /* Subroutine */ int cgghrd_(char *compq, char *compz, integer *n, integer * ilo, integer *ihi, complex *a, integer *lda, complex *b, integer *ldb, complex *q, integer *ldq, complex *z__, integer *ldz, integer *info) { /* -- LAPACK routine (version 3.1) -- Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. November 2006 Purpose ======= CGGHRD reduces a pair of complex matrices (A,B) to generalized upper Hessenberg form using unitary transformations, where A is a general matrix and B is upper triangular. The form of the generalized eigenvalue problem is A*x = lambda*B*x, and B is typically made upper triangular by computing its QR factorization and moving the unitary matrix Q to the left side of the equation. This subroutine simultaneously reduces A to a Hessenberg matrix H: Q**H*A*Z = H and transforms B to another upper triangular matrix T: Q**H*B*Z = T in order to reduce the problem to its standard form H*y = lambda*T*y where y = Z**H*x. The unitary matrices Q and Z are determined as products of Givens rotations. They may either be formed explicitly, or they may be postmultiplied into input matrices Q1 and Z1, so that Q1 * A * Z1**H = (Q1*Q) * H * (Z1*Z)**H Q1 * B * Z1**H = (Q1*Q) * T * (Z1*Z)**H If Q1 is the unitary matrix from the QR factorization of B in the original equation A*x = lambda*B*x, then CGGHRD reduces the original problem to generalized Hessenberg form. Arguments ========= COMPQ (input) CHARACTER*1 = 'N': do not compute Q; = 'I': Q is initialized to the unit matrix, and the unitary matrix Q is returned; = 'V': Q must contain a unitary matrix Q1 on entry, and the product Q1*Q is returned. COMPZ (input) CHARACTER*1 = 'N': do not compute Q; = 'I': Q is initialized to the unit matrix, and the unitary matrix Q is returned; = 'V': Q must contain a unitary matrix Q1 on entry, and the product Q1*Q is returned. N (input) INTEGER The order of the matrices A and B. N >= 0. ILO (input) INTEGER IHI (input) INTEGER ILO and IHI mark the rows and columns of A which are to be reduced. It is assumed that A is already upper triangular in rows and columns 1:ILO-1 and IHI+1:N. ILO and IHI are normally set by a previous call to CGGBAL; otherwise they should be set to 1 and N respectively. 1 <= ILO <= IHI <= N, if N > 0; ILO=1 and IHI=0, if N=0. A (input/output) COMPLEX array, dimension (LDA, N) On entry, the N-by-N general matrix to be reduced. On exit, the upper triangle and the first subdiagonal of A are overwritten with the upper Hessenberg matrix H, and the rest is set to zero. LDA (input) INTEGER The leading dimension of the array A. LDA >= max(1,N). B (input/output) COMPLEX array, dimension (LDB, N) On entry, the N-by-N upper triangular matrix B. On exit, the upper triangular matrix T = Q**H B Z. The elements below the diagonal are set to zero. LDB (input) INTEGER The leading dimension of the array B. LDB >= max(1,N). Q (input/output) COMPLEX array, dimension (LDQ, N) On entry, if COMPQ = 'V', the unitary matrix Q1, typically from the QR factorization of B. On exit, if COMPQ='I', the unitary matrix Q, and if COMPQ = 'V', the product Q1*Q. Not referenced if COMPQ='N'. LDQ (input) INTEGER The leading dimension of the array Q. LDQ >= N if COMPQ='V' or 'I'; LDQ >= 1 otherwise. Z (input/output) COMPLEX array, dimension (LDZ, N) On entry, if COMPZ = 'V', the unitary matrix Z1. On exit, if COMPZ='I', the unitary matrix Z, and if COMPZ = 'V', the product Z1*Z. Not referenced if COMPZ='N'. LDZ (input) INTEGER The leading dimension of the array Z. LDZ >= N if COMPZ='V' or 'I'; LDZ >= 1 otherwise. INFO (output) INTEGER = 0: successful exit. < 0: if INFO = -i, the i-th argument had an illegal value. Further Details =============== This routine reduces A to Hessenberg and B to triangular form by an unblocked reduction, as described in _Matrix_Computations_, by Golub and van Loan (Johns Hopkins Press). ===================================================================== Decode COMPQ Parameter adjustments */ /* Table of constant values */ static complex c_b1 = {1.f,0.f}; static complex c_b2 = {0.f,0.f}; static integer c__1 = 1; /* System generated locals */ integer a_dim1, a_offset, b_dim1, b_offset, q_dim1, q_offset, z_dim1, z_offset, i__1, i__2, i__3; complex q__1; /* Builtin functions */ void r_cnjg(complex *, complex *); /* Local variables */ static real c__; static complex s; static logical ilq, ilz; static integer jcol; extern /* Subroutine */ int crot_(integer *, complex *, integer *, complex *, integer *, real *, complex *); static integer jrow; extern logical lsame_(char *, char *); static complex ctemp; extern /* Subroutine */ int claset_(char *, integer *, integer *, complex *, complex *, complex *, integer *), clartg_(complex *, complex *, real *, complex *, complex *), xerbla_(char *, integer *); static integer icompq, icompz; a_dim1 = *lda; a_offset = 1 + a_dim1; a -= a_offset; b_dim1 = *ldb; b_offset = 1 + b_dim1; b -= b_offset; q_dim1 = *ldq; q_offset = 1 + q_dim1; q -= q_offset; z_dim1 = *ldz; z_offset = 1 + z_dim1; z__ -= z_offset; /* Function Body */ if (lsame_(compq, "N")) { ilq = FALSE_; icompq = 1; } else if (lsame_(compq, "V")) { ilq = TRUE_; icompq = 2; } else if (lsame_(compq, "I")) { ilq = TRUE_; icompq = 3; } else { icompq = 0; } /* Decode COMPZ */ if (lsame_(compz, "N")) { ilz = FALSE_; icompz = 1; } else if (lsame_(compz, "V")) { ilz = TRUE_; icompz = 2; } else if (lsame_(compz, "I")) { ilz = TRUE_; icompz = 3; } else { icompz = 0; } /* Test the input parameters. */ *info = 0; if (icompq <= 0) { *info = -1; } else if (icompz <= 0) { *info = -2; } else if (*n < 0) { *info = -3; } else if (*ilo < 1) { *info = -4; } else if (*ihi > *n || *ihi < *ilo - 1) { *info = -5; } else if (*lda < max(1,*n)) { *info = -7; } else if (*ldb < max(1,*n)) { *info = -9; } else if (ilq && *ldq < *n || *ldq < 1) { *info = -11; } else if (ilz && *ldz < *n || *ldz < 1) { *info = -13; } if (*info != 0) { i__1 = -(*info); xerbla_("CGGHRD", &i__1); return 0; } /* Initialize Q and Z if desired. */ if (icompq == 3) { claset_("Full", n, n, &c_b2, &c_b1, &q[q_offset], ldq); } if (icompz == 3) { claset_("Full", n, n, &c_b2, &c_b1, &z__[z_offset], ldz); } /* Quick return if possible */ if (*n <= 1) { return 0; } /* Zero out lower triangle of B */ i__1 = *n - 1; for (jcol = 1; jcol <= i__1; ++jcol) { i__2 = *n; for (jrow = jcol + 1; jrow <= i__2; ++jrow) { i__3 = jrow + jcol * b_dim1; b[i__3].r = 0.f, b[i__3].i = 0.f; /* L10: */ } /* L20: */ } /* Reduce A and B */ i__1 = *ihi - 2; for (jcol = *ilo; jcol <= i__1; ++jcol) { i__2 = jcol + 2; for (jrow = *ihi; jrow >= i__2; --jrow) { /* Step 1: rotate rows JROW-1, JROW to kill A(JROW,JCOL) */ i__3 = jrow - 1 + jcol * a_dim1; ctemp.r = a[i__3].r, ctemp.i = a[i__3].i; clartg_(&ctemp, &a[jrow + jcol * a_dim1], &c__, &s, &a[jrow - 1 + jcol * a_dim1]); i__3 = jrow + jcol * a_dim1; a[i__3].r = 0.f, a[i__3].i = 0.f; i__3 = *n - jcol; crot_(&i__3, &a[jrow - 1 + (jcol + 1) * a_dim1], lda, &a[jrow + ( jcol + 1) * a_dim1], lda, &c__, &s); i__3 = *n + 2 - jrow; crot_(&i__3, &b[jrow - 1 + (jrow - 1) * b_dim1], ldb, &b[jrow + ( jrow - 1) * b_dim1], ldb, &c__, &s); if (ilq) { r_cnjg(&q__1, &s); crot_(n, &q[(jrow - 1) * q_dim1 + 1], &c__1, &q[jrow * q_dim1 + 1], &c__1, &c__, &q__1); } /* Step 2: rotate columns JROW, JROW-1 to kill B(JROW,JROW-1) */ i__3 = jrow + jrow * b_dim1; ctemp.r = b[i__3].r, ctemp.i = b[i__3].i; clartg_(&ctemp, &b[jrow + (jrow - 1) * b_dim1], &c__, &s, &b[jrow + jrow * b_dim1]); i__3 = jrow + (jrow - 1) * b_dim1; b[i__3].r = 0.f, b[i__3].i = 0.f; crot_(ihi, &a[jrow * a_dim1 + 1], &c__1, &a[(jrow - 1) * a_dim1 + 1], &c__1, &c__, &s); i__3 = jrow - 1; crot_(&i__3, &b[jrow * b_dim1 + 1], &c__1, &b[(jrow - 1) * b_dim1 + 1], &c__1, &c__, &s); if (ilz) { crot_(n, &z__[jrow * z_dim1 + 1], &c__1, &z__[(jrow - 1) * z_dim1 + 1], &c__1, &c__, &s); } /* L30: */ } /* L40: */ } return 0; /* End of CGGHRD */ } /* cgghrd_ */