Doxygen Source Code Documentation
eis_htridi.c File Reference
#include "f2c.h"
Go to the source code of this file.
Functions | |
int | htridi_ (integer *nm, integer *n, doublereal *ar, doublereal *ai, doublereal *d__, doublereal *e, doublereal *e2, doublereal *tau) |
Function Documentation
|
Definition at line 8 of file eis_htridi.c. References abs, l, pythag_(), and scale. Referenced by ch_().
00011 { 00012 /* System generated locals */ 00013 integer ar_dim1, ar_offset, ai_dim1, ai_offset, i__1, i__2, i__3; 00014 doublereal d__1, d__2; 00015 00016 /* Builtin functions */ 00017 double sqrt(doublereal); 00018 00019 /* Local variables */ 00020 static doublereal f, g, h__; 00021 static integer i__, j, k, l; 00022 static doublereal scale, fi, gi, hh; 00023 static integer ii; 00024 static doublereal si; 00025 extern doublereal pythag_(doublereal *, doublereal *); 00026 static integer jp1; 00027 00028 00029 00030 /* THIS SUBROUTINE IS A TRANSLATION OF A COMPLEX ANALOGUE OF */ 00031 /* THE ALGOL PROCEDURE TRED1, NUM. MATH. 11, 181-195(1968) */ 00032 /* BY MARTIN, REINSCH, AND WILKINSON. */ 00033 /* HANDBOOK FOR AUTO. COMP., VOL.II-LINEAR ALGEBRA, 212-226(1971). */ 00034 00035 /* THIS SUBROUTINE REDUCES A COMPLEX HERMITIAN MATRIX */ 00036 /* TO A REAL SYMMETRIC TRIDIAGONAL MATRIX USING */ 00037 /* UNITARY SIMILARITY TRANSFORMATIONS. */ 00038 00039 /* ON INPUT */ 00040 00041 /* NM MUST BE SET TO THE ROW DIMENSION OF TWO-DIMENSIONAL */ 00042 /* ARRAY PARAMETERS AS DECLARED IN THE CALLING PROGRAM */ 00043 /* DIMENSION STATEMENT. */ 00044 00045 /* N IS THE ORDER OF THE MATRIX. */ 00046 00047 /* AR AND AI CONTAIN THE REAL AND IMAGINARY PARTS, */ 00048 /* RESPECTIVELY, OF THE COMPLEX HERMITIAN INPUT MATRIX. */ 00049 /* ONLY THE LOWER TRIANGLE OF THE MATRIX NEED BE SUPPLIED. */ 00050 00051 /* ON OUTPUT */ 00052 00053 /* AR AND AI CONTAIN INFORMATION ABOUT THE UNITARY TRANS- */ 00054 /* FORMATIONS USED IN THE REDUCTION IN THEIR FULL LOWER */ 00055 /* TRIANGLES. THEIR STRICT UPPER TRIANGLES AND THE */ 00056 /* DIAGONAL OF AR ARE UNALTERED. */ 00057 00058 /* D CONTAINS THE DIAGONAL ELEMENTS OF THE THE TRIDIAGONAL MATRIX. 00059 */ 00060 00061 /* E CONTAINS THE SUBDIAGONAL ELEMENTS OF THE TRIDIAGONAL */ 00062 /* MATRIX IN ITS LAST N-1 POSITIONS. E(1) IS SET TO ZERO. */ 00063 00064 /* E2 CONTAINS THE SQUARES OF THE CORRESPONDING ELEMENTS OF E. */ 00065 /* E2 MAY COINCIDE WITH E IF THE SQUARES ARE NOT NEEDED. */ 00066 00067 /* TAU CONTAINS FURTHER INFORMATION ABOUT THE TRANSFORMATIONS. */ 00068 00069 /* CALLS PYTHAG FOR DSQRT(A*A + B*B) . */ 00070 00071 /* QUESTIONS AND COMMENTS SHOULD BE DIRECTED TO BURTON S. GARBOW, */ 00072 /* MATHEMATICS AND COMPUTER SCIENCE DIV, ARGONNE NATIONAL LABORATORY 00073 */ 00074 00075 /* THIS VERSION DATED AUGUST 1983. */ 00076 00077 /* ------------------------------------------------------------------ 00078 */ 00079 00080 /* Parameter adjustments */ 00081 tau -= 3; 00082 --e2; 00083 --e; 00084 --d__; 00085 ai_dim1 = *nm; 00086 ai_offset = ai_dim1 + 1; 00087 ai -= ai_offset; 00088 ar_dim1 = *nm; 00089 ar_offset = ar_dim1 + 1; 00090 ar -= ar_offset; 00091 00092 /* Function Body */ 00093 tau[(*n << 1) + 1] = 1.; 00094 tau[(*n << 1) + 2] = 0.; 00095 00096 i__1 = *n; 00097 for (i__ = 1; i__ <= i__1; ++i__) { 00098 /* L100: */ 00099 d__[i__] = ar[i__ + i__ * ar_dim1]; 00100 } 00101 /* .......... FOR I=N STEP -1 UNTIL 1 DO -- .......... */ 00102 i__1 = *n; 00103 for (ii = 1; ii <= i__1; ++ii) { 00104 i__ = *n + 1 - ii; 00105 l = i__ - 1; 00106 h__ = 0.; 00107 scale = 0.; 00108 if (l < 1) { 00109 goto L130; 00110 } 00111 /* .......... SCALE ROW (ALGOL TOL THEN NOT NEEDED) .......... */ 00112 i__2 = l; 00113 for (k = 1; k <= i__2; ++k) { 00114 /* L120: */ 00115 scale = scale + (d__1 = ar[i__ + k * ar_dim1], abs(d__1)) + (d__2 00116 = ai[i__ + k * ai_dim1], abs(d__2)); 00117 } 00118 00119 if (scale != 0.) { 00120 goto L140; 00121 } 00122 tau[(l << 1) + 1] = 1.; 00123 tau[(l << 1) + 2] = 0.; 00124 L130: 00125 e[i__] = 0.; 00126 e2[i__] = 0.; 00127 goto L290; 00128 00129 L140: 00130 i__2 = l; 00131 for (k = 1; k <= i__2; ++k) { 00132 ar[i__ + k * ar_dim1] /= scale; 00133 ai[i__ + k * ai_dim1] /= scale; 00134 h__ = h__ + ar[i__ + k * ar_dim1] * ar[i__ + k * ar_dim1] + ai[ 00135 i__ + k * ai_dim1] * ai[i__ + k * ai_dim1]; 00136 /* L150: */ 00137 } 00138 00139 e2[i__] = scale * scale * h__; 00140 g = sqrt(h__); 00141 e[i__] = scale * g; 00142 f = pythag_(&ar[i__ + l * ar_dim1], &ai[i__ + l * ai_dim1]); 00143 /* .......... FORM NEXT DIAGONAL ELEMENT OF MATRIX T .......... */ 00144 if (f == 0.) { 00145 goto L160; 00146 } 00147 tau[(l << 1) + 1] = (ai[i__ + l * ai_dim1] * tau[(i__ << 1) + 2] - ar[ 00148 i__ + l * ar_dim1] * tau[(i__ << 1) + 1]) / f; 00149 si = (ar[i__ + l * ar_dim1] * tau[(i__ << 1) + 2] + ai[i__ + l * 00150 ai_dim1] * tau[(i__ << 1) + 1]) / f; 00151 h__ += f * g; 00152 g = g / f + 1.; 00153 ar[i__ + l * ar_dim1] = g * ar[i__ + l * ar_dim1]; 00154 ai[i__ + l * ai_dim1] = g * ai[i__ + l * ai_dim1]; 00155 if (l == 1) { 00156 goto L270; 00157 } 00158 goto L170; 00159 L160: 00160 tau[(l << 1) + 1] = -tau[(i__ << 1) + 1]; 00161 si = tau[(i__ << 1) + 2]; 00162 ar[i__ + l * ar_dim1] = g; 00163 L170: 00164 f = 0.; 00165 00166 i__2 = l; 00167 for (j = 1; j <= i__2; ++j) { 00168 g = 0.; 00169 gi = 0.; 00170 /* .......... FORM ELEMENT OF A*U .......... */ 00171 i__3 = j; 00172 for (k = 1; k <= i__3; ++k) { 00173 g = g + ar[j + k * ar_dim1] * ar[i__ + k * ar_dim1] + ai[j + 00174 k * ai_dim1] * ai[i__ + k * ai_dim1]; 00175 gi = gi - ar[j + k * ar_dim1] * ai[i__ + k * ai_dim1] + ai[j 00176 + k * ai_dim1] * ar[i__ + k * ar_dim1]; 00177 /* L180: */ 00178 } 00179 00180 jp1 = j + 1; 00181 if (l < jp1) { 00182 goto L220; 00183 } 00184 00185 i__3 = l; 00186 for (k = jp1; k <= i__3; ++k) { 00187 g = g + ar[k + j * ar_dim1] * ar[i__ + k * ar_dim1] - ai[k + 00188 j * ai_dim1] * ai[i__ + k * ai_dim1]; 00189 gi = gi - ar[k + j * ar_dim1] * ai[i__ + k * ai_dim1] - ai[k 00190 + j * ai_dim1] * ar[i__ + k * ar_dim1]; 00191 /* L200: */ 00192 } 00193 /* .......... FORM ELEMENT OF P .......... */ 00194 L220: 00195 e[j] = g / h__; 00196 tau[(j << 1) + 2] = gi / h__; 00197 f = f + e[j] * ar[i__ + j * ar_dim1] - tau[(j << 1) + 2] * ai[i__ 00198 + j * ai_dim1]; 00199 /* L240: */ 00200 } 00201 00202 hh = f / (h__ + h__); 00203 /* .......... FORM REDUCED A .......... */ 00204 i__2 = l; 00205 for (j = 1; j <= i__2; ++j) { 00206 f = ar[i__ + j * ar_dim1]; 00207 g = e[j] - hh * f; 00208 e[j] = g; 00209 fi = -ai[i__ + j * ai_dim1]; 00210 gi = tau[(j << 1) + 2] - hh * fi; 00211 tau[(j << 1) + 2] = -gi; 00212 00213 i__3 = j; 00214 for (k = 1; k <= i__3; ++k) { 00215 ar[j + k * ar_dim1] = ar[j + k * ar_dim1] - f * e[k] - g * ar[ 00216 i__ + k * ar_dim1] + fi * tau[(k << 1) + 2] + gi * ai[ 00217 i__ + k * ai_dim1]; 00218 ai[j + k * ai_dim1] = ai[j + k * ai_dim1] - f * tau[(k << 1) 00219 + 2] - g * ai[i__ + k * ai_dim1] - fi * e[k] - gi * 00220 ar[i__ + k * ar_dim1]; 00221 /* L260: */ 00222 } 00223 } 00224 00225 L270: 00226 i__3 = l; 00227 for (k = 1; k <= i__3; ++k) { 00228 ar[i__ + k * ar_dim1] = scale * ar[i__ + k * ar_dim1]; 00229 ai[i__ + k * ai_dim1] = scale * ai[i__ + k * ai_dim1]; 00230 /* L280: */ 00231 } 00232 00233 tau[(l << 1) + 2] = -si; 00234 L290: 00235 hh = d__[i__]; 00236 d__[i__] = ar[i__ + i__ * ar_dim1]; 00237 ar[i__ + i__ * ar_dim1] = hh; 00238 ai[i__ + i__ * ai_dim1] = scale * sqrt(h__); 00239 /* L300: */ 00240 } 00241 00242 return 0; 00243 } /* htridi_ */ |