1 /*****************************************************************************
2 *
3 * Elmer, A Finite Element Software for Multiphysical Problems
4 *
5 * Copyright 1st April 1995 - , CSC - IT Center for Science Ltd., Finland
6 *
7 * This library is free software; you can redistribute it and/or
8 * modify it under the terms of the GNU Lesser General Public
9 * License as published by the Free Software Foundation; either
10 * version 2.1 of the License, or (at your option) any later version.
11 *
12 * This library is distributed in the hope that it will be useful,
13 * but WITHOUT ANY WARRANTY; without even the implied warranty of
14 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
15 * Lesser General Public License for more details.
16 *
17 * You should have received a copy of the GNU Lesser General Public
18 * License along with this library (in file ../LGPL-2.1); if not, write
19 * to the Free Software Foundation, Inc., 51 Franklin Street,
20 * Fifth Floor, Boston, MA 02110-1301 USA
21 *
22 *****************************************************************************/
23
24 /*******************************************************************************
25 *
26 * The eigenvalue/eigenvector extraction routines.
27 *
28 *******************************************************************************
29 *
30 * Author: Juha Ruokolainen
31 *
32 * Address: CSC - IT Center for Science Ltd.
33 * Keilaranta 14, P.O. BOX 405
34 * 02101 Espoo, Finland
35 * Tel. +358 0 457 2723
36 * Telefax: +358 0 457 2302
37 * EMail: Juha.Ruokolainen@csc.fi
38 *
39 * Date: 30 May 1996
40 *
41 * Modified by:
42 *
43 * Date of modification:
44 *
45 ******************************************************************************/
46 /***********************************************************************
47 |
48 | EIG.C - Last Edited 8. 8. 1988
49 |
50 ***********************************************************************/
51
52 /*======================================================================
53 |Syntax of the manual pages:
54 |
55 |FUNCTION NAME(...) params ...
56 |
57 $ usage of the function and type of the parameters
58 ? explane the effects of the function
59 = return value and the type of value if not of type int
60 @ globals effected directly by this routine
61 ! current known bugs or limitations
62 & functions called by this function
63 ~ these functions may interest you as an alternative function or
64 | because they control this function somehow
65 ^=====================================================================*/
66
67
68 /*
69 * $Id: eig.c,v 1.1.1.1 2005/04/14 13:29:14 vierinen Exp $
70 *
71 * $Log: eig.c,v $
72 * Revision 1.1.1.1 2005/04/14 13:29:14 vierinen
73 * initial matc automake package
74 *
75 * Revision 1.2 1998/08/01 12:34:33 jpr
76 *
77 * Added Id, started Log.
78 *
79 *
80 */
81
82 #include "elmer/matc.h"
83
84 #define A(i,j) a[n * (i) + (j)]
85
86 #define MAXITER 1000
87 #define EPS 1e-16
88
mtr_hesse(var)89 VARIABLE *mtr_hesse(var)
90 VARIABLE *var;
91 {
92 VARIABLE *res;
93
94 double *a;
95
96 int n;
97
98 if (NCOL(var) != NROW(var))
99 {
100 error( "hesse: matrix must be square, current dimensions: [%d,%d]\n", NROW(var), NCOL(var) );
101 }
102
103 res = var_temp_copy(var);
104
105 a = MATR(res); n = NROW(res);
106 if (NROW(res) == 1) return res;
107
108 hesse(a, n, n);
109
110 return res;
111 }
112
mtr_eig(var)113 VARIABLE *mtr_eig(var)
114 VARIABLE *var;
115 {
116 VARIABLE *ptr, *res;
117
118 int iter, i, j, k, n;
119 double *a, b, s, t;
120
121 if (NCOL(var) != NROW(var))
122 {
123 error("eig: matrix must be square, current dimensions: [%d,%d]\n", NROW(var), NCOL(var));
124 }
125
126 ptr = var_temp_copy(var);
127
128 a = MATR(ptr); n = NROW(ptr);
129
130 if (NROW(ptr) == 1) return ptr;
131
132 hesse(a, n, n);
133
134 for(iter = 0; iter < MAXITER; iter++)
135 {
136
137 for (i = 0; i < n - 1; i++)
138 {
139 s = EPS*(abs(A(i,i))+abs(A(i+1,i+1)));
140 if (abs(A(i+1,i)) < s) A(i+1,i) = 0.0;
141 }
142
143 i = 0;
144 do
145 {
146 for(j = i; j < n - 1; j++)
147 if (A(j+1,j) != 0) break;
148
149 for(k = j; k < n - 1; k++)
150 if (A(k+1,k) == 0) break;
151
152 i = k;
153
154 } while(i < n - 1 && k - j + 1 < 3);
155
156 if (k - j + 1 < 3) break;
157
158 francis(&A(j,j), k - j + 1, n);
159 }
160
161 res = var_temp_new(TYPE_DOUBLE, n, 2);
162
163 for(i = 0, j = 0; i < n - 1; i++)
164 if (A(i+1,i) == 0)
165 M(res,j++,0) = A(i,i);
166 else
167 {
168 b = A(i,i) + A(i+1,i+1); s = b * b;
169 t = A(i,i) * A(i+1,i+1) - A(i,i+1)*A(i+1,i);
170 s = s - 4 * t;
171 if (s < 0)
172 {
173 M(res,j, 0) = b / 2;
174 M(res,j++, 1) = sqrt(-s) / 2;
175 M(res,j, 0) = b / 2;
176 M(res,j++,1) = -sqrt(-s) / 2;
177 }
178 else
179 {
180 M(res,j++, 0) = b / 2 + sqrt(s) / 2;
181 M(res,j++, 0) = b / 2 - sqrt(s) / 2;
182 }
183 i++;
184 }
185
186 if (A(n-1, n-2) == 0) M(res,j,0) = A(n-1, n-1);
187
188 var_delete_temp(ptr);
189
190 return res;
191 }
192
vbcalc(x,v,b,beg,end)193 void vbcalc(x,v,b,beg,end)
194 double x[],v[],*b;
195 int beg, end;
196 {
197 double alpha,m,mp1;
198 int i;
199
200 m = abs(x[beg]);
201 for(i = beg + 1; i <= end; i++)
202 m = max(m,abs(x[i]));
203
204 if (m < EPS)
205 {
206 /*
207 * for(i = beg; i <= end; i++) v[i] = 0;
208 */
209 memset(&v[beg], 0, (end-beg+1)<<3);
210 }
211 else
212 {
213 alpha = 0;
214 mp1 = 1 / m;
215 for(i = beg; i <= end; i++)
216 {
217 v[i] = x[i] * mp1;
218 alpha = alpha + v[i]*v[i];
219 }
220 alpha = sqrt(alpha);
221 *b = 1 / (alpha * (alpha + abs(v[beg])));
222 v[beg] = v[beg] + sign(v[beg]) * alpha;
223 }
224 return;
225 }
226
227 #define H(i,j) h[(i) * N + (j)]
228
hesse(h,DIM,N)229 void hesse(h, DIM, N)
230 int DIM, N;
231 double *h;
232 {
233 double *v, *x, b, s;
234 int i, j, k;
235
236 x = (double *)ALLOCMEM(DIM * sizeof(double));
237 v = (double *)ALLOCMEM(DIM * sizeof(double));
238
239 for (i = 0; i < DIM - 2; i++)
240 {
241
242 for(j = i + 1; j < DIM; j++) x[j] = H(j,i);
243
244 vbcalc(x, v, &b, i + 1, DIM - 1);
245
246 if (v[i+1] == 0) break;
247
248 for(j = i + 2; j < DIM; j++)
249 {
250 x[j] = v[j]/v[i+1];
251 v[j] = b * v[i+1] * v[j];
252 }
253 v[i+1] = b * v[i+1] * v[i+1];
254
255 for(j = 0; j < DIM; j++)
256 {
257 s = 0.0;
258 for(k = i + 1; k < DIM; k++)
259 s = s + H(j,k) * v[k];
260 H(j,i+1) = H(j,i+1) - s;
261 for(k = i + 2; k < DIM; k++)
262 H(j,k) = H(j,k) - s * x[k];
263 }
264
265 for(j = 0; j < DIM; j++)
266 {
267 s = H(i+1,j);
268 for(k = i + 2; k < DIM; k++)
269 s = s + H(k,j) * x[k];
270 for(k = i + 1; k < DIM; k++)
271 H(k,j) = H(k,j) - s * v[k];
272 }
273
274 for(j = i + 2; j < DIM; j++) H(j,i) = 0;
275 }
276 FREEMEM((char *)x); FREEMEM((char *)v);
277
278 return;
279 }
280
281
francis(h,DIM,N)282 void francis(h, DIM, N)
283 int DIM, N;
284 double *h;
285 {
286 double x[3], v[3], b, s, t, bv, v0i;
287 int i, i1, j, k, n, m, end;
288
289 n = DIM - 1; m = n - 1;
290
291 t = H(m,m) * H(n,n) - H(m,n) * H(n,m);
292 s = H(m,m) + H(n,n);
293
294 x[0] = H(0,0)*H(0,0)+H(0,1)*H(1,0)-s*H(0,0)+t;
295 x[1] = H(1,0) * (H(0,0) + H(1,1) - s);
296 x[2] = H(1,0) * H(2,1);
297
298 vbcalc(x, v, &b, 0, 2);
299
300 if (v[0] == 0) return;
301
302
303
304 /*
305 * for(i = 1; i < 3; i++)
306 * {
307 * x[i] = v[i]/v[0];
308 * v[i] = b * v[0] * v[i];
309 * }
310 */
311 bv = b * v[0];
312
313 x[1] = v[1] / v[0];
314 v[1] = bv * v[1];
315
316 x[2] = v[2] / v[0];
317 v[2] = bv * v[2];
318
319
320
321 x[0] = 1; v[0] = b * v[0] * v[0];
322
323 for(i = 0; i < DIM; i++)
324 {
325 /*
326 * s = 0.0;
327 * for(j = 0; j < 3; j++)
328 * s = s + H(i,j) * v[j];
329 */
330 i1 = i * N;
331 s = h[i1] * v[0] + h[i1+1]*v[1] + h[i1+2]*v[2];
332
333
334 H(i,0) = H(i,0) - s;
335 /*
336 * for(j = 1; j < 3; j++)
337 * H(i,j) = H(i,j) - s * x[j];
338 */
339 h[i1+1] = h[i1+1] - s * x[1];
340 h[i1+2] = h[i1+2] - s * x[2];
341 }
342
343 for(i = 0; i < DIM; i++)
344 {
345 /*
346 * s = H(0,i);
347 * for(j = 1; j < 3; j++)
348 * s = s + H(j,i) * x[j];
349 */
350 s = h[i] + h[N+i]*x[1] + h[(N<<1)+i]*x[2];
351 /*
352 * for(j = 0; j < 3; j++)
353 * H(j,i) = H(j,i) - s * v[j];
354 */
355 h[i] = h[i] - s * v[0];
356 h[N+i] = h[N+i] - s * v[1];
357 h[(N<<1)+i] = h[(N<<1)+i] - s * v[2];
358 }
359
360 for (i = 0; i < DIM - 2; i++)
361 {
362
363 end = min(2, DIM - i - 2);
364 for(j = 0; j <= end; j++) x[j] = H(i+j+1,i);
365
366 vbcalc(x, v, &b, 0, end);
367
368 if (v[0] == 0) break;
369
370 for(j = 1; j <= end; j++)
371 {
372 x[j] = v[j]/v[0];
373 v[j] = b * v[0] * v[j];
374 }
375 x[0] = 1; v[0] = b * v[0] * v[0];
376
377 for(j = 0; j < DIM; j++)
378 {
379 s = 0.0;
380 for(k = 0; k <= end; k++)
381 s = s + H(j,i+k+1) * v[k];
382 H(j,i+1) = H(j,i+1) - s;
383 for(k = 1; k <= end; k++)
384 H(j,i+k+1) = H(j,i+k+1) - s * x[k];
385 }
386
387 for(j = 0; j < DIM; j++)
388 {
389 s = H(i+1,j);
390 for(k = 1; k <= end; k++)
391 s = s + H(i+k+1,j) * x[k];
392 for(k = 0; k <= end; k++)
393 H(i+k+1,j) = H(i+k+1,j) - s * v[k];
394 }
395
396 for(j = i + 2; j < DIM; j++) H(j,i) = 0;
397 }
398 return;
399 }
400