00001 /* comhes.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 comhes_(integer *nm, integer *n, integer *low, integer *
00009         igh, doublereal *ar, doublereal *ai, integer *int__)
00010 {
00011     /* System generated locals */
00012     integer ar_dim1, ar_offset, ai_dim1, ai_offset, i__1, i__2, i__3;
00013     doublereal d__1, d__2;
00014 
00015     /* Local variables */
00016     extern /* Subroutine */ int cdiv_(doublereal *, doublereal *, doublereal *
00017             , doublereal *, doublereal *, doublereal *);
00018     static integer i__, j, m, la;
00019     static doublereal xi, yi, xr, yr;
00020     static integer mm1, kp1, mp1;
00021 
00022 
00023 
00024 /*     THIS SUBROUTINE IS A TRANSLATION OF THE ALGOL PROCEDURE COMHES, */
00025 /*     NUM. MATH. 12, 349-368(1968) BY MARTIN AND WILKINSON. */
00026 /*     HANDBOOK FOR AUTO. COMP., VOL.II-LINEAR ALGEBRA, 339-358(1971). */
00027 
00028 /*     GIVEN A COMPLEX GENERAL MATRIX, THIS SUBROUTINE */
00029 /*     REDUCES A SUBMATRIX SITUATED IN ROWS AND COLUMNS */
00030 /*     LOW THROUGH IGH TO UPPER HESSENBERG FORM BY */
00031 /*     STABILIZED ELEMENTARY SIMILARITY TRANSFORMATIONS. */
00032 
00033 /*     ON INPUT */
00034 
00035 /*        NM MUST BE SET TO THE ROW DIMENSION OF TWO-DIMENSIONAL */
00036 /*          ARRAY PARAMETERS AS DECLARED IN THE CALLING PROGRAM */
00037 /*          DIMENSION STATEMENT. */
00038 
00039 /*        N IS THE ORDER OF THE MATRIX. */
00040 
00041 /*        LOW AND IGH ARE INTEGERS DETERMINED BY THE BALANCING */
00042 /*          SUBROUTINE  CBAL.  IF  CBAL  HAS NOT BEEN USED, */
00043 /*          SET LOW=1, IGH=N. */
00044 
00045 /*        AR AND AI CONTAIN THE REAL AND IMAGINARY PARTS, */
00046 /*          RESPECTIVELY, OF THE COMPLEX INPUT MATRIX. */
00047 
00048 /*     ON OUTPUT */
00049 
00050 /*        AR AND AI CONTAIN THE REAL AND IMAGINARY PARTS, */
00051 /*          RESPECTIVELY, OF THE HESSENBERG MATRIX.  THE */
00052 /*          MULTIPLIERS WHICH WERE USED IN THE REDUCTION */
00053 /*          ARE STORED IN THE REMAINING TRIANGLES UNDER THE */
00054 /*          HESSENBERG MATRIX. */
00055 
00056 /*        INT CONTAINS INFORMATION ON THE ROWS AND COLUMNS */
00057 /*          INTERCHANGED IN THE REDUCTION. */
00058 /*          ONLY ELEMENTS LOW THROUGH IGH ARE USED. */
00059 
00060 /*     CALLS CDIV FOR COMPLEX DIVISION. */
00061 
00062 /*     QUESTIONS AND COMMENTS SHOULD BE DIRECTED TO BURTON S. GARBOW, */
00063 /*     MATHEMATICS AND COMPUTER SCIENCE DIV, ARGONNE NATIONAL LABORATORY 
00064 */
00065 
00066 /*     THIS VERSION DATED AUGUST 1983. */
00067 
00068 /*     ------------------------------------------------------------------ 
00069 */
00070 
00071     /* Parameter adjustments */
00072     ai_dim1 = *nm;
00073     ai_offset = ai_dim1 + 1;
00074     ai -= ai_offset;
00075     ar_dim1 = *nm;
00076     ar_offset = ar_dim1 + 1;
00077     ar -= ar_offset;
00078     --int__;
00079 
00080     /* Function Body */
00081     la = *igh - 1;
00082     kp1 = *low + 1;
00083     if (la < kp1) {
00084         goto L200;
00085     }
00086 
00087     i__1 = la;
00088     for (m = kp1; m <= i__1; ++m) {
00089         mm1 = m - 1;
00090         xr = 0.;
00091         xi = 0.;
00092         i__ = m;
00093 
00094         i__2 = *igh;
00095         for (j = m; j <= i__2; ++j) {
00096             if ((d__1 = ar[j + mm1 * ar_dim1], abs(d__1)) + (d__2 = ai[j + 
00097                     mm1 * ai_dim1], abs(d__2)) <= abs(xr) + abs(xi)) {
00098                 goto L100;
00099             }
00100             xr = ar[j + mm1 * ar_dim1];
00101             xi = ai[j + mm1 * ai_dim1];
00102             i__ = j;
00103 L100:
00104             ;
00105         }
00106 
00107         int__[m] = i__;
00108         if (i__ == m) {
00109             goto L130;
00110         }
00111 /*     .......... INTERCHANGE ROWS AND COLUMNS OF AR AND AI ..........
00112  */
00113         i__2 = *n;
00114         for (j = mm1; j <= i__2; ++j) {
00115             yr = ar[i__ + j * ar_dim1];
00116             ar[i__ + j * ar_dim1] = ar[m + j * ar_dim1];
00117             ar[m + j * ar_dim1] = yr;
00118             yi = ai[i__ + j * ai_dim1];
00119             ai[i__ + j * ai_dim1] = ai[m + j * ai_dim1];
00120             ai[m + j * ai_dim1] = yi;
00121 /* L110: */
00122         }
00123 
00124         i__2 = *igh;
00125         for (j = 1; j <= i__2; ++j) {
00126             yr = ar[j + i__ * ar_dim1];
00127             ar[j + i__ * ar_dim1] = ar[j + m * ar_dim1];
00128             ar[j + m * ar_dim1] = yr;
00129             yi = ai[j + i__ * ai_dim1];
00130             ai[j + i__ * ai_dim1] = ai[j + m * ai_dim1];
00131             ai[j + m * ai_dim1] = yi;
00132 /* L120: */
00133         }
00134 /*     .......... END INTERCHANGE .......... */
00135 L130:
00136         if (xr == 0. && xi == 0.) {
00137             goto L180;
00138         }
00139         mp1 = m + 1;
00140 
00141         i__2 = *igh;
00142         for (i__ = mp1; i__ <= i__2; ++i__) {
00143             yr = ar[i__ + mm1 * ar_dim1];
00144             yi = ai[i__ + mm1 * ai_dim1];
00145             if (yr == 0. && yi == 0.) {
00146                 goto L160;
00147             }
00148             cdiv_(&yr, &yi, &xr, &xi, &yr, &yi);
00149             ar[i__ + mm1 * ar_dim1] = yr;
00150             ai[i__ + mm1 * ai_dim1] = yi;
00151 
00152             i__3 = *n;
00153             for (j = m; j <= i__3; ++j) {
00154                 ar[i__ + j * ar_dim1] = ar[i__ + j * ar_dim1] - yr * ar[m + j 
00155                         * ar_dim1] + yi * ai[m + j * ai_dim1];
00156                 ai[i__ + j * ai_dim1] = ai[i__ + j * ai_dim1] - yr * ai[m + j 
00157                         * ai_dim1] - yi * ar[m + j * ar_dim1];
00158 /* L140: */
00159             }
00160 
00161             i__3 = *igh;
00162             for (j = 1; j <= i__3; ++j) {
00163                 ar[j + m * ar_dim1] = ar[j + m * ar_dim1] + yr * ar[j + i__ * 
00164                         ar_dim1] - yi * ai[j + i__ * ai_dim1];
00165                 ai[j + m * ai_dim1] = ai[j + m * ai_dim1] + yr * ai[j + i__ * 
00166                         ai_dim1] + yi * ar[j + i__ * ar_dim1];
00167 /* L150: */
00168             }
00169 
00170 L160:
00171             ;
00172         }
00173 
00174 L180:
00175         ;
00176     }
00177 
00178 L200:
00179     return 0;
00180 } /* comhes_ */
00181