cgghrd.c
Go to the documentation of this file.
00001 /* cgghrd.f -- translated by f2c (version 20061008).
00002    You must link the resulting object file with libf2c:
00003         on Microsoft Windows system, link with libf2c.lib;
00004         on Linux or Unix systems, link with .../path/to/libf2c.a -lm
00005         or, if you install libf2c.a in a standard place, with -lf2c -lm
00006         -- in that order, at the end of the command line, as in
00007                 cc *.o -lf2c -lm
00008         Source for libf2c is in /netlib/f2c/libf2c.zip, e.g.,
00009 
00010                 http://www.netlib.org/f2c/libf2c.zip
00011 */
00012 
00013 #include "f2c.h"
00014 #include "blaswrap.h"
00015 
00016 /* Table of constant values */
00017 
00018 static complex c_b1 = {1.f,0.f};
00019 static complex c_b2 = {0.f,0.f};
00020 static integer c__1 = 1;
00021 
00022 /* Subroutine */ int cgghrd_(char *compq, char *compz, integer *n, integer *
00023         ilo, integer *ihi, complex *a, integer *lda, complex *b, integer *ldb, 
00024          complex *q, integer *ldq, complex *z__, integer *ldz, integer *info)
00025 {
00026     /* System generated locals */
00027     integer a_dim1, a_offset, b_dim1, b_offset, q_dim1, q_offset, z_dim1, 
00028             z_offset, i__1, i__2, i__3;
00029     complex q__1;
00030 
00031     /* Builtin functions */
00032     void r_cnjg(complex *, complex *);
00033 
00034     /* Local variables */
00035     real c__;
00036     complex s;
00037     logical ilq, ilz;
00038     integer jcol;
00039     extern /* Subroutine */ int crot_(integer *, complex *, integer *, 
00040             complex *, integer *, real *, complex *);
00041     integer jrow;
00042     extern logical lsame_(char *, char *);
00043     complex ctemp;
00044     extern /* Subroutine */ int claset_(char *, integer *, integer *, complex 
00045             *, complex *, complex *, integer *), clartg_(complex *, 
00046             complex *, real *, complex *, complex *), xerbla_(char *, integer 
00047             *);
00048     integer icompq, icompz;
00049 
00050 
00051 /*  -- LAPACK routine (version 3.2) -- */
00052 /*     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. */
00053 /*     November 2006 */
00054 
00055 /*     .. Scalar Arguments .. */
00056 /*     .. */
00057 /*     .. Array Arguments .. */
00058 /*     .. */
00059 
00060 /*  Purpose */
00061 /*  ======= */
00062 
00063 /*  CGGHRD reduces a pair of complex matrices (A,B) to generalized upper */
00064 /*  Hessenberg form using unitary transformations, where A is a */
00065 /*  general matrix and B is upper triangular.  The form of the generalized */
00066 /*  eigenvalue problem is */
00067 /*     A*x = lambda*B*x, */
00068 /*  and B is typically made upper triangular by computing its QR */
00069 /*  factorization and moving the unitary matrix Q to the left side */
00070 /*  of the equation. */
00071 
00072 /*  This subroutine simultaneously reduces A to a Hessenberg matrix H: */
00073 /*     Q**H*A*Z = H */
00074 /*  and transforms B to another upper triangular matrix T: */
00075 /*     Q**H*B*Z = T */
00076 /*  in order to reduce the problem to its standard form */
00077 /*     H*y = lambda*T*y */
00078 /*  where y = Z**H*x. */
00079 
00080 /*  The unitary matrices Q and Z are determined as products of Givens */
00081 /*  rotations.  They may either be formed explicitly, or they may be */
00082 /*  postmultiplied into input matrices Q1 and Z1, so that */
00083 /*       Q1 * A * Z1**H = (Q1*Q) * H * (Z1*Z)**H */
00084 /*       Q1 * B * Z1**H = (Q1*Q) * T * (Z1*Z)**H */
00085 /*  If Q1 is the unitary matrix from the QR factorization of B in the */
00086 /*  original equation A*x = lambda*B*x, then CGGHRD reduces the original */
00087 /*  problem to generalized Hessenberg form. */
00088 
00089 /*  Arguments */
00090 /*  ========= */
00091 
00092 /*  COMPQ   (input) CHARACTER*1 */
00093 /*          = 'N': do not compute Q; */
00094 /*          = 'I': Q is initialized to the unit matrix, and the */
00095 /*                 unitary matrix Q is returned; */
00096 /*          = 'V': Q must contain a unitary matrix Q1 on entry, */
00097 /*                 and the product Q1*Q is returned. */
00098 
00099 /*  COMPZ   (input) CHARACTER*1 */
00100 /*          = 'N': do not compute Q; */
00101 /*          = 'I': Q is initialized to the unit matrix, and the */
00102 /*                 unitary matrix Q is returned; */
00103 /*          = 'V': Q must contain a unitary matrix Q1 on entry, */
00104 /*                 and the product Q1*Q is returned. */
00105 
00106 /*  N       (input) INTEGER */
00107 /*          The order of the matrices A and B.  N >= 0. */
00108 
00109 /*  ILO     (input) INTEGER */
00110 /*  IHI     (input) INTEGER */
00111 /*          ILO and IHI mark the rows and columns of A which are to be */
00112 /*          reduced.  It is assumed that A is already upper triangular */
00113 /*          in rows and columns 1:ILO-1 and IHI+1:N.  ILO and IHI are */
00114 /*          normally set by a previous call to CGGBAL; otherwise they */
00115 /*          should be set to 1 and N respectively. */
00116 /*          1 <= ILO <= IHI <= N, if N > 0; ILO=1 and IHI=0, if N=0. */
00117 
00118 /*  A       (input/output) COMPLEX array, dimension (LDA, N) */
00119 /*          On entry, the N-by-N general matrix to be reduced. */
00120 /*          On exit, the upper triangle and the first subdiagonal of A */
00121 /*          are overwritten with the upper Hessenberg matrix H, and the */
00122 /*          rest is set to zero. */
00123 
00124 /*  LDA     (input) INTEGER */
00125 /*          The leading dimension of the array A.  LDA >= max(1,N). */
00126 
00127 /*  B       (input/output) COMPLEX array, dimension (LDB, N) */
00128 /*          On entry, the N-by-N upper triangular matrix B. */
00129 /*          On exit, the upper triangular matrix T = Q**H B Z.  The */
00130 /*          elements below the diagonal are set to zero. */
00131 
00132 /*  LDB     (input) INTEGER */
00133 /*          The leading dimension of the array B.  LDB >= max(1,N). */
00134 
00135 /*  Q       (input/output) COMPLEX array, dimension (LDQ, N) */
00136 /*          On entry, if COMPQ = 'V', the unitary matrix Q1, typically */
00137 /*          from the QR factorization of B. */
00138 /*          On exit, if COMPQ='I', the unitary matrix Q, and if */
00139 /*          COMPQ = 'V', the product Q1*Q. */
00140 /*          Not referenced if COMPQ='N'. */
00141 
00142 /*  LDQ     (input) INTEGER */
00143 /*          The leading dimension of the array Q. */
00144 /*          LDQ >= N if COMPQ='V' or 'I'; LDQ >= 1 otherwise. */
00145 
00146 /*  Z       (input/output) COMPLEX array, dimension (LDZ, N) */
00147 /*          On entry, if COMPZ = 'V', the unitary matrix Z1. */
00148 /*          On exit, if COMPZ='I', the unitary matrix Z, and if */
00149 /*          COMPZ = 'V', the product Z1*Z. */
00150 /*          Not referenced if COMPZ='N'. */
00151 
00152 /*  LDZ     (input) INTEGER */
00153 /*          The leading dimension of the array Z. */
00154 /*          LDZ >= N if COMPZ='V' or 'I'; LDZ >= 1 otherwise. */
00155 
00156 /*  INFO    (output) INTEGER */
00157 /*          = 0:  successful exit. */
00158 /*          < 0:  if INFO = -i, the i-th argument had an illegal value. */
00159 
00160 /*  Further Details */
00161 /*  =============== */
00162 
00163 /*  This routine reduces A to Hessenberg and B to triangular form by */
00164 /*  an unblocked reduction, as described in _Matrix_Computations_, */
00165 /*  by Golub and van Loan (Johns Hopkins Press). */
00166 
00167 /*  ===================================================================== */
00168 
00169 /*     .. Parameters .. */
00170 /*     .. */
00171 /*     .. Local Scalars .. */
00172 /*     .. */
00173 /*     .. External Functions .. */
00174 /*     .. */
00175 /*     .. External Subroutines .. */
00176 /*     .. */
00177 /*     .. Intrinsic Functions .. */
00178 /*     .. */
00179 /*     .. Executable Statements .. */
00180 
00181 /*     Decode COMPQ */
00182 
00183     /* Parameter adjustments */
00184     a_dim1 = *lda;
00185     a_offset = 1 + a_dim1;
00186     a -= a_offset;
00187     b_dim1 = *ldb;
00188     b_offset = 1 + b_dim1;
00189     b -= b_offset;
00190     q_dim1 = *ldq;
00191     q_offset = 1 + q_dim1;
00192     q -= q_offset;
00193     z_dim1 = *ldz;
00194     z_offset = 1 + z_dim1;
00195     z__ -= z_offset;
00196 
00197     /* Function Body */
00198     if (lsame_(compq, "N")) {
00199         ilq = FALSE_;
00200         icompq = 1;
00201     } else if (lsame_(compq, "V")) {
00202         ilq = TRUE_;
00203         icompq = 2;
00204     } else if (lsame_(compq, "I")) {
00205         ilq = TRUE_;
00206         icompq = 3;
00207     } else {
00208         icompq = 0;
00209     }
00210 
00211 /*     Decode COMPZ */
00212 
00213     if (lsame_(compz, "N")) {
00214         ilz = FALSE_;
00215         icompz = 1;
00216     } else if (lsame_(compz, "V")) {
00217         ilz = TRUE_;
00218         icompz = 2;
00219     } else if (lsame_(compz, "I")) {
00220         ilz = TRUE_;
00221         icompz = 3;
00222     } else {
00223         icompz = 0;
00224     }
00225 
00226 /*     Test the input parameters. */
00227 
00228     *info = 0;
00229     if (icompq <= 0) {
00230         *info = -1;
00231     } else if (icompz <= 0) {
00232         *info = -2;
00233     } else if (*n < 0) {
00234         *info = -3;
00235     } else if (*ilo < 1) {
00236         *info = -4;
00237     } else if (*ihi > *n || *ihi < *ilo - 1) {
00238         *info = -5;
00239     } else if (*lda < max(1,*n)) {
00240         *info = -7;
00241     } else if (*ldb < max(1,*n)) {
00242         *info = -9;
00243     } else if (ilq && *ldq < *n || *ldq < 1) {
00244         *info = -11;
00245     } else if (ilz && *ldz < *n || *ldz < 1) {
00246         *info = -13;
00247     }
00248     if (*info != 0) {
00249         i__1 = -(*info);
00250         xerbla_("CGGHRD", &i__1);
00251         return 0;
00252     }
00253 
00254 /*     Initialize Q and Z if desired. */
00255 
00256     if (icompq == 3) {
00257         claset_("Full", n, n, &c_b2, &c_b1, &q[q_offset], ldq);
00258     }
00259     if (icompz == 3) {
00260         claset_("Full", n, n, &c_b2, &c_b1, &z__[z_offset], ldz);
00261     }
00262 
00263 /*     Quick return if possible */
00264 
00265     if (*n <= 1) {
00266         return 0;
00267     }
00268 
00269 /*     Zero out lower triangle of B */
00270 
00271     i__1 = *n - 1;
00272     for (jcol = 1; jcol <= i__1; ++jcol) {
00273         i__2 = *n;
00274         for (jrow = jcol + 1; jrow <= i__2; ++jrow) {
00275             i__3 = jrow + jcol * b_dim1;
00276             b[i__3].r = 0.f, b[i__3].i = 0.f;
00277 /* L10: */
00278         }
00279 /* L20: */
00280     }
00281 
00282 /*     Reduce A and B */
00283 
00284     i__1 = *ihi - 2;
00285     for (jcol = *ilo; jcol <= i__1; ++jcol) {
00286 
00287         i__2 = jcol + 2;
00288         for (jrow = *ihi; jrow >= i__2; --jrow) {
00289 
00290 /*           Step 1: rotate rows JROW-1, JROW to kill A(JROW,JCOL) */
00291 
00292             i__3 = jrow - 1 + jcol * a_dim1;
00293             ctemp.r = a[i__3].r, ctemp.i = a[i__3].i;
00294             clartg_(&ctemp, &a[jrow + jcol * a_dim1], &c__, &s, &a[jrow - 1 + 
00295                     jcol * a_dim1]);
00296             i__3 = jrow + jcol * a_dim1;
00297             a[i__3].r = 0.f, a[i__3].i = 0.f;
00298             i__3 = *n - jcol;
00299             crot_(&i__3, &a[jrow - 1 + (jcol + 1) * a_dim1], lda, &a[jrow + (
00300                     jcol + 1) * a_dim1], lda, &c__, &s);
00301             i__3 = *n + 2 - jrow;
00302             crot_(&i__3, &b[jrow - 1 + (jrow - 1) * b_dim1], ldb, &b[jrow + (
00303                     jrow - 1) * b_dim1], ldb, &c__, &s);
00304             if (ilq) {
00305                 r_cnjg(&q__1, &s);
00306                 crot_(n, &q[(jrow - 1) * q_dim1 + 1], &c__1, &q[jrow * q_dim1 
00307                         + 1], &c__1, &c__, &q__1);
00308             }
00309 
00310 /*           Step 2: rotate columns JROW, JROW-1 to kill B(JROW,JROW-1) */
00311 
00312             i__3 = jrow + jrow * b_dim1;
00313             ctemp.r = b[i__3].r, ctemp.i = b[i__3].i;
00314             clartg_(&ctemp, &b[jrow + (jrow - 1) * b_dim1], &c__, &s, &b[jrow 
00315                     + jrow * b_dim1]);
00316             i__3 = jrow + (jrow - 1) * b_dim1;
00317             b[i__3].r = 0.f, b[i__3].i = 0.f;
00318             crot_(ihi, &a[jrow * a_dim1 + 1], &c__1, &a[(jrow - 1) * a_dim1 + 
00319                     1], &c__1, &c__, &s);
00320             i__3 = jrow - 1;
00321             crot_(&i__3, &b[jrow * b_dim1 + 1], &c__1, &b[(jrow - 1) * b_dim1 
00322                     + 1], &c__1, &c__, &s);
00323             if (ilz) {
00324                 crot_(n, &z__[jrow * z_dim1 + 1], &c__1, &z__[(jrow - 1) * 
00325                         z_dim1 + 1], &c__1, &c__, &s);
00326             }
00327 /* L30: */
00328         }
00329 /* L40: */
00330     }
00331 
00332     return 0;
00333 
00334 /*     End of CGGHRD */
00335 
00336 } /* cgghrd_ */


swiftnav
Author(s):
autogenerated on Sat Jun 8 2019 18:55:27