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