00001 /* cgehd2.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 integer c__1 = 1; 00019 00020 /* Subroutine */ int cgehd2_(integer *n, integer *ilo, integer *ihi, complex * 00021 a, integer *lda, complex *tau, complex *work, integer *info) 00022 { 00023 /* System generated locals */ 00024 integer a_dim1, a_offset, i__1, i__2, i__3; 00025 complex q__1; 00026 00027 /* Builtin functions */ 00028 void r_cnjg(complex *, complex *); 00029 00030 /* Local variables */ 00031 integer i__; 00032 complex alpha; 00033 extern /* Subroutine */ int clarf_(char *, integer *, integer *, complex * 00034 , integer *, complex *, complex *, integer *, complex *), 00035 clarfg_(integer *, complex *, complex *, integer *, complex *), 00036 xerbla_(char *, integer *); 00037 00038 00039 /* -- LAPACK routine (version 3.2) -- */ 00040 /* Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. */ 00041 /* November 2006 */ 00042 00043 /* .. Scalar Arguments .. */ 00044 /* .. */ 00045 /* .. Array Arguments .. */ 00046 /* .. */ 00047 00048 /* Purpose */ 00049 /* ======= */ 00050 00051 /* CGEHD2 reduces a complex general matrix A to upper Hessenberg form H */ 00052 /* by a unitary similarity transformation: Q' * A * Q = H . */ 00053 00054 /* Arguments */ 00055 /* ========= */ 00056 00057 /* N (input) INTEGER */ 00058 /* The order of the matrix A. N >= 0. */ 00059 00060 /* ILO (input) INTEGER */ 00061 /* IHI (input) INTEGER */ 00062 /* It is assumed that A is already upper triangular in rows */ 00063 /* and columns 1:ILO-1 and IHI+1:N. ILO and IHI are normally */ 00064 /* set by a previous call to CGEBAL; otherwise they should be */ 00065 /* set to 1 and N respectively. See Further Details. */ 00066 /* 1 <= ILO <= IHI <= max(1,N). */ 00067 00068 /* A (input/output) COMPLEX array, dimension (LDA,N) */ 00069 /* On entry, the n by n general matrix to be reduced. */ 00070 /* On exit, the upper triangle and the first subdiagonal of A */ 00071 /* are overwritten with the upper Hessenberg matrix H, and the */ 00072 /* elements below the first subdiagonal, with the array TAU, */ 00073 /* represent the unitary matrix Q as a product of elementary */ 00074 /* reflectors. See Further Details. */ 00075 00076 /* LDA (input) INTEGER */ 00077 /* The leading dimension of the array A. LDA >= max(1,N). */ 00078 00079 /* TAU (output) COMPLEX array, dimension (N-1) */ 00080 /* The scalar factors of the elementary reflectors (see Further */ 00081 /* Details). */ 00082 00083 /* WORK (workspace) COMPLEX array, dimension (N) */ 00084 00085 /* INFO (output) INTEGER */ 00086 /* = 0: successful exit */ 00087 /* < 0: if INFO = -i, the i-th argument had an illegal value. */ 00088 00089 /* Further Details */ 00090 /* =============== */ 00091 00092 /* The matrix Q is represented as a product of (ihi-ilo) elementary */ 00093 /* reflectors */ 00094 00095 /* Q = H(ilo) H(ilo+1) . . . H(ihi-1). */ 00096 00097 /* Each H(i) has the form */ 00098 00099 /* H(i) = I - tau * v * v' */ 00100 00101 /* where tau is a complex scalar, and v is a complex vector with */ 00102 /* v(1:i) = 0, v(i+1) = 1 and v(ihi+1:n) = 0; v(i+2:ihi) is stored on */ 00103 /* exit in A(i+2:ihi,i), and tau in TAU(i). */ 00104 00105 /* The contents of A are illustrated by the following example, with */ 00106 /* n = 7, ilo = 2 and ihi = 6: */ 00107 00108 /* on entry, on exit, */ 00109 00110 /* ( a a a a a a a ) ( a a h h h h a ) */ 00111 /* ( a a a a a a ) ( a h h h h a ) */ 00112 /* ( a a a a a a ) ( h h h h h h ) */ 00113 /* ( a a a a a a ) ( v2 h h h h h ) */ 00114 /* ( a a a a a a ) ( v2 v3 h h h h ) */ 00115 /* ( a a a a a a ) ( v2 v3 v4 h h h ) */ 00116 /* ( a ) ( a ) */ 00117 00118 /* where a denotes an element of the original matrix A, h denotes a */ 00119 /* modified element of the upper Hessenberg matrix H, and vi denotes an */ 00120 /* element of the vector defining H(i). */ 00121 00122 /* ===================================================================== */ 00123 00124 /* .. Parameters .. */ 00125 /* .. */ 00126 /* .. Local Scalars .. */ 00127 /* .. */ 00128 /* .. External Subroutines .. */ 00129 /* .. */ 00130 /* .. Intrinsic Functions .. */ 00131 /* .. */ 00132 /* .. Executable Statements .. */ 00133 00134 /* Test the input parameters */ 00135 00136 /* Parameter adjustments */ 00137 a_dim1 = *lda; 00138 a_offset = 1 + a_dim1; 00139 a -= a_offset; 00140 --tau; 00141 --work; 00142 00143 /* Function Body */ 00144 *info = 0; 00145 if (*n < 0) { 00146 *info = -1; 00147 } else if (*ilo < 1 || *ilo > max(1,*n)) { 00148 *info = -2; 00149 } else if (*ihi < min(*ilo,*n) || *ihi > *n) { 00150 *info = -3; 00151 } else if (*lda < max(1,*n)) { 00152 *info = -5; 00153 } 00154 if (*info != 0) { 00155 i__1 = -(*info); 00156 xerbla_("CGEHD2", &i__1); 00157 return 0; 00158 } 00159 00160 i__1 = *ihi - 1; 00161 for (i__ = *ilo; i__ <= i__1; ++i__) { 00162 00163 /* Compute elementary reflector H(i) to annihilate A(i+2:ihi,i) */ 00164 00165 i__2 = i__ + 1 + i__ * a_dim1; 00166 alpha.r = a[i__2].r, alpha.i = a[i__2].i; 00167 i__2 = *ihi - i__; 00168 /* Computing MIN */ 00169 i__3 = i__ + 2; 00170 clarfg_(&i__2, &alpha, &a[min(i__3, *n)+ i__ * a_dim1], &c__1, &tau[ 00171 i__]); 00172 i__2 = i__ + 1 + i__ * a_dim1; 00173 a[i__2].r = 1.f, a[i__2].i = 0.f; 00174 00175 /* Apply H(i) to A(1:ihi,i+1:ihi) from the right */ 00176 00177 i__2 = *ihi - i__; 00178 clarf_("Right", ihi, &i__2, &a[i__ + 1 + i__ * a_dim1], &c__1, &tau[ 00179 i__], &a[(i__ + 1) * a_dim1 + 1], lda, &work[1]); 00180 00181 /* Apply H(i)' to A(i+1:ihi,i+1:n) from the left */ 00182 00183 i__2 = *ihi - i__; 00184 i__3 = *n - i__; 00185 r_cnjg(&q__1, &tau[i__]); 00186 clarf_("Left", &i__2, &i__3, &a[i__ + 1 + i__ * a_dim1], &c__1, &q__1, 00187 &a[i__ + 1 + (i__ + 1) * a_dim1], lda, &work[1]); 00188 00189 i__2 = i__ + 1 + i__ * a_dim1; 00190 a[i__2].r = alpha.r, a[i__2].i = alpha.i; 00191 /* L10: */ 00192 } 00193 00194 return 0; 00195 00196 /* End of CGEHD2 */ 00197 00198 } /* cgehd2_ */