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