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