00001 /* ortbak.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 ortbak_(integer *nm, integer *low, integer *igh, 
00009         doublereal *a, doublereal *ort, integer *m, 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, la, mm, mp, kp1, mp1;
00017 
00018 
00019 
00020 /*     THIS SUBROUTINE IS A TRANSLATION OF THE ALGOL PROCEDURE ORTBAK, */
00021 /*     NUM. MATH. 12, 349-368(1968) BY MARTIN AND WILKINSON. */
00022 /*     HANDBOOK FOR AUTO. COMP., VOL.II-LINEAR ALGEBRA, 339-358(1971). */
00023 
00024 /*     THIS SUBROUTINE FORMS THE EIGENVECTORS OF A REAL GENERAL */
00025 /*     MATRIX BY BACK TRANSFORMING THOSE OF THE CORRESPONDING */
00026 /*     UPPER HESSENBERG MATRIX DETERMINED 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 /*        LOW AND IGH ARE INTEGERS DETERMINED BY THE BALANCING */
00035 /*          SUBROUTINE  BALANC.  IF  BALANC  HAS NOT BEEN USED, */
00036 /*          SET LOW=1 AND IGH EQUAL TO THE ORDER OF THE MATRIX. */
00037 
00038 /*        A CONTAINS INFORMATION ABOUT THE ORTHOGONAL TRANS- */
00039 /*          FORMATIONS USED IN THE REDUCTION BY  ORTHES */
00040 /*          IN ITS STRICT LOWER TRIANGLE. */
00041 
00042 /*        ORT CONTAINS FURTHER INFORMATION ABOUT THE TRANS- */
00043 /*          FORMATIONS USED IN THE REDUCTION BY  ORTHES. */
00044 /*          ONLY ELEMENTS LOW THROUGH IGH ARE USED. */
00045 
00046 /*        M IS THE NUMBER OF COLUMNS OF Z TO BE BACK TRANSFORMED. */
00047 
00048 /*        Z CONTAINS THE REAL AND IMAGINARY PARTS OF THE EIGEN- */
00049 /*          VECTORS TO BE BACK TRANSFORMED IN ITS FIRST M COLUMNS. */
00050 
00051 /*     ON OUTPUT */
00052 
00053 /*        Z CONTAINS THE REAL AND IMAGINARY PARTS OF THE */
00054 /*          TRANSFORMED EIGENVECTORS IN ITS FIRST M COLUMNS. */
00055 
00056 /*        ORT HAS BEEN ALTERED. */
00057 
00058 /*     NOTE THAT ORTBAK PRESERVES VECTOR EUCLIDEAN NORMS. */
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     --ort;
00071     a_dim1 = *nm;
00072     a_offset = a_dim1 + 1;
00073     a -= a_offset;
00074     z_dim1 = *nm;
00075     z_offset = z_dim1 + 1;
00076     z__ -= z_offset;
00077 
00078     /* Function Body */
00079     if (*m == 0) {
00080         goto L200;
00081     }
00082     la = *igh - 1;
00083     kp1 = *low + 1;
00084     if (la < kp1) {
00085         goto L200;
00086     }
00087 /*     .......... FOR MP=IGH-1 STEP -1 UNTIL LOW+1 DO -- .......... */
00088     i__1 = la;
00089     for (mm = kp1; mm <= i__1; ++mm) {
00090         mp = *low + *igh - mm;
00091         if (a[mp + (mp - 1) * a_dim1] == 0.) {
00092             goto L140;
00093         }
00094         mp1 = mp + 1;
00095 
00096         i__2 = *igh;
00097         for (i__ = mp1; i__ <= i__2; ++i__) {
00098 /* L100: */
00099             ort[i__] = a[i__ + (mp - 1) * a_dim1];
00100         }
00101 
00102         i__2 = *m;
00103         for (j = 1; j <= i__2; ++j) {
00104             g = 0.;
00105 
00106             i__3 = *igh;
00107             for (i__ = mp; i__ <= i__3; ++i__) {
00108 /* L110: */
00109                 g += ort[i__] * z__[i__ + j * z_dim1];
00110             }
00111 /*     .......... DIVISOR BELOW IS NEGATIVE OF H FORMED IN ORTHES.
00112  */
00113 /*                DOUBLE DIVISION AVOIDS POSSIBLE UNDERFLOW ......
00114 .... */
00115             g = g / ort[mp] / a[mp + (mp - 1) * a_dim1];
00116 
00117             i__3 = *igh;
00118             for (i__ = mp; i__ <= i__3; ++i__) {
00119 /* L120: */
00120                 z__[i__ + j * z_dim1] += g * ort[i__];
00121             }
00122 
00123 /* L130: */
00124         }
00125 
00126 L140:
00127         ;
00128     }
00129 
00130 L200:
00131     return 0;
00132 } /* ortbak_ */
00133