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 * Solve symmetric positive definite eigenvalue problem.
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 * $Id: jacobi.c,v 1.1.1.1 2005/04/14 13:29:14 vierinen Exp $
49 *
50 * $Log: jacobi.c,v $
51 * Revision 1.1.1.1 2005/04/14 13:29:14 vierinen
52 * initial matc automake package
53 *
54 * Revision 1.2 1998/08/01 12:34:43 jpr
55 *
56 * Added Id, started Log.
57 *
58 *
59 */
60
61 #include "elmer/matc.h"
62
mtr_jacob(a)63 VARIABLE *mtr_jacob(a) VARIABLE *a;
64 {
65 VARIABLE *x, *ev;
66
67 double *b, *d, rtol;
68 int i, n;
69
70 if (NROW(a) != NCOL(a))
71 error("Jacob: Matrix must be square.\n");
72
73 b = MATR(NEXT(a));
74 n = NROW(a);
75
76 if (NROW(NEXT(a)) != NCOL(NEXT(a)) || n != NROW(NEXT(a)))
77 error("Jacob: Matrix dimensions incompatible.\n");
78
79 rtol = *MATR(NEXT(NEXT(a)));
80
81 x = var_new("eigv", TYPE_DOUBLE, NROW(a), NCOL(a));
82 d = (double *)ALLOCMEM(n * sizeof(double));
83 ev = var_temp_new(TYPE_DOUBLE, 1, n);
84
85 matc_jacobi(MATR(a), b, MATR(x), MATR(ev), d, n, rtol);
86 FREEMEM((char *)d);
87
88 return ev;
89 };
90
91 /************************************************************************
92 P R O G R A M
93 To solve the generalized eigenproblem using the
94 generalized Jacobi iteration
95
96 INPUT
97
98 a(n,n) = Stiffness matrix (assumed positive definite)
99 b(n,n) = Mass matrix (assumed positive definite)
100 x(n,n) = Matrix storing eigenvectors on solution exit
101 eigv(n) = Vector storing eigenvalues on solution exit
102 d(n) = Working vector
103 n = Order of the matrices a and b
104 rtol = Convergence tolerance (usually set to 10.**-12)
105 nsmax = Maximum number of sweeps allowed (usually set to 15)
106
107 OUTPUT
108
109 a(n,n) = Diagonalized stiffness matrix
110 b(n,n) = Diagonalized mass matrix
111 x(n,n) = Eigenvectors stored columnwise
112 eigv(n) = Eigenvalues
113
114 *************************************************************************/
115
matc_jacobi(a,b,x,eigv,d,n,rtol)116 int matc_jacobi(a,b,x,eigv,d,n,rtol)
117 double a[],b[],x[],eigv[],d[],rtol;
118 int n;
119 {
120 register int i,j,k,ii,jj;
121 int nsmax=50, /* Max number of sweeps */
122 nsweep, /* Current sweep number */
123 nr,
124 jp1,jm1,kp1,km1,
125 convergence;
126
127 double eps,
128 eptola,eptolb,
129 akk,ajj,ab,
130 check,sqch,
131 d1,d2,den,
132 ca,cg,
133 aj,bj,ak,bk,
134 xj,xk,
135 tol,dif,
136 epsa,epsb,
137 bb; /* Scale mass matrix */
138
139 /************************************************************************
140 Initialize eigenvalue and eigenvector matrices
141 *************************************************************************/
142 for( i=0 ; i<n ; i++)
143 {
144 ii = i*n+i; /* address of diagonal element */
145 if (a[ii]<=0.0 || b[ii]<=0.0) return 0;
146 eigv[i] = d[i] = a[ii]/b[ii];
147 x[ii] = 1.0; /* make an unit matrix */
148 }
149 if(n==1) return 1; /* Return if single degree of freedom system */
150
151 /************************************************************************
152 Initialize sweep counter and begin iteration
153 *************************************************************************/
154 nr=n - 1;
155 for ( nsweep = 0; nsweep < nsmax; nsweep++)
156 {
157 /************************************************************************
158 Check if present off-diagonal element is large enough to require zeroing
159 *************************************************************************/
160
161 eps = pow(0.01, 2.0 * (nsweep + 1));
162 for(j=0 ; j<nr ; j++)
163 {
164 jj=j+1;
165 for(k=jj ; k<n ; k++)
166 {
167 eptola=( a[j*n+k]*a[j*n+k] / ( a[j*n+j]*a[k*n+k] ) );
168 eptolb=( b[j*n+k]*b[j*n+k] / ( b[j*n+j]*b[k*n+k] ) );
169 if ( eptola >= eps || eptolb >=eps )
170 {
171
172 /************************************************************************
173 if zeroing is required, calculate the rotation matrix elements
174 *************************************************************************/
175
176 akk=a[k*n+k]*b[j*n+k] - b[k*n+k]*a[j*n+k];
177 ajj=a[j*n+j]*b[j*n+k] - b[j*n+j]*a[j*n+k];
178 ab =a[j*n+j]*b[k*n+k] - a[k*n+k]*b[j*n+j];
179 check=(ab*ab+4.0*akk*ajj)/4.0;
180
181 if (check<=0.0)
182 {
183 printf("***Error solution stop in *jacobi*\n");
184 printf(" check = %20.14e\n", check);
185 return 1;
186 }
187 sqch= sqrt(check);
188 d1 = ab/2.0+sqch;
189 d2 = ab/2.0-sqch;
190 den = d1;
191 if ( abs(d2) > abs(d1) ) den=d2;
192
193 if (den == 0.0)
194 {
195 ca=0.0;
196 cg= -a[j*n+k] / a[k*n+k];
197 }
198 else
199 {
200 ca= akk/den;
201 cg= -ajj/den;
202 }
203 /************************************************************************
204 Perform the generalized rotation to zero the present off-diagonal element
205 *************************************************************************/
206 if (n != 2)
207 {
208 jp1=j+1;
209 jm1=j-1;
210 kp1=k+1;
211 km1=k-1;
212 /**************************************/
213 if ( jm1 >= 0)
214 for (i=0 ; i<=jm1 ; i++)
215 {
216 aj = a[i*n+j];
217 bj = b[i*n+j];
218 ak = a[i*n+k];
219 bk = b[i*n+k];
220 a[i*n+j] = aj+cg*ak;
221 b[i*n+j] = bj+cg*bk;
222 a[i*n+k] = ak+ca*aj;
223 b[i*n+k] = bk+ca*bj;
224 };
225 /**************************************/
226 if ((kp1-n+1) <= 0)
227 for (i=kp1 ; i<n ; i++)
228 {
229 aj = a[j*n+i];
230 bj = b[j*n+i];
231 ak = a[k*n+i];
232 bk = b[k*n+i];
233 a[j*n+i] = aj+cg*ak;
234 b[j*n+i] = bj+cg*bk;
235 a[k*n+i] = ak+ca*aj;
236 b[k*n+i] = bk+ca*bj;
237 };
238 /**************************************/
239 if ((jp1-km1) <= 0.0)
240 for (i=jp1 ; i<=km1 ; i++)
241 {
242 aj = a[j*n+i];
243 bj = b[j*n+i];
244 ak = a[i*n+k];
245 bk = b[i*n+k];
246 a[j*n+i] = aj+cg*ak;
247 b[j*n+i] = bj+cg*bk;
248 a[i*n+k] = ak+ca*aj;
249 b[i*n+k] = bk+ca*bj;
250 };
251 };
252
253 ak = a[k*n+k];
254 bk = b[k*n+k];
255 a[k*n+k] = ak+2.0*ca*a[j*n+k]+ca*ca*a[j*n+j];
256 b[k*n+k] = bk+2.0*ca*b[j*n+k]+ca*ca*b[j*n+j];
257 a[j*n+j] = a[j*n+j]+2.0*cg*a[j*n+k]+cg*cg*ak;
258 b[j*n+j] = b[j*n+j]+2.0*cg*b[j*n+k]+cg*cg*bk;
259 a[j*n+k] = 0.0;
260 b[j*n+k] = 0.0;
261
262 /************************************************************************
263 Update the eigenvector matrix after each rotation
264 *************************************************************************/
265 for (i=0 ; i<n ; i++)
266 {
267 xj = x[i*n+j];
268 xk = x[i*n+k];
269 x[i*n+j] = xj+cg*xk;
270 x[i*n+k] = xk+ca*xj;
271 };
272 };
273 };
274 };
275 /************************************************************************
276 Update the eigenvalues after each sweep
277 *************************************************************************/
278 for (i=0 ; i<n ; i++)
279 {
280 ii=i*n+i;
281 if (a[ii]<=0.0 || b[ii]<=0.0)
282 {
283 error( "*** Error solution stop in *jacobi*\n Matrix not positive definite.");
284 };
285 eigv[i] = a[ii] / b[ii];
286 };
287
288 /************************************************************************
289 check for convergence
290 *************************************************************************/
291 convergence = 1;
292 for(i=0 ; i<n ; i++)
293 {
294 tol = rtol*d[i];
295 dif = abs(eigv[i]-d[i]);
296 if (dif > tol) convergence = 0;
297 if (! convergence) break;
298 };
299 /************************************************************************
300 check if all off-diag elements are satisfactorily small
301 *************************************************************************/
302 if (convergence)
303 {
304 eps=rtol*rtol;
305 for (j=0 ; j<nr ; j++)
306 {
307 jj=j+1;
308 for (k=jj ; k<n ; k++)
309 {
310 epsa = ( a[j*n+k]*a[j*n+k] / ( a[j*n+j]*a[k*n+k] ) );
311 epsb = ( b[j*n+k]*b[j*n+k] / ( b[j*n+j]*b[k*n+k] ) );
312 if ( epsa >= eps || epsb >=eps ) convergence = 0;
313 if (! convergence) break;
314 };
315 if (! convergence) break;
316 };
317 };
318 if (! convergence)
319 {
320 for (i=0 ; i<n ; i++)
321 d[i] = eigv[i];
322 };
323 if (convergence) break;
324 };
325 /************************************************************************
326 Fill out bottom triangle of resultant matrices and scale eigenvectors
327 *************************************************************************/
328 for (i=0 ; i<n ; i++)
329 for (j=i ; j<n ; j++)
330 {
331 b[j*n+i] = b[i*n+j];
332 a[j*n+i] = a[i*n+j];
333 };
334
335 for (j=0 ; j<n ; j++)
336 {
337 bb = sqrt( b[j*n+j] );
338 for (k=0 ; k<n ; k++)
339 x[k*n+j] /= bb;
340 };
341 PrintOut( "jacobi: nsweeps %d\n", nsweep );
342 return 1 ;
343 }
344