Doxygen Source Code Documentation
eis_reduc.c File Reference
#include "f2c.h"
Go to the source code of this file.
Functions | |
int | reduc_ (integer *nm, integer *n, doublereal *a, doublereal *b, doublereal *dl, integer *ierr) |
Function Documentation
|
Definition at line 8 of file eis_reduc.c. Referenced by rsg_().
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_ */ |