00001 /* ortran.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 ortran_(integer *nm, integer *n, integer *low, integer *
00009         igh, doublereal *a, doublereal *ort, doublereal *z__)
00010 {
00011     /* System generated locals */
00012     integer a_dim1, a_offset, z_dim1, z_offset, i__1, i__2, i__3;
00013 
00014     /* Local variables */
00015     static doublereal g;
00016     static integer i__, j, kl, mm, mp, mp1;
00017 
00018 
00019 
00020 /*     THIS SUBROUTINE IS A TRANSLATION OF THE ALGOL PROCEDURE ORTRANS, */
00021 /*     NUM. MATH. 16, 181-204(1970) BY PETERS AND WILKINSON. */
00022 /*     HANDBOOK FOR AUTO. COMP., VOL.II-LINEAR ALGEBRA, 372-395(1971). */
00023 
00024 /*     THIS SUBROUTINE ACCUMULATES THE ORTHOGONAL SIMILARITY */
00025 /*     TRANSFORMATIONS USED IN THE REDUCTION OF A REAL GENERAL */
00026 /*     MATRIX TO UPPER HESSENBERG FORM BY  ORTHES. */
00027 
00028 /*     ON INPUT */
00029 
00030 /*        NM MUST BE SET TO THE ROW DIMENSION OF TWO-DIMENSIONAL */
00031 /*          ARRAY PARAMETERS AS DECLARED IN THE CALLING PROGRAM */
00032 /*          DIMENSION STATEMENT. */
00033 
00034 /*        N IS THE ORDER OF THE MATRIX. */
00035 
00036 /*        LOW AND IGH ARE INTEGERS DETERMINED BY THE BALANCING */
00037 /*          SUBROUTINE  BALANC.  IF  BALANC  HAS NOT BEEN USED, */
00038 /*          SET LOW=1, IGH=N. */
00039 
00040 /*        A CONTAINS INFORMATION ABOUT THE ORTHOGONAL TRANS- */
00041 /*          FORMATIONS USED IN THE REDUCTION BY  ORTHES */
00042 /*          IN ITS STRICT LOWER TRIANGLE. */
00043 
00044 /*        ORT CONTAINS FURTHER INFORMATION ABOUT THE TRANS- */
00045 /*          FORMATIONS USED IN THE REDUCTION BY  ORTHES. */
00046 /*          ONLY ELEMENTS LOW THROUGH IGH ARE USED. */
00047 
00048 /*     ON OUTPUT */
00049 
00050 /*        Z CONTAINS THE TRANSFORMATION MATRIX PRODUCED IN THE */
00051 /*          REDUCTION BY  ORTHES. */
00052 
00053 /*        ORT HAS BEEN ALTERED. */
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 /*     .......... INITIALIZE Z TO IDENTITY MATRIX .......... */
00065     /* Parameter adjustments */
00066     z_dim1 = *nm;
00067     z_offset = z_dim1 + 1;
00068     z__ -= z_offset;
00069     --ort;
00070     a_dim1 = *nm;
00071     a_offset = a_dim1 + 1;
00072     a -= a_offset;
00073 
00074     /* Function Body */
00075     i__1 = *n;
00076     for (j = 1; j <= i__1; ++j) {
00077 
00078         i__2 = *n;
00079         for (i__ = 1; i__ <= i__2; ++i__) {
00080 /* L60: */
00081             z__[i__ + j * z_dim1] = 0.;
00082         }
00083 
00084         z__[j + j * z_dim1] = 1.;
00085 /* L80: */
00086     }
00087 
00088     kl = *igh - *low - 1;
00089     if (kl < 1) {
00090         goto L200;
00091     }
00092 /*     .......... FOR MP=IGH-1 STEP -1 UNTIL LOW+1 DO -- .......... */
00093     i__1 = kl;
00094     for (mm = 1; mm <= i__1; ++mm) {
00095         mp = *igh - mm;
00096         if (a[mp + (mp - 1) * a_dim1] == 0.) {
00097             goto L140;
00098         }
00099         mp1 = mp + 1;
00100 
00101         i__2 = *igh;
00102         for (i__ = mp1; i__ <= i__2; ++i__) {
00103 /* L100: */
00104             ort[i__] = a[i__ + (mp - 1) * a_dim1];
00105         }
00106 
00107         i__2 = *igh;
00108         for (j = mp; j <= i__2; ++j) {
00109             g = 0.;
00110 
00111             i__3 = *igh;
00112             for (i__ = mp; i__ <= i__3; ++i__) {
00113 /* L110: */
00114                 g += ort[i__] * z__[i__ + j * z_dim1];
00115             }
00116 /*     .......... DIVISOR BELOW IS NEGATIVE OF H FORMED IN ORTHES.
00117  */
00118 /*                DOUBLE DIVISION AVOIDS POSSIBLE UNDERFLOW ......
00119 .... */
00120             g = g / ort[mp] / a[mp + (mp - 1) * a_dim1];
00121 
00122             i__3 = *igh;
00123             for (i__ = mp; i__ <= i__3; ++i__) {
00124 /* L120: */
00125                 z__[i__ + j * z_dim1] += g * ort[i__];
00126             }
00127 
00128 /* L130: */
00129         }
00130 
00131 L140:
00132         ;
00133     }
00134 
00135 L200:
00136     return 0;
00137 } /* ortran_ */
00138