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