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