/[ascend]/trunk/base/generic/solver/example/example2.c
ViewVC logotype

Contents of /trunk/base/generic/solver/example/example2.c

Parent Directory Parent Directory | Revision Log Revision Log


Revision 525 - (show annotations) (download) (as text)
Fri Apr 21 10:41:45 2006 UTC (18 years, 1 month ago) by johnpye
File MIME type: text/x-csrc
File size: 7014 byte(s)
Rolled 'ascPrintType.h' into 'ascPrint.h' for simplicity.
1 #include <utilities/ascConfig.h>
2 #include <solver/mtx.h>
3 #include <solver/linsol.h>
4
5 #ifdef __WIN32__
6 #include <utilities/ascPrint.h>
7 #define f_vtable_name "example_vtable"
8 static struct Asc_PrintVTable f_vtable = {f_vtable_name, vfprintf, fflush, NULL};
9 #endif
10
11 static linsol_system_t sys;
12 static mtx_matrix_t N, NTN;
13 static real64 D1[10], D2[10];
14
15 static make_coefficient_matrix_1()
16 {
17 mtx_coord_t coord;
18
19 printf("Create N\n");
20 N = mtx_create();
21 mtx_set_order(N,10);
22
23 /* eqn 0 */
24
25 coord.row = 0;
26 coord.col = 0;
27 mtx_set_value(N,&coord,1.);
28
29 /* eqn 1 */
30
31 coord.row = 1;
32 coord.col = 0;
33 mtx_set_value(N,&coord,0.007);
34 coord.col = 1;
35 mtx_set_value(N,&coord,0.571);
36 coord.col = 2;
37 mtx_set_value(N,&coord,0.422);
38
39 /* eqn 2 */
40
41 coord.row = 2;
42 coord.col = 1;
43 mtx_set_value(N,&coord,0.125);
44 coord.col = 2;
45 mtx_set_value(N,&coord,0.75);
46 coord.col = 3;
47 mtx_set_value(N,&coord,0.125);
48
49 /* eqn 3 */
50
51 coord.row = 3;
52 coord.col = 2;
53 mtx_set_value(N,&coord,0.422);
54 coord.col = 3;
55 mtx_set_value(N,&coord,0.571);
56 coord.col = 4;
57 mtx_set_value(N,&coord,0.007);
58
59 /* eqn 4 */
60
61 coord.row = 4;
62 coord.col = 4;
63 mtx_set_value(N,&coord,1.);
64 }
65
66 static make_coefficient_matrix_2()
67 {
68 mtx_coord_t coord;
69
70 printf("Create N\n");
71 N = mtx_create();
72 mtx_set_order(N,10);
73
74 /* eqn 0 */
75
76 coord.row = 0;
77 coord.col = 0;
78 mtx_set_value(N,&coord,1.0);
79
80 /* eqn 1 */
81
82 coord.row = 1;
83 coord.col = 0;
84 mtx_set_value(N,&coord,0.15);
85 coord.col = 1;
86 mtx_set_value(N,&coord,0.662);
87 coord.col = 2;
88 mtx_set_value(N,&coord,0.188);
89
90 /* eqn 2 */
91
92 coord.row = 2;
93 coord.col = 1;
94 mtx_set_value(N,&coord,0.5);
95 coord.col = 2;
96 mtx_set_value(N,&coord,0.5);
97
98 /* eqn 3 */
99
100 coord.row = 3;
101 coord.col = 1;
102 mtx_set_value(N,&coord,0.188);
103 coord.col = 2;
104 mtx_set_value(N,&coord,0.662);
105 coord.col = 3;
106 mtx_set_value(N,&coord,0.15);
107
108 /* eqn 4 */
109
110 coord.row = 4;
111 coord.col = 3;
112 mtx_set_value(N,&coord,1.0);
113 }
114
115 static make_rhs()
116 {
117 printf("\tD1 = { 0.0 } D2 = { 0.0 }\n"
118 "\t { 1.5 } { 2.0 }\n"
119 "\t { 3.0 } { 2.5 }\n"
120 "\t { 4.5 } { 2.0 }\n"
121 "\t { 6.0 } { 0.0 }\n");
122 D1[0] = 0.0;
123 D1[1] = 1.5;
124 D1[2] = 3.0;
125 D1[3] = 4.5;
126 D1[4] = 6.0;
127
128 linsol_add_rhs(sys,D1,FALSE);
129
130 D2[0] = 0.0;
131 D2[1] = 2.0;
132 D2[2] = 2.5;
133 D2[3] = 2.0;
134 D2[4] = 0.0;
135
136 linsol_add_rhs(sys,D2,FALSE);
137 }
138
139 mtx_matrix_t ATA(mtx_matrix_t N)
140 {
141 mtx_matrix_t NTN;
142 int32 i;
143 int32 v1,v2,v3;
144 mtx_coord_t c1,c2,c3,c4;
145
146 NTN = mtx_create();
147 mtx_set_order(NTN,mtx_order(N));
148
149 for (i = 0; i < mtx_order(N); i++) {
150 c1.row = i;
151 c1.col = mtx_FIRST;
152 c2.row = c1.row;
153 while (mtx_next_in_row(N,&c1,mtx_ALL_COLS)) {
154 v1 = mtx_value(N,&c1);
155 c2.col = mtx_FIRST;
156 while (mtx_next_in_row(N,&c2,mtx_ALL_COLS)) {
157 if (c2.col >= c1.col) {
158 v2 = v1*mtx_value(N,&c2);
159 c3.row = c2.col;
160 c3.col = c1.col;
161 v3 = v2 + mtx_value(NTN,&c3);
162 mtx_set_value(NTN,&c3,v3);
163 if (c3.row != c3.col) {
164 c4.row = c3.col;
165 c4.col = c3.row;
166 v3 = v2 + mtx_value(NTN,&c4);
167 mtx_set_value(NTN,&c4,v3);
168 }
169 }
170 }
171 }
172 }
173 return(NTN);
174 }
175
176 static ATv(mtx_matrix_t A, real64 v[], real64 r[])
177 {
178 mtx_coord_t c;
179 int32 i;
180
181 for (i = 0; i < mtx_order(A); i++) r[i] = 0.0;
182
183 for (i = 0; i < mtx_order(A); i++) {
184 c.row = i;
185 c.col = mtx_FIRST;
186 while (mtx_next_in_row(A,&c,mtx_ALL_COLS)) {
187 r[c.col] = r[c.col] + mtx_value(A,&c)*v[i];
188 }
189 }
190 }
191
192 /* print out non-zero elements of a matrix in text by row */
193 static print_matrix(mtx_matrix_t N, mtx_region_t region)
194 {
195 mtx_coord_t c;
196 int32 i;
197 int has_row;
198
199 for (i = region.col.low; i <= region.col.high; i++) {
200 c.row = i;
201 c.col = mtx_FIRST;
202 has_row = FALSE;
203 while (mtx_next_in_row(N,&c,mtx_ALL_COLS)) {
204 has_row = TRUE;
205 printf("\tmtx[%d][%d] = %f\n",c.row,c.col,mtx_value(N,&c));
206 }
207 if (has_row) printf("\n");
208 }
209 }
210
211 static print_vector(real64 *v, mtx_range_t range)
212 {
213 int32 i;
214
215 for (i = range.low; i <= range.high; i++) {
216 printf("\tvec[%d] = %f\n",i,v[i]);
217 }
218 }
219
220 static print_solution(real64 *rhs, mtx_range_t range)
221 {
222 int32 i;
223
224 printf("Solution:\n");
225 for (i = range.low; i <= range.high; i++) {
226 printf("\ts[%d] = %f\n",i,linsol_var_value(sys,rhs,i));
227 }
228 }
229
230 static print_residuals(real64 *rhs, mtx_range_t range)
231 {
232 int32 i;
233
234 printf("Residuals:\n");
235 for (i = range.low; i <= range.high; i++) {
236 printf("\tr[%d] = %f\n",i,linsol_eqn_residual(sys,rhs,i));
237 }
238 }
239
240 int main()
241 {
242 real64 rhs[10];
243 mtx_range_t range;
244 mtx_region_t region = {0,4,0,4};
245
246 #ifdef __WIN32__
247 Asc_PrintPushVTable(&f_vtable);
248 #endif
249
250 sys = linsol_create();
251
252 printf("\n1st Case - Square\n");
253
254 printf("Make coefficient matrix N\n");
255 make_coefficient_matrix_1();
256 print_matrix(N, region);
257 printf("Make right hand side(s) D1 & D2\n");
258 make_rhs();
259
260 printf("Number of right hand sides = %d\n"
261 ,linsol_number_of_rhs(sys));
262
263 printf("Set coefficient matrix\n");
264 linsol_set_matrix(sys,N);
265 printf("Reorder equations\n");
266 linsol_reorder(sys,&region);
267
268 printf("Invert matrix\n");
269 linsol_invert(sys,&region);
270
271 printf("Solve D1\n");
272 linsol_solve(sys,D1);
273 printf("Rank = %d\n",linsol_rank(sys));
274 range.low = 0;
275 range.high = 4;
276 print_residuals(D1,range);
277 print_solution(D1,range);
278
279 printf("Solve D2\n");
280 linsol_solve(sys,D2);
281 printf("Rank = %d\n",linsol_rank(sys));
282 print_residuals(D2,range);
283 print_solution(D2,range);
284
285 mtx_destroy(N);
286 linsol_destroy(sys);
287
288 /* 2nd case */
289
290 printf("\n2nd case - overspecified\n");
291 printf("Make coefficient matrix N\n");
292 make_coefficient_matrix_2();
293 print_matrix(N, region);
294 printf("Computing NTN\n");
295 NTN = ATA(N);
296 print_matrix(NTN, region);
297 printf("Compute rhs\n");
298 ATv(N,D1,rhs);
299 range.low = 0;
300 range.high = 3;
301 print_vector(rhs,range);
302 sys = linsol_create();
303 linsol_set_matrix(sys,NTN);
304 linsol_add_rhs(sys,rhs,FALSE);
305 printf("Reorder matrix\n");
306 linsol_reorder(sys, &region);
307 printf("Invert matrix\n");
308 linsol_invert(sys, &region);
309 printf("Solve with inverse rhs\n");
310 linsol_solve(sys,rhs);
311
312 printf("Rank = %d\n",linsol_rank(sys));
313 print_residuals(rhs,range);
314 print_solution(rhs,range);
315
316 printf("Solve with inverse rhs\n");
317 ATv(N,D2,rhs);
318 print_vector(rhs,range);
319 linsol_rhs_was_changed(sys,rhs);
320 linsol_solve(sys,rhs);
321 printf("Rank = %d\n",linsol_rank(sys));
322 print_residuals(rhs,range);
323 print_solution(rhs,range);
324
325 #ifdef __WIN32__
326 Asc_PrintRemoveVTable(f_vtable_name);
327 #endif
328
329 return 0;
330 }

john.pye@anu.edu.au
ViewVC Help
Powered by ViewVC 1.1.22