1 /****************************************************************************/
2 /* extrsolv.c */
3 /****************************************************************************/
4 /* */
5 /* EXTeRn SOLVers for systems of linear equations */
6 /* */
7 /* Copyright (C) 1992-1995 Tomas Skalicky. All rights reserved. */
8 /* */
9 /****************************************************************************/
10 /* */
11 /* ANY USE OF THIS CODE CONSTITUTES ACCEPTANCE OF THE TERMS */
12 /* OF THE COPYRIGHT NOTICE (SEE FILE COPYRGHT.H) */
13 /* */
14 /****************************************************************************/
15
16 #include <math.h>
17 #include <stdio.h>
18
19 #include <laspack/errhandl.h>
20 #include <laspack/operats.h>
21 #include <laspack/precond.h>
22 #include <laspack/rtc.h>
23 #include <laspack/itersolv.h>
24 #include <laspack/copyrght.h>
25
26 #include "extrsolv.h"
27 #include "mlstest.h"
28
29 extern double RTCEps;
30
31 #define XSOR
32
33 #ifdef SYM_STOR
34
TestIter(QMatrix * a,Vector * x,Vector * b,int maxit,PrecondProcType dummy,double rlx)35 Vector *TestIter(QMatrix *a, Vector *x, Vector *b, int maxit,
36 PrecondProcType dummy, double rlx)
37 {
38 LASBreak();
39 fprintf(stderr, "\n");
40 fprintf(stderr, "No test iterative method available for nonsymmetric matrices.\n");
41 fprintf(stderr, "\n");
42
43 return(x);
44 }
45
46 #else
47
48 #ifdef XSOR
49
TestIter(QMatrix * a,Vector * x,Vector * b,int maxit,PrecondProcType dummy,double rlx)50 Vector *TestIter(QMatrix *a, Vector *x, Vector *b, int maxit,
51 PrecondProcType dummy, double rlx)
52 /* SOR lexicographicaly forward iteration ... implementation by Andreas Auge */
53 {
54 int it;
55 size_t i, j, col, dim;
56 double x_old_i, sum_diff_2,
57 sum_ax, sum_b_2 = 0.0, diag_el;
58
59 dim = a->Dim;
60
61 if (LASResult() == LASOK) {
62 for (i = 1; i <= dim; i++)
63 sum_b_2 += b->Cmp[i] * b->Cmp[i];
64
65 it = 0;
66 do {
67 sum_diff_2 = 0.0;
68
69 for (i = 1; i <= dim; i++) {
70 x_old_i = x->Cmp[i];
71 sum_ax = 0.0;
72 diag_el = 0.0;
73
74 for (j = 0; j < a->Len[i]; j++) {
75 col = a->El[i][j].Pos;
76 if (col != i) {
77 sum_ax += a->El[i][j].Val * x->Cmp[col];
78 } else {
79 diag_el = a->El[i][j].Val;
80 } /*else*/
81 } /*j*/
82 x->Cmp[i] = rlx * ((b->Cmp[i] - sum_ax) / diag_el - x_old_i) + x_old_i;
83
84 } /*i*/
85
86 for (i = 1; i <= dim; i++) {
87 sum_ax = 0.0;
88
89 for (j = 0; j < a->Len[i]; j++) {
90 sum_ax += a->El[i][j].Val * x->Cmp[a->El[i][j].Pos];
91 } /*j*/
92 sum_diff_2 += (b->Cmp[i] - sum_ax) * (b->Cmp[i] - sum_ax);
93 } /*i*/
94
95 it++;
96 } while (!RTCResult(it, sqrt(sum_diff_2), sqrt(sum_b_2), SORForwIterId + 100)
97 && it < maxit);
98 }
99
100 return(x);
101 }
102
103 #endif /* XSOR */
104
105 #ifdef XSSOR
106
TestIter(QMatrix * a,Vector * x,Vector * b,int maxit,PrecondProcType dummy,double rlx)107 Vector *TestIter(QMatrix *a, Vector *x, Vector *b, int maxit,
108 PrecondProcType dummy, double rlx)
109 /* SSOR iteration ... implementation by Andreas Auge */
110 {
111 int it;
112 size_t i, j, col, dim;
113 double x_old_i, sum_diff_2,
114 sum_ax, sum_b_2 = 0.0, diag_el;
115
116 dim = a->Dim;
117
118 if (LASResult() == LASOK) {
119 for (i = 1; i <= dim; i++)
120 sum_b_2 += b->Cmp[i] * b->Cmp[i];
121
122 it = 0;
123 do {
124 sum_diff_2 = 0.0;
125
126 for (i = 1; i <= dim; i++) {
127 x_old_i = x->Cmp[i];
128 sum_ax = 0.0;
129 diag_el = 0.0;
130
131 for (j = 0; j < a->Len[i]; j++) {
132 col = a->El[i][j].Pos;
133 if (col != i) {
134 sum_ax += a->El[i][j].Val * x->Cmp[col];
135 } else {
136 diag_el = a->El[i][j].Val;
137 } /*else*/
138 } /*j*/
139 x->Cmp[i] = rlx * ((b->Cmp[i] - sum_ax) / diag_el - x_old_i) + x_old_i;
140
141 } /*i*/
142
143 for (i = dim; i >= 1; i--) {
144 x_old_i = x->Cmp[i];
145 sum_ax = 0.0;
146 diag_el = 0.0;
147
148 for (j = 0; j < a->Len[i]; j++) {
149 col = a->El[i][j].Pos;
150 if (col != i) {
151 sum_ax += a->El[i][j].Val * x->Cmp[col];
152 } else {
153 diag_el = a->El[i][j].Val;
154 } /*else*/
155 } /*j*/
156 x->Cmp[i] = rlx * ((b->Cmp[i] - sum_ax) / diag_el - x_old_i) + x_old_i;
157
158 } /*i*/
159
160 for (i = 1; i <= dim; i++) {
161 sum_ax = 0.0;
162
163 for (j = 0; j < a->Len[i]; j++) {
164 sum_ax += a->El[i][j].Val * x->Cmp[a->El[i][j].Pos];
165 } /*j*/
166 sum_diff_2 += (b->Cmp[i] - sum_ax) * (b->Cmp[i] - sum_ax);
167 } /*i*/
168
169 it++;
170 } while (!RTCResult(it, sqrt(sum_diff_2), sqrt(sum_b_2), SSORIterId + 100)
171 && it < maxit);
172 }
173
174 return(x);
175 }
176
177 #endif /* XSSOR */
178
179 #endif /* SYM_STOR */
180