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