"Fossies" - the Fresh Open Source Software Archive 
Member "laspack/eigenval.c" (8 Aug 1995, 12568 Bytes) of package /linux/privat/old/laspack.tgz:
As a special service "Fossies" has tried to format the requested source page into HTML format using (guessed) C and C++ source code syntax highlighting (style:
standard) with prefixed line numbers and
code folding option.
Alternatively you can here
view or
download the uninterpreted source code file.
1 /****************************************************************************/
2 /* eigenval.c */
3 /****************************************************************************/
4 /* */
5 /* estimation of extremal EIGENVALues */
6 /* */
7 /* Copyright (C) 1992-1995 Tomas Skalicky. All rights reserved. */
8 /* */
9 /****************************************************************************/
10 /* */
11 /* ANY USE OF THIS CODE CONSTITUTES ACCEPTANCE OF THE TERMS */
12 /* OF THE COPYRIGHT NOTICE (SEE FILE COPYRGHT.H) */
13 /* */
14 /****************************************************************************/
15
16 #include <stddef.h>
17 #include <math.h>
18
19 #include "laspack/eigenval.h"
20 #include "laspack/elcmp.h"
21 #include "laspack/errhandl.h"
22 #include "laspack/operats.h"
23 #include "laspack/rtc.h"
24 #include "laspack/copyrght.h"
25
26 typedef struct {
27 double MinEigenval;
28 double MaxEigenval;
29 PrecondProcType PrecondProcUsed;
30 double OmegaPrecondUsed;
31 } EigenvalInfoType;
32
33 /* accuracy for the estimation of extremal eigenvalues */
34 static double EigenvalEps = 1e-4;
35
36 static void EstimEigenvals(QMatrix *A, PrecondProcType PrecondProc, double OmegaPrecond);
37 static void SearchEigenval(size_t n, double *Alpha, double *Beta, size_t k,
38 double BoundMin, double BoundMax, Boolean *Found, double *Lambda);
39 static size_t NoSmallerEigenvals(size_t n, double *Alpha, double *Beta, double Lambda);
40
41 #define max(x, y) ((x) > (y) ? (x) : (y))
42 #define min(x, y) ((x) < (y) ? (x) : (y))
43
44 void SetEigenvalAccuracy(double Eps)
45 /* set accuracy for the estimation of extremal eigenvalues */
46 {
47 EigenvalEps = Eps;
48 }
49
50 double GetMinEigenval(QMatrix *A, PrecondProcType PrecondProc, double OmegaPrecond)
51 /* returns estimate for minimum eigenvalue of the matrix A */
52 {
53 double MinEigenval;
54
55 EigenvalInfoType *EigenvalInfo;
56
57 Q_Lock(A);
58
59 if (LASResult() == LASOK) {
60 EigenvalInfo = (EigenvalInfoType *)*(Q_EigenvalInfo(A));
61 /* if eigenvalues not estimated yet, ... */
62 if (EigenvalInfo == NULL) {
63 EigenvalInfo = (EigenvalInfoType *)malloc(sizeof(EigenvalInfoType));
64 if (EigenvalInfo != NULL) {
65 *(Q_EigenvalInfo(A)) = (void *)EigenvalInfo;
66 EstimEigenvals(A, PrecondProc, OmegaPrecond);
67 } else {
68 LASError(LASMemAllocErr, "GetMinEigenval", Q_GetName(A), NULL, NULL);
69 }
70 }
71
72 /* if eigenvalues estimated with an other preconditioner, ... */
73 if (EigenvalInfo->PrecondProcUsed != PrecondProc
74 || EigenvalInfo->OmegaPrecondUsed != OmegaPrecond) {
75 EstimEigenvals(A, PrecondProc, OmegaPrecond);
76 }
77
78 if (LASResult() == LASOK)
79 MinEigenval = EigenvalInfo->MinEigenval;
80 else
81 MinEigenval = 1.0;
82 } else {
83 MinEigenval = 1.0;
84 }
85
86 return(MinEigenval);
87 }
88
89 double GetMaxEigenval(QMatrix *A, PrecondProcType PrecondProc, double OmegaPrecond)
90 /* returns estimate for maximum eigenvalue of the matrix A */
91 {
92 double MaxEigenval;
93
94 EigenvalInfoType *EigenvalInfo;
95
96 Q_Lock(A);
97
98 if (LASResult() == LASOK) {
99 EigenvalInfo = (EigenvalInfoType *)*(Q_EigenvalInfo(A));
100 /* if eigenvalues not estimated yet, ... */
101 if (EigenvalInfo == NULL) {
102 EigenvalInfo = (EigenvalInfoType *)malloc(sizeof(EigenvalInfoType));
103 if (EigenvalInfo != NULL) {
104 *(Q_EigenvalInfo(A)) = (void *)EigenvalInfo;
105 EstimEigenvals(A, PrecondProc, OmegaPrecond);
106 } else {
107 LASError(LASMemAllocErr, "GetMaxEigenval", Q_GetName(A), NULL, NULL);
108 }
109 }
110
111 /* if eigenvalues estimated with an other preconditioner, ... */
112 if (EigenvalInfo->PrecondProcUsed != PrecondProc
113 || EigenvalInfo->OmegaPrecondUsed != OmegaPrecond) {
114 EstimEigenvals(A, PrecondProc, OmegaPrecond);
115 }
116
117 if (LASResult() == LASOK)
118 MaxEigenval = EigenvalInfo->MaxEigenval;
119 else
120 MaxEigenval = 1.0;
121 } else {
122 MaxEigenval = 1.0;
123 }
124
125 return(MaxEigenval);
126 }
127
128 static void EstimEigenvals(QMatrix *A, PrecondProcType PrecondProc, double OmegaPrecond)
129 /* estimates extremal eigenvalues of the matrix A by means of the Lanczos method */
130 {
131 /*
132 * for details to the Lanczos algorithm see
133 *
134 * G. H. Golub, Ch. F. van Loan:
135 * Matrix Computations;
136 * North Oxford Academic, Oxford, 1986
137 *
138 * (for modification for preconditioned matrices compare with sec. 10.3)
139 *
140 */
141
142 double LambdaMin = 0.0, LambdaMax = 0.0;
143 double LambdaMinOld, LambdaMaxOld;
144 double GershBoundMin = 0.0, GershBoundMax = 0.0;
145 double *Alpha, *Beta;
146 size_t Dim, j;
147 Boolean Found;
148 Vector q, qOld, h, p;
149
150 Q_Lock(A);
151
152 Dim = Q_GetDim(A);
153 V_Constr(&q, "q", Dim, Normal, True);
154 V_Constr(&qOld, "qOld", Dim, Normal, True);
155 V_Constr(&h, "h", Dim, Normal, True);
156 if (PrecondProc != NULL)
157 V_Constr(&p, "p", Dim, Normal, True);
158
159 if (LASResult() == LASOK) {
160 Alpha = (double *)malloc((Dim + 1) * sizeof(double));
161 Beta = (double *)malloc((Dim + 1) * sizeof(double));
162 if (Alpha != NULL && Beta != NULL) {
163 j = 0;
164
165 V_SetAllCmp(&qOld, 0.0);
166 V_SetRndCmp(&q);
167 if (Q_KerDefined(A))
168 OrthoRightKer_VQ(&q, A);
169 if (Q_GetSymmetry(A) && PrecondProc != NULL) {
170 (*PrecondProc)(A, &p, &q, OmegaPrecond);
171 MulAsgn_VS(&q, 1.0 / sqrt(Mul_VV(&q, &p)));
172 } else {
173 MulAsgn_VS(&q, 1.0 / l2Norm_V(&q));
174 }
175
176 Beta[0] = 1.0;
177 do {
178 j++;
179 if (Q_GetSymmetry(A) && PrecondProc != NULL) {
180 /* p = M^(-1) q */
181 (*PrecondProc)(A, &p, &q, OmegaPrecond);
182 /* h = A p */
183 Asgn_VV(&h, Mul_QV(A, &p));
184 if (Q_KerDefined(A))
185 OrthoRightKer_VQ(&h, A);
186 /* Alpha = p . h */
187 Alpha[j] = Mul_VV(&p, &h);
188 /* r = h - Alpha q - Beta qOld */
189 SubAsgn_VV(&h, Add_VV(Mul_SV(Alpha[j], &q), Mul_SV(Beta[j-1], &qOld)));
190 /* z = M^(-1) r */
191 (*PrecondProc)(A, &p, &h, OmegaPrecond);
192 /* Beta = sqrt(r . z) */
193 Beta[j] = sqrt(Mul_VV(&h, &p));
194 Asgn_VV(&qOld, &q);
195 /* q = r / Beta */
196 Asgn_VV(&q, Mul_SV(1.0 / Beta[j], &h));
197 } else {
198 /* h = A p */
199 if (Q_GetSymmetry(A)) {
200 Asgn_VV(&h, Mul_QV(A, &q));
201 } else {
202 if (PrecondProc != NULL) {
203 (*PrecondProc)(A, &h, Mul_QV(A, &q), OmegaPrecond);
204 (*PrecondProc)(Transp_Q(A), &h, &h, OmegaPrecond);
205 Asgn_VV(&h, Mul_QV(Transp_Q(A), &h));
206 } else {
207 Asgn_VV(&h, Mul_QV(Transp_Q(A), Mul_QV(A, &q)));
208 }
209 }
210 if (Q_KerDefined(A))
211 OrthoRightKer_VQ(&h, A);
212 /* Alpha = q . h */
213 Alpha[j] = Mul_VV(&q, &h);
214 /* r = h - Alpha q - Beta qOld */
215 SubAsgn_VV(&h, Add_VV(Mul_SV(Alpha[j], &q), Mul_SV(Beta[j-1], &qOld)));
216 /* Beta = || r || */
217 Beta[j] = l2Norm_V(&h);
218 Asgn_VV(&qOld, &q);
219 /* q = r / Beta */
220 Asgn_VV(&q, Mul_SV(1.0 / Beta[j], &h));
221 }
222
223 LambdaMaxOld = LambdaMax;
224 LambdaMinOld = LambdaMin;
225
226 /* determination of extremal eigenvalues of the tridiagonal matrix
227 (Beta[i-1] Alpha[i] Beta[i]) (where 1 <= i <= j)
228 by means of the method of bisection; bounds for eigenvalues
229 are determined after Gershgorin circle theorem */
230 if (j == 1) {
231 GershBoundMin = Alpha[1] - fabs(Beta[1]);
232 GershBoundMax = Alpha[1] + fabs(Beta[1]);
233
234 LambdaMin = Alpha[1];
235 LambdaMax = Alpha[1];
236 } else {
237 GershBoundMin = min(Alpha[j] - fabs(Beta[j]) - fabs(Beta[j - 1]),
238 GershBoundMin);
239 GershBoundMax = max(Alpha[j] + fabs(Beta[j]) + fabs(Beta[j - 1]),
240 GershBoundMax);
241
242 SearchEigenval(j, Alpha, Beta, 1, GershBoundMin, LambdaMin,
243 &Found, &LambdaMin);
244 if (!Found)
245 SearchEigenval(j, Alpha, Beta, 1, GershBoundMin, GershBoundMax,
246 &Found, &LambdaMin);
247
248 SearchEigenval(j, Alpha, Beta, j, LambdaMax, GershBoundMax,
249 &Found, &LambdaMax);
250 if (!Found)
251 SearchEigenval(j, Alpha, Beta, j, GershBoundMin, GershBoundMax,
252 &Found, &LambdaMax);
253 }
254 } while (!IsZero(Beta[j]) && j < Dim
255 && (fabs(LambdaMin - LambdaMinOld) > EigenvalEps * LambdaMin
256 || fabs(LambdaMax - LambdaMaxOld) > EigenvalEps * LambdaMax)
257 && LASResult() == LASOK);
258
259 if (Q_GetSymmetry(A)) {
260 LambdaMin = (1.0 - j * EigenvalEps) * LambdaMin;
261 } else {
262 LambdaMin = (1.0 - sqrt(j) * EigenvalEps) * sqrt(LambdaMin);
263 }
264 if (Alpha != NULL)
265 free(Alpha);
266 if (Beta != NULL)
267 free(Beta);
268 } else {
269 LASError(LASMemAllocErr, "EstimEigenvals", Q_GetName(A), NULL, NULL);
270 }
271
272 }
273
274 V_Destr(&q);
275 V_Destr(&qOld);
276 V_Destr(&h);
277 if (PrecondProc != NULL)
278 V_Destr(&p);
279
280 if (LASResult() == LASOK) {
281 ((EigenvalInfoType *)*(Q_EigenvalInfo(A)))->MinEigenval = LambdaMin;
282 ((EigenvalInfoType *)*(Q_EigenvalInfo(A)))->MaxEigenval = LambdaMax;
283 ((EigenvalInfoType *)*(Q_EigenvalInfo(A)))->PrecondProcUsed = PrecondProc;
284 ((EigenvalInfoType *)*(Q_EigenvalInfo(A)))->OmegaPrecondUsed = OmegaPrecond;
285 } else {
286 ((EigenvalInfoType *)*(Q_EigenvalInfo(A)))->MinEigenval = 1.0;
287 ((EigenvalInfoType *)*(Q_EigenvalInfo(A)))->MaxEigenval = 1.0;
288 ((EigenvalInfoType *)*(Q_EigenvalInfo(A)))->PrecondProcUsed = NULL;
289 ((EigenvalInfoType *)*(Q_EigenvalInfo(A)))->OmegaPrecondUsed = 1.0;
290 }
291
292 Q_Unlock(A);
293 }
294
295 static void SearchEigenval(size_t n, double *Alpha, double *Beta, size_t k,
296 double BoundMin, double BoundMax, Boolean *Found, double *Lambda)
297 /* search the k-th eigenvalue of the tridiagonal matrix
298 (Beta[i-1] Alpha[i] Beta[i]) (where 1 <= i <= n)
299 by means of the method of bisection */
300 {
301 /*
302 * for details to the method of bisection see
303 *
304 * G. H. Golub, Ch. F. van Loan:
305 * Matrix Computations;
306 * North Oxford Academic, Oxford, 1986
307 *
308 */
309
310 if (NoSmallerEigenvals(n, Alpha, Beta, BoundMin) < k
311 && NoSmallerEigenvals(n, Alpha, Beta, BoundMax) >= k) {
312 while (fabs(BoundMax - BoundMin) > 0.01 * EigenvalEps
313 * (fabs(BoundMin) + fabs(BoundMax))) {
314 *Lambda = 0.5 * (BoundMin + BoundMax);
315 if (NoSmallerEigenvals(n, Alpha, Beta, *Lambda) >= k)
316 BoundMax = *Lambda;
317 else
318 BoundMin = *Lambda;
319 }
320 *Lambda = BoundMax;
321
322 *Found = True;
323 } else {
324 *Found = False;
325 }
326 }
327
328 static size_t NoSmallerEigenvals(size_t n, double *Alpha, double *Beta, double Lambda)
329 /* returns number of eigenvalues of the tridiagonal matrix
330 (Beta[i-1] Alpha[i] Beta[i]) (where 1 <= i <= n)
331 which are less then Lambda */
332 {
333 size_t No;
334
335 double p, pNew, pOld, Sign;
336 size_t i;
337
338 No = 0;
339
340 pOld = 1.0;
341 p = (Alpha[1] - Lambda) / fabs(Beta[1]);
342 /* check for change of sign */
343 if (IsZero(p) || p * pOld < 0)
344 No++;
345
346 for (i = 2; i <= n; i++) {
347 Sign = Beta[i-1] / fabs(Beta[i-1]);
348 pNew = ((Alpha[i] - Lambda) * p - Beta[i-1] * Sign * pOld) / fabs(Beta[i]);
349 pOld = p;
350 p = pNew;
351
352 /* check for change of sign */
353 if (p * pOld < 0 || (IsZero(p) && !IsZero(pOld)))
354 No++;
355 }
356
357 return(No);
358 }