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

Contents of /trunk/ascend4/solver/example/example2.c

Parent Directory Parent Directory | Revision Log Revision Log


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

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