00001 /* reduc.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 reduc_(integer *nm, integer *n, doublereal *a, 
00009         doublereal *b, doublereal *dl, integer *ierr)
00010 {
00011     /* System generated locals */
00012     integer a_dim1, a_offset, b_dim1, b_offset, i__1, i__2, i__3;
00013 
00014     /* Builtin functions */
00015     double sqrt(doublereal);
00016 
00017     /* Local variables */
00018     static integer i__, j, k;
00019     static doublereal x, y;
00020     static integer i1, j1, nn;
00021 
00022 
00023 
00024 /*     THIS SUBROUTINE IS A TRANSLATION OF THE ALGOL PROCEDURE REDUC1, */
00025 /*     NUM. MATH. 11, 99-110(1968) BY MARTIN AND WILKINSON. */
00026 /*     HANDBOOK FOR AUTO. COMP., VOL.II-LINEAR ALGEBRA, 303-314(1971). */
00027 
00028 /*     THIS SUBROUTINE REDUCES THE GENERALIZED SYMMETRIC EIGENPROBLEM */
00029 /*     AX=(LAMBDA)BX, WHERE B IS POSITIVE DEFINITE, TO THE STANDARD */
00030 /*     SYMMETRIC EIGENPROBLEM USING THE CHOLESKY FACTORIZATION OF B. */
00031 
00032 /*     ON INPUT */
00033 
00034 /*        NM MUST BE SET TO THE ROW DIMENSION OF TWO-DIMENSIONAL */
00035 /*          ARRAY PARAMETERS AS DECLARED IN THE CALLING PROGRAM */
00036 /*          DIMENSION STATEMENT. */
00037 
00038 /*        N IS THE ORDER OF THE MATRICES A AND B.  IF THE CHOLESKY */
00039 /*          FACTOR L OF B IS ALREADY AVAILABLE, N SHOULD BE PREFIXED */
00040 /*          WITH A MINUS SIGN. */
00041 
00042 /*        A AND B CONTAIN THE REAL SYMMETRIC INPUT MATRICES.  ONLY THE */
00043 /*          FULL UPPER TRIANGLES OF THE MATRICES NEED BE SUPPLIED.  IF */
00044 /*          N IS NEGATIVE, THE STRICT LOWER TRIANGLE OF B CONTAINS, */
00045 /*          INSTEAD, THE STRICT LOWER TRIANGLE OF ITS CHOLESKY FACTOR L. 
00046 */
00047 
00048 /*        DL CONTAINS, IF N IS NEGATIVE, THE DIAGONAL ELEMENTS OF L. */
00049 
00050 /*     ON OUTPUT */
00051 
00052 /*        A CONTAINS IN ITS FULL LOWER TRIANGLE THE FULL LOWER TRIANGLE */
00053 /*          OF THE SYMMETRIC MATRIX DERIVED FROM THE REDUCTION TO THE */
00054 /*          STANDARD FORM.  THE STRICT UPPER TRIANGLE OF A IS UNALTERED. 
00055 */
00056 
00057 /*        B CONTAINS IN ITS STRICT LOWER TRIANGLE THE STRICT LOWER */
00058 /*          TRIANGLE OF ITS CHOLESKY FACTOR L.  THE FULL UPPER */
00059 /*          TRIANGLE OF B IS UNALTERED. */
00060 
00061 /*        DL CONTAINS THE DIAGONAL ELEMENTS OF L. */
00062 
00063 /*        IERR IS SET TO */
00064 /*          ZERO       FOR NORMAL RETURN, */
00065 /*          7*N+1      IF B IS NOT POSITIVE DEFINITE. */
00066 
00067 /*     QUESTIONS AND COMMENTS SHOULD BE DIRECTED TO BURTON S. GARBOW, */
00068 /*     MATHEMATICS AND COMPUTER SCIENCE DIV, ARGONNE NATIONAL LABORATORY 
00069 */
00070 
00071 /*     THIS VERSION DATED AUGUST 1983. */
00072 
00073 /*     ------------------------------------------------------------------ 
00074 */
00075 
00076     /* Parameter adjustments */
00077     --dl;
00078     b_dim1 = *nm;
00079     b_offset = b_dim1 + 1;
00080     b -= b_offset;
00081     a_dim1 = *nm;
00082     a_offset = a_dim1 + 1;
00083     a -= a_offset;
00084 
00085     /* Function Body */
00086     *ierr = 0;
00087     nn = abs(*n);
00088     if (*n < 0) {
00089         goto L100;
00090     }
00091 /*     .......... FORM L IN THE ARRAYS B AND DL .......... */
00092     i__1 = *n;
00093     for (i__ = 1; i__ <= i__1; ++i__) {
00094         i1 = i__ - 1;
00095 
00096         i__2 = *n;
00097         for (j = i__; j <= i__2; ++j) {
00098             x = b[i__ + j * b_dim1];
00099             if (i__ == 1) {
00100                 goto L40;
00101             }
00102 
00103             i__3 = i1;
00104             for (k = 1; k <= i__3; ++k) {
00105 /* L20: */
00106                 x -= b[i__ + k * b_dim1] * b[j + k * b_dim1];
00107             }
00108 
00109 L40:
00110             if (j != i__) {
00111                 goto L60;
00112             }
00113             if (x <= 0.) {
00114                 goto L1000;
00115             }
00116             y = sqrt(x);
00117             dl[i__] = y;
00118             goto L80;
00119 L60:
00120             b[j + i__ * b_dim1] = x / y;
00121 L80:
00122             ;
00123         }
00124     }
00125 /*     .......... FORM THE TRANSPOSE OF THE UPPER TRIANGLE OF INV(L)*A */
00126 /*                IN THE LOWER TRIANGLE OF THE ARRAY A .......... */
00127 L100:
00128     i__2 = nn;
00129     for (i__ = 1; i__ <= i__2; ++i__) {
00130         i1 = i__ - 1;
00131         y = dl[i__];
00132 
00133         i__1 = nn;
00134         for (j = i__; j <= i__1; ++j) {
00135             x = a[i__ + j * a_dim1];
00136             if (i__ == 1) {
00137                 goto L180;
00138             }
00139 
00140             i__3 = i1;
00141             for (k = 1; k <= i__3; ++k) {
00142 /* L160: */
00143                 x -= b[i__ + k * b_dim1] * a[j + k * a_dim1];
00144             }
00145 
00146 L180:
00147             a[j + i__ * a_dim1] = x / y;
00148 /* L200: */
00149         }
00150     }
00151 /*     .......... PRE-MULTIPLY BY INV(L) AND OVERWRITE .......... */
00152     i__1 = nn;
00153     for (j = 1; j <= i__1; ++j) {
00154         j1 = j - 1;
00155 
00156         i__2 = nn;
00157         for (i__ = j; i__ <= i__2; ++i__) {
00158             x = a[i__ + j * a_dim1];
00159             if (i__ == j) {
00160                 goto L240;
00161             }
00162             i1 = i__ - 1;
00163 
00164             i__3 = i1;
00165             for (k = j; k <= i__3; ++k) {
00166 /* L220: */
00167                 x -= a[k + j * a_dim1] * b[i__ + k * b_dim1];
00168             }
00169 
00170 L240:
00171             if (j == 1) {
00172                 goto L280;
00173             }
00174 
00175             i__3 = j1;
00176             for (k = 1; k <= i__3; ++k) {
00177 /* L260: */
00178                 x -= a[j + k * a_dim1] * b[i__ + k * b_dim1];
00179             }
00180 
00181 L280:
00182             a[i__ + j * a_dim1] = x / dl[i__];
00183 /* L300: */
00184         }
00185     }
00186 
00187     goto L1001;
00188 /*     .......... SET ERROR -- B IS NOT POSITIVE DEFINITE .......... */
00189 L1000:
00190     *ierr = *n * 7 + 1;
00191 L1001:
00192     return 0;
00193 } /* reduc_ */
00194