00001 /* htribk.f -- translated by f2c (version 19961017).
00002    You must link the resulting object file with the libraries:
00003         -lf2c -lm   (in that order)
00004 */
00005 
00006 #include "f2c.h"
00007 
00008 /* Subroutine */ int htribk_(integer *nm, integer *n, doublereal *ar, 
00009         doublereal *ai, doublereal *tau, integer *m, doublereal *zr, 
00010         doublereal *zi)
00011 {
00012     /* System generated locals */
00013     integer ar_dim1, ar_offset, ai_dim1, ai_offset, zr_dim1, zr_offset, 
00014             zi_dim1, zi_offset, i__1, i__2, i__3;
00015 
00016     /* Local variables */
00017     static doublereal h__;
00018     static integer i__, j, k, l;
00019     static doublereal s, si;
00020 
00021 
00022 
00023 /*     THIS SUBROUTINE IS A TRANSLATION OF A COMPLEX ANALOGUE OF */
00024 /*     THE ALGOL PROCEDURE TRBAK1, NUM. MATH. 11, 181-195(1968) */
00025 /*     BY MARTIN, REINSCH, AND WILKINSON. */
00026 /*     HANDBOOK FOR AUTO. COMP., VOL.II-LINEAR ALGEBRA, 212-226(1971). */
00027 
00028 /*     THIS SUBROUTINE FORMS THE EIGENVECTORS OF A COMPLEX HERMITIAN */
00029 /*     MATRIX BY BACK TRANSFORMING THOSE OF THE CORRESPONDING */
00030 /*     REAL SYMMETRIC TRIDIAGONAL MATRIX DETERMINED BY  HTRIDI. */
00031 
00032 /*     ON INPUT */
00033 
00034 /*        NM MUST BE SET TO THE ROW DIMENSION OF TWO-DIMENSIONAL */
00035 /*          ARRAY PARAMETERS AS DECLARED IN THE CALLING PROGRAM */
00036 /*          DIMENSION STATEMENT. */
00037 
00038 /*        N IS THE ORDER OF THE MATRIX. */
00039 
00040 /*        AR AND AI CONTAIN INFORMATION ABOUT THE UNITARY TRANS- */
00041 /*          FORMATIONS USED IN THE REDUCTION BY  HTRIDI  IN THEIR */
00042 /*          FULL LOWER TRIANGLES EXCEPT FOR THE DIAGONAL OF AR. */
00043 
00044 /*        TAU CONTAINS FURTHER INFORMATION ABOUT THE TRANSFORMATIONS. */
00045 
00046 /*        M IS THE NUMBER OF EIGENVECTORS TO BE BACK TRANSFORMED. */
00047 
00048 /*        ZR CONTAINS THE EIGENVECTORS TO BE BACK TRANSFORMED */
00049 /*          IN ITS FIRST M COLUMNS. */
00050 
00051 /*     ON OUTPUT */
00052 
00053 /*        ZR AND ZI CONTAIN THE REAL AND IMAGINARY PARTS, */
00054 /*          RESPECTIVELY, OF THE TRANSFORMED EIGENVECTORS */
00055 /*          IN THEIR FIRST M COLUMNS. */
00056 
00057 /*     NOTE THAT THE LAST COMPONENT OF EACH RETURNED VECTOR */
00058 /*     IS REAL AND THAT VECTOR EUCLIDEAN NORMS ARE PRESERVED. */
00059 
00060 /*     QUESTIONS AND COMMENTS SHOULD BE DIRECTED TO BURTON S. GARBOW, */
00061 /*     MATHEMATICS AND COMPUTER SCIENCE DIV, ARGONNE NATIONAL LABORATORY 
00062 */
00063 
00064 /*     THIS VERSION DATED AUGUST 1983. */
00065 
00066 /*     ------------------------------------------------------------------ 
00067 */
00068 
00069     /* Parameter adjustments */
00070     tau -= 3;
00071     ai_dim1 = *nm;
00072     ai_offset = ai_dim1 + 1;
00073     ai -= ai_offset;
00074     ar_dim1 = *nm;
00075     ar_offset = ar_dim1 + 1;
00076     ar -= ar_offset;
00077     zi_dim1 = *nm;
00078     zi_offset = zi_dim1 + 1;
00079     zi -= zi_offset;
00080     zr_dim1 = *nm;
00081     zr_offset = zr_dim1 + 1;
00082     zr -= zr_offset;
00083 
00084     /* Function Body */
00085     if (*m == 0) {
00086         goto L200;
00087     }
00088 /*     .......... TRANSFORM THE EIGENVECTORS OF THE REAL SYMMETRIC */
00089 /*                TRIDIAGONAL MATRIX TO THOSE OF THE HERMITIAN */
00090 /*                TRIDIAGONAL MATRIX. .......... */
00091     i__1 = *n;
00092     for (k = 1; k <= i__1; ++k) {
00093 
00094         i__2 = *m;
00095         for (j = 1; j <= i__2; ++j) {
00096             zi[k + j * zi_dim1] = -zr[k + j * zr_dim1] * tau[(k << 1) + 2];
00097             zr[k + j * zr_dim1] *= tau[(k << 1) + 1];
00098 /* L50: */
00099         }
00100     }
00101 
00102     if (*n == 1) {
00103         goto L200;
00104     }
00105 /*     .......... RECOVER AND APPLY THE HOUSEHOLDER MATRICES .......... */
00106     i__2 = *n;
00107     for (i__ = 2; i__ <= i__2; ++i__) {
00108         l = i__ - 1;
00109         h__ = ai[i__ + i__ * ai_dim1];
00110         if (h__ == 0.) {
00111             goto L140;
00112         }
00113 
00114         i__1 = *m;
00115         for (j = 1; j <= i__1; ++j) {
00116             s = 0.;
00117             si = 0.;
00118 
00119             i__3 = l;
00120             for (k = 1; k <= i__3; ++k) {
00121                 s = s + ar[i__ + k * ar_dim1] * zr[k + j * zr_dim1] - ai[i__ 
00122                         + k * ai_dim1] * zi[k + j * zi_dim1];
00123                 si = si + ar[i__ + k * ar_dim1] * zi[k + j * zi_dim1] + ai[
00124                         i__ + k * ai_dim1] * zr[k + j * zr_dim1];
00125 /* L110: */
00126             }
00127 /*     .......... DOUBLE DIVISIONS AVOID POSSIBLE UNDERFLOW ......
00128 .... */
00129             s = s / h__ / h__;
00130             si = si / h__ / h__;
00131 
00132             i__3 = l;
00133             for (k = 1; k <= i__3; ++k) {
00134                 zr[k + j * zr_dim1] = zr[k + j * zr_dim1] - s * ar[i__ + k * 
00135                         ar_dim1] - si * ai[i__ + k * ai_dim1];
00136                 zi[k + j * zi_dim1] = zi[k + j * zi_dim1] - si * ar[i__ + k * 
00137                         ar_dim1] + s * ai[i__ + k * ai_dim1];
00138 /* L120: */
00139             }
00140 
00141 /* L130: */
00142         }
00143 
00144 L140:
00145         ;
00146     }
00147 
00148 L200:
00149     return 0;
00150 } /* htribk_ */
00151