Line data Source code
1 : /*
2 : * Copyright (c) 2010-2015: G-CSC, Goethe University Frankfurt
3 : * Author: Daniel Jungblut
4 : *
5 : * This file is part of UG4.
6 : *
7 : * UG4 is free software: you can redistribute it and/or modify it under the
8 : * terms of the GNU Lesser General Public License version 3 (as published by the
9 : * Free Software Foundation) with the following additional attribution
10 : * requirements (according to LGPL/GPL v3 §7):
11 : *
12 : * (1) The following notice must be displayed in the Appropriate Legal Notices
13 : * of covered and combined works: "Based on UG4 (www.ug4.org/license)".
14 : *
15 : * (2) The following notice must be displayed at a prominent place in the
16 : * terminal output of covered works: "Based on UG4 (www.ug4.org/license)".
17 : *
18 : * (3) The following bibliography is recommended for citation and must be
19 : * preserved in all covered files:
20 : * "Reiter, S., Vogel, A., Heppner, I., Rupp, M., and Wittum, G. A massively
21 : * parallel geometric multigrid solver on hierarchically distributed grids.
22 : * Computing and visualization in science 16, 4 (2013), 151-164"
23 : * "Vogel, A., Reiter, S., Rupp, M., Nägel, A., and Wittum, G. UG4 -- a novel
24 : * flexible software system for simulating pde based models on high performance
25 : * computers. Computing and visualization in science 16, 4 (2013), 165-179"
26 : *
27 : * This program is distributed in the hope that it will be useful,
28 : * but WITHOUT ANY WARRANTY; without even the implied warranty of
29 : * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
30 : * GNU Lesser General Public License for more details.
31 : */
32 :
33 : #include <cmath>
34 : #include "../ugmath_types.h"
35 : #include "eigenvalues.h"
36 :
37 : namespace ug
38 : {
39 :
40 : static void rot(number A[3][3], const number s, const number tau,
41 : const int i, const int j, const int k, const int l)
42 : {
43 :
44 0 : number g = A[i][j];
45 0 : number h = A[k][l];
46 :
47 0 : A[i][j] = g - s * (h + g * tau);
48 0 : A[k][l] = h + s * (g - h * tau);
49 : }
50 :
51 : // For symmetrische Matrizen:
52 0 : bool CalculateEigenvalues(const ug::matrix33& mat, number& lambdaMinOut,
53 : number& lambdaMedOut, number& lambdaMaxOut,
54 : ug::vector3& evMinOut, ug::vector3& evMedOut,
55 : ug::vector3& evMaxOut)
56 : {
57 :
58 : number A[3][3];
59 :
60 : int i, j, ip, iq;
61 : number tresh, theta, tau, t, sm, s, h, g, c;
62 :
63 : int n = 3;
64 :
65 : number b[3];
66 : number d[3];
67 : number z[3];
68 :
69 : number V[3][3];
70 :
71 : // v wird als Einheitsmatrix initialisiert:
72 : // mat wird in A wird kopiert
73 0 : for(ip = 0; ip < 3; ip++) {
74 0 : for(iq = 0; iq < 3; iq++) {
75 0 : A[ip][iq] = mat[ip][iq];
76 0 : V[ip][iq] = 0.0f;
77 : }
78 0 : V[ip][ip] = 1.0f;
79 : }
80 :
81 : // b und d werden mit der Diagonale von A initialisiert:
82 : // z wird mit 0 initialisiert:
83 0 : for(ip = 0; ip < 3; ip++) {
84 0 : b[ip] = d[ip] = A[ip][ip];
85 0 : z[ip] = 0;
86 : }
87 :
88 : //int nrot = 0;
89 :
90 0 : for(i = 1; i <= 50; i++) {
91 : sm = 0.0f;
92 :
93 0 : for(ip = 0; ip < 2; ip++)
94 0 : for(iq = ip+1; iq < 3; iq ++)
95 0 : sm += fabs(A[ip][iq]);
96 :
97 0 : if(fabs(sm) <= 0.00001f) {
98 :
99 : // Berechnung fertig:
100 :
101 : // Code zum sortieren der Eigenwerte einf�gen:
102 : int max_index = 0;
103 : int med_index = 0;
104 : int min_index = 0;
105 : number max = 0.0f;
106 :
107 0 : for(int ii = 0; ii < 3; ii++)
108 0 : if(fabs(d[ii]) > max) {
109 : max_index = ii;
110 : max = fabs(d[ii]);
111 : }
112 :
113 0 : switch(max_index) {
114 :
115 0 : case 0: if(fabs(d[1]) > fabs(d[2])) {med_index = 1; min_index = 2;} else {med_index = 2; min_index = 1;} break;
116 0 : case 1: if(fabs(d[0]) > fabs(d[2])) {med_index = 0; min_index = 2;} else {med_index = 2; min_index = 0;} break;
117 0 : case 2: if(fabs(d[0]) > fabs(d[1])) {med_index = 0; min_index = 1;} else {med_index = 1; min_index = 0;} break;
118 : default: break;
119 :
120 : }
121 :
122 : // Gr��ter Eigenwert:
123 0 : lambdaMaxOut = d[max_index];
124 0 : evMaxOut[0] = V[0][max_index];
125 0 : evMaxOut[1] = V[1][max_index];
126 0 : evMaxOut[2] = V[2][max_index];
127 :
128 : // Mittlerer Eigenwert:
129 0 : lambdaMedOut = d[med_index];
130 0 : evMedOut[0] = V[0][med_index];
131 0 : evMedOut[1] = V[1][med_index];
132 0 : evMedOut[2] = V[2][med_index];
133 :
134 : // Kleinster Eigenwert:
135 0 : lambdaMinOut = d[min_index];
136 0 : evMinOut[0] = V[0][min_index];
137 0 : evMinOut[1] = V[1][min_index];
138 0 : evMinOut[2] = V[2][min_index];
139 :
140 0 : return true; // Fertig!!!
141 :
142 : }
143 :
144 :
145 0 : if(i < 4)
146 0 : tresh = 0.2f * sm / 9.0f;
147 :
148 : else
149 : tresh = 0.0f;
150 :
151 0 : for(ip = 0; ip < (n-1); ip++) {
152 0 : for(iq = ip+1; iq < 3; iq++) {
153 0 : g = 100.0f * fabs(A[ip][iq]);
154 :
155 0 : if((i > 4) && ((fabs(d[ip])+g) == fabs(d[ip])) && ((fabs(d[iq])+g) == fabs(d[iq])))
156 0 : A[ip][iq] = 0.0f;
157 :
158 0 : else if(fabs(A[ip][iq]) > tresh) {
159 0 : h = d[iq] - d[ip];
160 :
161 0 : if((fabs(h)+g) == fabs(h))
162 0 : t = A[ip][iq] / h;
163 : else {
164 0 : theta = 0.5f * h / A[ip][iq];
165 0 : t = 1.0f / (fabs(theta) + sqrt(1.0f + theta * theta));
166 0 : if(theta < 0.0f) t *= -1.0f;
167 : }
168 :
169 0 : c = 1.0f / sqrt(1.0f + t * t);
170 0 : s = t * c;
171 0 : tau = s / (1.0f + c);
172 0 : h = t * A[ip][iq];
173 0 : z[ip] -= h;
174 0 : z[iq] += h;
175 0 : d[ip] -= h;
176 0 : d[iq] += h;
177 0 : A[ip][iq] = 0.0f;
178 :
179 0 : for(j = 0; j < ip; j++)
180 : rot(A, s, tau, j, ip, j, iq);
181 :
182 0 : for(j = ip + 1; j < iq; j++)
183 : rot(A, s, tau, ip, j, j, iq);
184 :
185 0 : for(j = iq + 1; j < 3; j++)
186 : rot(A, s, tau, ip, j, iq, j);
187 :
188 0 : for(j = 0; j < 3; j++)
189 : rot(V, s, tau, j, ip, j, iq);
190 :
191 : //nrot ++;
192 :
193 : }
194 : }
195 : }
196 :
197 0 : for(ip = 0; ip < 3; ip++) {
198 0 : b[ip] += z[ip];
199 0 : d[ip] = b[ip];
200 0 : z[ip] = 0.0f;
201 : }
202 :
203 : }
204 :
205 : return false;
206 : }
207 :
208 : }// end of namespace
209 :
|
