1 /*--------------------------------------------------------------------------*/
2 /* ALBERTA: an Adaptive multi Level finite element toolbox using */
3 /* Bisectioning refinement and Error control by Residual */
4 /* Techniques for scientific Applications */
5 /* */
6 /* file: sor.c */
7 /* */
8 /* description: SSOR-method for scalar and decoupled vector valued */
9 /* problems */
10 /* */
11 /*--------------------------------------------------------------------------*/
12 /* */
13 /* authors: Alfred Schmidt */
14 /* Zentrum fuer Technomathematik */
15 /* Fachbereich 3 Mathematik/Informatik */
16 /* Univesitaet Bremen */
17 /* Bibliothekstr. 2 */
18 /* D-28359 Bremen, Germany */
19 /* */
20 /* Kunibert G. Siebert */
21 /* Institut fuer Mathematik */
22 /* Universitaet Augsburg */
23 /* Universitaetsstr. 14 */
24 /* D-86159 Augsburg, Germany */
25 /* */
26 /* Claus-Justus Heine */
27 /* Abteilung fuer Angewandte Mathematik */
28 /* Alberta-Ludwigs-Universitaet Freiburg */
29 /* Hermann-Herder-Str. 10 */
30 /* D-79104 Freiburg im Breisgau, Germany */
31 /* */
32 /* http://www.mathematik.uni-freiburg.de/IAM/ALBERTA */
33 /* */
34 /* (c) by A. Schmidt and K.G. Siebert (1996-2003) */
35 /* */
36 /*--------------------------------------------------------------------------*/
37
38 #include "alberta.h"
39
ssor_s(DOF_MATRIX * a,const DOF_REAL_VEC * f,const DOF_SCHAR_VEC * bound,DOF_REAL_VEC * u,REAL omega,REAL tol,int max_iter,int info)40 int ssor_s(DOF_MATRIX *a, const DOF_REAL_VEC *f, const DOF_SCHAR_VEC *bound,
41 DOF_REAL_VEC *u, REAL omega, REAL tol, int max_iter, int info)
42 {
43 FUNCNAME("ssor_s");
44 const REAL *fvec;
45 REAL *uvec;
46 S_CHAR *bvec;
47 REAL max = 0.0, omega1, accu, unew;
48 int i, iter, dim;
49 MATRIX_ROW *row;
50
51 fvec = f->vec;
52 uvec = u->vec;
53 bvec = bound ? bound->vec : NULL;
54
55 TEST_EXIT(a->row_fe_space->admin == a->col_fe_space->admin,
56 "Row and column FE_SPACEs don't match!\n");
57
58 if (a->row_fe_space->admin->hole_count > 0)
59 dof_compress(a->row_fe_space->mesh);
60
61 if (omega <= 0 || omega > 2)
62 {
63 ERROR("omega %le not in (0,2], setting omega = 1.0\n", omega);
64 omega = 1.0;
65 }
66 omega1 = 1.0 - omega;
67
68 if (info >= 2)
69 {
70 MSG("omega = %.3lf, tol = %.3le, max_iter = %d\n",
71 omega, tol, max_iter);
72 }
73
74 for (iter = 0; iter < max_iter; iter++)
75 {
76 max = 0.0;
77 dim = u->fe_space->admin->size_used;
78
79 for (i = 0; i < dim; i++)
80 {
81 if (!bvec || bvec[i] < DIRICHLET)
82 {
83 if (a->matrix_row[i])
84 {
85 accu = 0.0;
86 FOR_ALL_MAT_COLS(REAL, a->matrix_row[i], {
87 accu += row->entry[col_idx] * uvec[col_dof];
88 });
89 if ((row = a->matrix_row[i])) {
90 unew =
91 omega1 * uvec[i] + omega * (fvec[i] - accu) / row->entry.real[0];
92 } else {
93 unew = 0.0;
94 }
95 max = MAX(max, ABS(uvec[i] - unew));
96 uvec[i] = unew;
97 }
98 }
99 }
100
101 for (i = dim-1; i >= 0; i--)
102 {
103 if (!bvec || bvec[i] < DIRICHLET)
104 {
105 if (a->matrix_row[i])
106 {
107 accu = 0.0;
108 FOR_ALL_MAT_COLS(REAL, a->matrix_row[i], {
109 accu += row->entry[col_idx] * uvec[col_dof];
110 });
111 if ((row = a->matrix_row[i])) {
112 unew =
113 omega1 * uvec[i]
114 +
115 omega * (fvec[i] - accu) / row->entry.real[0];
116 } else {
117 unew = 0.0;
118 }
119 max = MAX(max, ABS(uvec[i] - unew));
120 uvec[i] = unew;
121 }
122 }
123 }
124 if (info >= 4)
125 MSG("iter %3d: max = %.3le\n",iter,max);
126
127 if (max < tol) break;
128 }
129
130 if (info >= 2)
131 {
132 if (iter < max_iter)
133 MSG("convergence after iter %3d: max = %.3le\n", iter, max);
134 else
135 MSG("NO CONVERGENCE after iter %3d: max = %.3le\n", iter, max);
136 }
137 return(iter);
138 }
139
ssor_d(DOF_MATRIX * a,const DOF_REAL_D_VEC * f,const DOF_SCHAR_VEC * bound,DOF_REAL_D_VEC * u,REAL omega,REAL tol,int max_iter,int info)140 int ssor_d(DOF_MATRIX *a,
141 const DOF_REAL_D_VEC *f, const DOF_SCHAR_VEC *bound,
142 DOF_REAL_D_VEC *u, REAL omega, REAL tol, int max_iter, int info)
143 {
144 FUNCNAME("ssor_d");
145 const REAL_D *fvec;
146 REAL_D *uvec;
147 S_CHAR *bvec;
148 REAL max = 0.0, omega1;
149 REAL_D accu, unew;
150 int n, i, iter = 0, dim;
151 MATRIX_ROW *row;
152
153 fvec = (const REAL_D *) f->vec;
154 uvec = u->vec;
155 bvec = bound ? bound->vec : NULL;
156
157 TEST_EXIT(a->row_fe_space->admin == a->col_fe_space->admin,
158 "Row and column FE_SPACEs don't match!\n");
159
160 if (a->row_fe_space->admin->hole_count > 0)
161 dof_compress(a->row_fe_space->mesh);
162
163 if (omega <= 0 || omega > 2) {
164 WARNING("omega %le not in (0,2], setting omega = 1.0\n", omega);
165 omega = 1.0;
166 }
167 omega1 = 1.0 - omega;
168
169 if (info >= 2) {
170 MSG("omega = %.3lf, tol = %.3le, max_iter = %d\n",
171 omega, tol, max_iter);
172 }
173
174 #undef MAT_BODY
175 #define MAT_BODY(F, CC, C, S, TYPE) \
176 for (iter = 0; iter < max_iter; iter++) \
177 { \
178 max = 0.0; \
179 dim = u->fe_space->admin->size_used; \
180 \
181 for (i = 0; i < dim; i++) { \
182 if ((row = a->matrix_row[i]) == NULL) { \
183 continue; \
184 } \
185 \
186 if (!bvec || bvec[i] < DIRICHLET) { \
187 COPY_DOW(fvec[i], accu); \
188 \
189 FOR_ALL_MAT_COLS(TYPE, a->matrix_row[i], { \
190 if (col_dof != i) { \
191 F##GEMV_DOW(-1.0, CC row->entry[col_idx], uvec[col_dof], \
192 1.0, accu); \
193 } else { \
194 F##GEMV_ND_DOW(-1.0, CC row->entry[col_idx], uvec[col_dof], \
195 1.0, accu); \
196 } \
197 }); \
198 F##DIV_DOW(CC row->entry.S[0], accu, accu); \
199 AXPBY_DOW(omega, accu, omega1, uvec[i], unew); \
200 \
201 for (n = 0; n < DIM_OF_WORLD; n++) { \
202 max = MAX(max, ABS(uvec[i][n] - unew[n])); \
203 uvec[i][n] = unew[n]; \
204 } \
205 } \
206 } \
207 \
208 for (i = dim-1; i >= 0; i--) { \
209 if ((row = a->matrix_row[i]) == NULL) { \
210 continue; \
211 } \
212 \
213 if (!bvec || bvec[i] < DIRICHLET) { \
214 COPY_DOW(fvec[i], accu); \
215 \
216 FOR_ALL_MAT_COLS(TYPE, a->matrix_row[i], { \
217 if (col_dof != i) { \
218 F##GEMV_DOW(-1.0, CC row->entry[col_idx], uvec[col_dof], \
219 1.0, accu); \
220 } else { \
221 F##GEMV_ND_DOW(-1.0, CC row->entry[col_idx], uvec[col_dof], \
222 1.0, accu); \
223 } \
224 }); \
225 F##DIV_DOW(CC row->entry.S[0], accu, accu); \
226 AXPBY_DOW(omega, accu, omega1, uvec[i], unew); \
227 \
228 for (n = 0; n < DIM_OF_WORLD; n++) { \
229 max = MAX(max, ABS(uvec[i][n] - unew[n])); \
230 uvec[i][n] = unew[n]; \
231 } \
232 } \
233 } \
234 if (info >= 4) \
235 MSG("iter %3d: max = %.3le\n",iter,max); \
236 \
237 if (max < tol) break; \
238 }
239
240 MAT_EMIT_BODY_SWITCH(a->type);
241
242 if (info >= 2) {
243 if (iter < max_iter)
244 MSG("convergence after iter %3d: max = %.3le\n", iter, max);
245 else
246 MSG("NO CONVERGENCE after iter %3d: max = %.3le\n", iter, max);
247 }
248 return iter;
249 }
250