OOFEM  2.4 OOFEM.org - Object Oriented Finite Element Solver
gjacobi.C
Go to the documentation of this file.
1 /*
2  *
3  * ##### ##### ###### ###### ### ###
4  * ## ## ## ## ## ## ## ### ##
5  * ## ## ## ## #### #### ## # ##
6  * ## ## ## ## ## ## ## ##
7  * ## ## ## ## ## ## ## ##
8  * ##### ##### ## ###### ## ##
9  *
10  *
11  * OOFEM : Object Oriented Finite Element Code
12  *
13  * Copyright (C) 1993 - 2013 Borek Patzak
14  *
15  *
16  *
17  * Czech Technical University, Faculty of Civil Engineering,
18  * Department of Structural Mechanics, 166 29 Prague, Czech Republic
19  *
20  * This library is free software; you can redistribute it and/or
21  * modify it under the terms of the GNU Lesser General Public
22  * License as published by the Free Software Foundation; either
23  * version 2.1 of the License, or (at your option) any later version.
24  *
25  * This program is distributed in the hope that it will be useful,
26  * but WITHOUT ANY WARRANTY; without even the implied warranty of
27  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
28  * Lesser General Public License for more details.
29  *
30  * You should have received a copy of the GNU Lesser General Public
31  * License along with this library; if not, write to the Free Software
32  * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
33  */
34
35
36 #include "gjacobi.h"
37 #include "mathfem.h"
38 #include "floatmatrix.h"
39 #include "nmstatus.h"
40
41 namespace oofem {
43  NumericalMethod(d, m)
44 {
45  nsmax = 15; // default maximum number of sweeps allowed
46  rtol = 10.E-12; // convergence tolerance
47 }
48
50
51
52 #define GJacobi_ZERO_CHECK_TOL 1.e-40
53
56 //
57 // this function solve the generalized eigenproblem using the Generalized
58 // jacobi iteration
59 //
60 //
61 {
62  int nsweep, nr;
63  double eps, eptola, eptolb, akk, ajj, ab, check, sqch, d1, d2, den, ca, cg, ak, bk, xj, xk;
64  double aj, bj;
65  int jm1, kp1, km1, jp1;
66
67  if ( a.giveNumberOfRows() != b.giveNumberOfRows() ||
68  !a.isSquare() || !b.isSquare() ) {
69  OOFEM_ERROR("A matrix, B mtrix -> size mismatch");
70  }
71
72  int n = a.giveNumberOfRows();
73  eigv.resize(n);
74  x.resize(n, n);
75
76  //
77  // Create temporary arrays
78  //
79  FloatArray d(n);
80  //
81  // Initialize EigenValue and EigenVector Matrices
82  //
83  for ( int i = 1; i <= n; i++ ) {
84  // if((a.at(i,i) <= 0. ) && (b.at(i,i) <= 0.))
85  // OOFEM_ERROR("Matrices are not positive definite");
86  d.at(i) = a.at(i, i) / b.at(i, i);
87  eigv.at(i) = d.at(i);
88  }
89
90  x.beUnitMatrix();
91
92  if ( n == 1 ) {
93  return NM_Success;
94  }
95
96  //
97  // Initialize sweep counter and begin iteration
98  //
99  nsweep = 0;
100  nr = n - 1;
101
102  do {
103  nsweep++;
104 # ifdef DETAILED_REPORT
105  OOFEM_LOG_DEBUG("*GJacobi*: sweep number %4d\n", nsweep);
106 #endif
107  //
108  // check if present off-diagonal element is large enough to require zeroing
109  //
110  eps = pow(0.01, ( double ) nsweep);
111  eps *= eps;
112  for ( int j = 1; j <= nr; j++ ) {
113  int jj = j + 1;
114  for ( int k = jj; k <= n; k++ ) {
115  eptola = ( a.at(j, k) * a.at(j, k) ) / ( a.at(j, j) * a.at(k, k) );
116  eptolb = ( b.at(j, k) * b.at(j, k) ) / ( b.at(j, j) * b.at(k, k) );
117  if ( ( eptola < eps ) && ( eptolb < eps ) ) {
118  continue;
119  }
120
121  //
122  // if zeroing is required, calculate the rotation matrix elements ca and cg
123  //
124  akk = a.at(k, k) * b.at(j, k) - b.at(k, k) * a.at(j, k);
125  ajj = a.at(j, j) * b.at(j, k) - b.at(j, j) * a.at(j, k);
126  ab = a.at(j, j) * b.at(k, k) - a.at(k, k) * b.at(j, j);
127  check = ( ab * ab + 4.0 * akk * ajj ) / 4.0;
128  if ( fabs(check) < GJacobi_ZERO_CHECK_TOL ) {
129  check = fabs(check);
130  } else if ( check < 0.0 ) {
131  OOFEM_ERROR("Matrices are not positive definite");
132  }
133
134  sqch = sqrt(check);
135  d1 = ab / 2. + sqch;
136  d2 = ab / 2. - sqch;
137  den = d1;
138  if ( fabs(d2) > fabs(d1) ) {
139  den = d2;
140  }
141
142  if ( den != 0.0 ) { // strange !
143  ca = akk / den;
144  cg = -ajj / den;
145  } else {
146  ca = 0.0;
147  cg = -a.at(j, k) / a.at(k, k);
148  }
149
150  //
151  // perform the generalized rotation to zero
152  //
153  if ( ( n - 2 ) != 0 ) {
154  jp1 = j + 1;
155  jm1 = j - 1;
156  kp1 = k + 1;
157  km1 = k - 1;
158  if ( ( jm1 - 1 ) >= 0 ) {
159  for ( int i = 1; i <= jm1; i++ ) {
160  aj = a.at(i, j);
161  bj = b.at(i, j);
162  ak = a.at(i, k);
163  bk = b.at(i, k);
164  a.at(i, j) = aj + cg * ak;
165  b.at(i, j) = bj + cg * bk;
166  a.at(i, k) = ak + ca * aj;
167  b.at(i, k) = bk + ca * bj;
168  }
169  }
170
171  if ( ( kp1 - n ) <= 0 ) {
172  for ( int i = kp1; i <= n; i++ ) { // label 140
173  aj = a.at(j, i);
174  bj = b.at(j, i);
175  ak = a.at(k, i);
176  bk = b.at(k, i);
177  a.at(j, i) = aj + cg * ak;
178  b.at(j, i) = bj + cg * bk;
179  a.at(k, i) = ak + ca * aj;
180  b.at(k, i) = bk + ca * bj;
181  }
182  }
183
184  if ( ( jp1 - km1 ) <= 0 ) { // label 160
185  for ( int i = jp1; i <= km1; i++ ) {
186  aj = a.at(j, i);
187  bj = b.at(j, i);
188  ak = a.at(i, k);
189  bk = b.at(i, k);
190  a.at(j, i) = aj + cg * ak;
191  b.at(j, i) = bj + cg * bk;
192  a.at(i, k) = ak + ca * aj;
193  b.at(i, k) = bk + ca * bj;
194  }
195  }
196  } // label 190
197
198  ak = a.at(k, k);
199  bk = b.at(k, k);
200  a.at(k, k) = ak + 2.0 *ca *a.at(j, k) + ca *ca *a.at(j, j);
201  b.at(k, k) = bk + 2.0 *ca *b.at(j, k) + ca *ca *b.at(j, j);
202  a.at(j, j) = a.at(j, j) + 2.0 *cg *a.at(j, k) + cg * cg * ak;
203  b.at(j, j) = b.at(j, j) + 2.0 *cg *b.at(j, k) + cg * cg * bk;
204  a.at(j, k) = 0.0;
205  b.at(j, k) = 0.0;
206  //
207  // update the eigenvector matrix after each rotation
208  //
209  for ( int i = 1; i <= n; i++ ) {
210  xj = x.at(i, j);
211  xk = x.at(i, k);
212  x.at(i, j) = xj + cg * xk;
213  x.at(i, k) = xk + ca * xj;
214  } // label 200
215  }
216  } // label 210
217
218  //
219  // update the eigenvalues after each sweep
220  //
221 #ifdef DETAILED_REPORT
222  OOFEM_LOG_DEBUG("GJacobi: a,b after sweep\n");
223  a.printYourself();
224  b.printYourself();
225 #endif
226  for ( int i = 1; i <= n; i++ ) {
227  // in original uncommented
228  // if ((a.at(i,i) <= 0.) || (b.at(i,i) <= 0.))
229  // error ("solveYourselfAt: Matrices are not positive definite");
230  eigv.at(i) = a.at(i, i) / b.at(i, i);
231  } // label 220
232
233 # ifdef DETAILED_REPORT
234  OOFEM_LOG_DEBUG("GJacobi: current eigenvalues are:\n");
235  eigv.printYourself();
236  OOFEM_LOG_DEBUG("GJacobi: current eigenvectors are:\n");
237  x.printYourself();
238 # endif
239  //
240  // check for convergence
241  //
242  for ( int i = 1; i <= n; i++ ) { // label 230
243  double tol = rtol * d.at(i);
244  double dif = ( eigv.at(i) - d.at(i) );
245  if ( fabs(dif) > tol ) {
246  goto label280;
247  }
248  } // label 240
249
250  //
251  // check all off-diagonal elements to see if another sweep is required
252  //
253  eps = rtol * rtol;
254  for ( int j = 1; j <= nr; j++ ) {
255  int jj = j + 1;
256  for ( int k = jj; k <= n; k++ ) {
257  double epsa = ( a.at(j, k) * a.at(j, k) ) / ( a.at(j, j) * a.at(k, k) );
258  double epsb = ( b.at(j, k) * b.at(j, k) ) / ( b.at(j, j) * b.at(k, k) );
259  if ( ( epsa < eps ) && ( epsb < eps ) ) {
260  continue;
261  }
262
263  goto label280;
264  }
265  } // label 250
266
267  //
268  // fill out bottom triangle of resultant matrices and scale eigenvectors
269  //
270  break;
271  // goto label255 ;
272  //
273  // update d matrix and start new sweep, if allowed
274  //
275 label280:
276  d = eigv;
277  } while ( nsweep < nsmax );
278
279  // label255:
280  for ( int i = 1; i <= n; i++ ) {
281  for ( int j = 1; j <= n; j++ ) {
282  a.at(j, i) = a.at(i, j);
283  b.at(j, i) = b.at(i, j);
284  } // label 260
285  }
286
287  for ( int j = 1; j <= n; j++ ) {
288  double bb = sqrt( fabs( b.at(j, j) ) );
289  for ( int k = 1; k <= n; k++ ) {
290  x.at(k, j) /= bb;
291  }
292  } // label 270
293
294  return NM_Success;
295 }
296 } // end namespace oofem
#define NM_Success
Numerical method exited with success.
Definition: nmstatus.h:47
Class and object Domain.
Definition: domain.h:115
#define GJacobi_ZERO_CHECK_TOL
Definition: gjacobi.C:52
double & at(int i)
Coefficient access function.
Definition: floatarray.h:131
This base class is an abstraction for numerical algorithm.
Definition: nummet.h:80
double rtol
Definition: gjacobi.h:57
bool isSquare() const
Returns nonzero if receiver is square matrix.
Definition: floatmatrix.h:160
unsigned long NM_Status
Mask defining NumMetod Status; which can be asked after finishing computation by Numerical Method...
Definition: nmstatus.h:44
#define OOFEM_LOG_DEBUG(...)
Definition: logger.h:128
virtual ~GJacobi()
Definition: gjacobi.C:49
GJacobi(Domain *d, EngngModel *m)
Definition: gjacobi.C:42
#define OOFEM_ERROR(...)
Definition: error.h:61
double at(int i, int j) const
Coefficient access function.
Definition: floatmatrix.h:176
Class representing vector of real numbers.
Definition: floatarray.h:82
Implementation of matrix containing floating point numbers.
Definition: floatmatrix.h:94
void resize(int rows, int cols)
Checks size of receiver towards requested bounds.
Definition: floatmatrix.C:1358
virtual void printYourself() const
Print receiver on stdout.
Definition: floatarray.C:747
void printYourself() const
Prints matrix to stdout. Useful for debugging.
Definition: floatmatrix.C:1458
virtual NM_Status solve(FloatMatrix &K, FloatMatrix &M, FloatArray &w, FloatMatrix &x)
Solves the given sparse generalized eigenvalue system of equations .
Definition: gjacobi.C:55
void beUnitMatrix()
Sets receiver to unity matrix.
Definition: floatmatrix.C:1332
Abstract base class representing the "problem" under consideration.
Definition: engngm.h:181
the oofem namespace is to define a context or scope in which all oofem names are defined.
int giveNumberOfRows() const
Returns number of rows of receiver.
Definition: floatmatrix.h:156
void resize(int s)
Resizes receiver towards requested size.
Definition: floatarray.C:631

This page is part of the OOFEM documentation. Copyright (c) 2011 Borek Patzak
Project e-mail: info@oofem.org
Generated at Tue Jan 2 2018 20:07:28 for OOFEM by doxygen 1.8.11 written by Dimitri van Heesch, © 1997-2011