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

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

Parent Directory Parent Directory | Revision Log Revision Log | View Patch Patch

revision 60 by aw0a, Sat Nov 13 16:45:56 2004 UTC revision 61 by jds, Mon Nov 14 02:37:20 2005 UTC
# Line 1  Line 1 
1  #include "base.h"  #include "utilities/ascConfig.h"
2  #include "linsol.h"  #include "solver/mtx.h"
3    #include "solver/linsol.h"
4    
5    #ifdef __WIN32__
6    #include "utilities/ascPrintType.h"
7    #include "utilities/ascPrint.h"
8    #define f_vtable_name "example_vtable"
9    static struct Asc_PrintVTable f_vtable = {f_vtable_name, vfprintf, fflush, NULL};
10    #endif
11    
12  static linsol_system_t sys;  static linsol_system_t sys;
13  static mtx_matrix_t N,NTN;  static mtx_matrix_t N, NTN;
14  static mtx_number_t D1[10],D2[10];  static real64 D1[10], D2[10];
15    
16  static make_coefficient_matrix_1()  static make_coefficient_matrix_1()
17  {  {
# Line 59  static make_coefficient_matrix_1() Line 67  static make_coefficient_matrix_1()
67  static make_coefficient_matrix_2()  static make_coefficient_matrix_2()
68  {  {
69     mtx_coord_t coord;     mtx_coord_t coord;
70    
71     printf("Create N\n");     printf("Create N\n");
72     N = mtx_create();     N = mtx_create();
73     mtx_set_order(N,10);     mtx_set_order(N,10);
# Line 107  static make_coefficient_matrix_2() Line 115  static make_coefficient_matrix_2()
115    
116  static make_rhs()  static make_rhs()
117  {  {
118       printf("\tD1 = { 0.0 }   D2 = { 0.0 }\n"
119              "\t     { 1.5 }        { 2.0 }\n"
120              "\t     { 3.0 }        { 2.5 }\n"
121              "\t     { 4.5 }        { 2.0 }\n"
122              "\t     { 6.0 }        { 0.0 }\n");
123     D1[0] = 0.0;     D1[0] = 0.0;
124     D1[1] = 1.5;     D1[1] = 1.5;
125     D1[2] = 3.0;     D1[2] = 3.0;
126     D1[3] = 4.5;     D1[3] = 4.5;
127     D1[4] = 6.0;     D1[4] = 6.0;
128    
129     linsol_add_rhs(sys,D1);     linsol_add_rhs(sys,D1,FALSE);
130    
131     D2[0] = 0.0;     D2[0] = 0.0;
132     D2[1] = 2.0;     D2[1] = 2.0;
# Line 121  static make_rhs() Line 134  static make_rhs()
134     D2[3] = 2.0;     D2[3] = 2.0;
135     D2[4] = 0.0;     D2[4] = 0.0;
136    
137     linsol_add_rhs(sys,D2);     linsol_add_rhs(sys,D2,FALSE);
138  }  }
139    
140  mtx_matrix_t ATA(N)  mtx_matrix_t ATA(mtx_matrix_t N)
 mtx_matrix_t N;  
141  {  {
142      mtx_matrix_t NTN;     mtx_matrix_t NTN;
143      mtx_index_t i;     int32 i;
144      mtx_number_t v1,v2,v3;     int32 v1,v2,v3;
145      mtx_coord_t c1,c2,c3,c4;     mtx_coord_t c1,c2,c3,c4;
146        
147       NTN = mtx_create();
148      NTN = mtx_create();     mtx_set_order(NTN,mtx_order(N));
149      mtx_set_order(NTN,mtx_order(N));  
150       for (i = 0; i < mtx_order(N); i++) {
151      for (i = 0; i < mtx_order(N); i++) {        c1.row = i;
152         c1.row = i;        c1.col = mtx_FIRST;
153         c1.col = mtx_FIRST;        c2.row = c1.row;
154         c2.row = c1.row;        while (mtx_next_in_row(N,&c1,mtx_ALL_COLS)) {
155         while (mtx_next_in_row(N,&c1,mtx_ALL_COLS)) {           v1 = mtx_value(N,&c1);
156            v1 = mtx_value(N,&c1);           c2.col = mtx_FIRST;
157        c2.col = mtx_FIRST;           while (mtx_next_in_row(N,&c2,mtx_ALL_COLS)) {
158        while (mtx_next_in_row(N,&c2,mtx_ALL_COLS)) {              if (c2.col >= c1.col) {
159           if (c2.col >= c1.col) {                 v2 = v1*mtx_value(N,&c2);
160              v2 = v1*mtx_value(N,&c2);                 c3.row = c2.col;
161              c3.row = c2.col;                 c3.col = c1.col;
162              c3.col = c1.col;                 v3 = v2 + mtx_value(NTN,&c3);
163              v3 = v2 + mtx_value(NTN,&c3);                 mtx_set_value(NTN,&c3,v3);
164          mtx_set_value(NTN,&c3,v3);                 if (c3.row != c3.col) {
165          if (c3.row != c3.col) {                    c4.row = c3.col;
166             c4.row = c3.col;                    c4.col = c3.row;
167             c4.col = c3.row;                    v3 = v2 + mtx_value(NTN,&c4);
168             v3 = v2 + mtx_value(NTN,&c4);                    mtx_set_value(NTN,&c4,v3);
169             mtx_set_value(NTN,&c4,v3);                 }
170           }              }
171            }           }
172        }          }
173         }     }
     }  
174     return(NTN);     return(NTN);
175  }  }
176    
177  static ATv(A,v,r)  static ATv(mtx_matrix_t A, real64 v[], real64 r[])
 mtx_matrix_t A;  
 mtx_number_t v[],r[];  
178  {  {
179     mtx_coord_t c;     mtx_coord_t c;
180     mtx_index_t i;     int32 i;
181    
182     for (i = 0; i < mtx_order(A); i++) r[i] = 0.0;     for (i = 0; i < mtx_order(A); i++) r[i] = 0.0;
183    
# Line 176  mtx_number_t v[],r[]; Line 185  mtx_number_t v[],r[];
185        c.row = i;        c.row = i;
186        c.col = mtx_FIRST;        c.col = mtx_FIRST;
187        while (mtx_next_in_row(A,&c,mtx_ALL_COLS)) {        while (mtx_next_in_row(A,&c,mtx_ALL_COLS)) {
188        r[c.col] = r[c.col] + mtx_value(A,&c)*v[i];           r[c.col] = r[c.col] + mtx_value(A,&c)*v[i];
189        }        }
190     }     }
191  }  }
192    
193  static print_matrix(N)  /* print out non-zero elements of a matrix in text by row */
194  mtx_matrix_t N;  static print_matrix(mtx_matrix_t N, mtx_region_t region)
195  {  {
196     mtx_coord_t c;     mtx_coord_t c;
197     mtx_index_t i;     int32 i;
198       int has_row;
199    
200     for (i = 0; i < mtx_order(N); i++) {     for (i = region.col.low; i <= region.col.high; i++) {
201        c.row = i;        c.row = i;
202        c.col = mtx_FIRST;        c.col = mtx_FIRST;
203          has_row = FALSE;
204        while (mtx_next_in_row(N,&c,mtx_ALL_COLS)) {        while (mtx_next_in_row(N,&c,mtx_ALL_COLS)) {
205           printf("x[%d][%d] = %f\n",c.row,c.col,mtx_value(N,&c));           has_row = TRUE;
206             printf("\tmtx[%d][%d] = %f\n",c.row,c.col,mtx_value(N,&c));
207        }        }
208        printf("\n");        if (has_row) printf("\n");
209     }     }
210  }  }
211    
212  static print_vector(v, range)  static print_vector(real64 *v, mtx_range_t range)
 mtx_number_t *v;  
 mtx_range_t range;  
   
213  {  {
214      mtx_index_t i;      int32 i;
215    
216      for (i = range.low; i <= range.high; i++) {      for (i = range.low; i <= range.high; i++) {
217        printf("x[%d] = %f\n",i,v[i]);        printf("\tvec[%d] = %f\n",i,v[i]);
218      }      }
219  }  }
220    
221  static print_solution(rhs, range)  static print_solution(real64 *rhs, mtx_range_t range)
 mtx_number_t *rhs;  
 mtx_range_t range;  
222  {  {
223      mtx_index_t i;      int32 i;
224    
225        printf("Solution:\n");
226      for (i = range.low; i <= range.high; i++) {      for (i = range.low; i <= range.high; i++) {
227        printf("x[%d] = %f\n",i,linsol_var_value(sys,rhs,i));        printf("\ts[%d] = %f\n",i,linsol_var_value(sys,rhs,i));
228      }      }
229  }  }
230    
231  static print_residuals(rhs, range)  static print_residuals(real64 *rhs, mtx_range_t range)
 mtx_number_t *rhs;  
 mtx_range_t range;  
232  {  {
233      mtx_index_t i;      int32 i;
234    
235        printf("Residuals:\n");
236      for (i = range.low; i <= range.high; i++) {      for (i = range.low; i <= range.high; i++) {
237        printf("Residual[%d] = %f\n",i,linsol_eqn_residual(sys,rhs,i));        printf("\tr[%d] = %f\n",i,linsol_eqn_residual(sys,rhs,i));
238      }      }
239  }  }
240    
241  main()  int main()
242  {  {
243     mtx_number_t rhs[10];     real64 rhs[10];
244     mtx_range_t range;     mtx_range_t range;
245       mtx_region_t region = {0,4,0,4};
246    
247    #ifdef __WIN32__
248       Asc_PrintPushVTable(&f_vtable);
249    #endif
250    
251     sys = linsol_create();     sys = linsol_create();
252    
253       printf("\n1st Case - Square\n");
254    
255     printf("Make coefficient matrix N\n");     printf("Make coefficient matrix N\n");
256     make_coefficient_matrix_1();     make_coefficient_matrix_1();
257       print_matrix(N, region);
258     printf("Make right hand side(s) D1 & D2\n");     printf("Make right hand side(s) D1 & D2\n");
259     make_rhs();     make_rhs();
260    
# Line 247  main() Line 262  main()
262           ,linsol_number_of_rhs(sys));           ,linsol_number_of_rhs(sys));
263    
264     printf("Set coefficient matrix\n");     printf("Set coefficient matrix\n");
265     linsol_set_coef_matrix(sys,N);     linsol_set_matrix(sys,N);
266     printf("Reorder equations\n");     printf("Reorder equations\n");
267     linsol_reorder(sys);     linsol_reorder(sys,&region);
268    
269       printf("Invert matrix\n");
270       linsol_invert(sys,&region);
271    
272     printf("Solve with inverse D1\n");     printf("Solve D1\n");
273     linsol_solve(sys,D1);     linsol_solve(sys,D1);
274     printf("Rank = %d\n",linsol_rank(sys));     printf("Rank = %d\n",linsol_rank(sys));
275     range.low = 0;     range.low = 0;
# Line 259  main() Line 277  main()
277     print_residuals(D1,range);     print_residuals(D1,range);
278     print_solution(D1,range);     print_solution(D1,range);
279    
280     printf("Solve with inverse D2\n");     printf("Solve D2\n");
281     linsol_solve(sys,D2);     linsol_solve(sys,D2);
282     printf("Rank = %d\n",linsol_rank(sys));     printf("Rank = %d\n",linsol_rank(sys));
283     print_residuals(D2,range);     print_residuals(D2,range);
# Line 271  main() Line 289  main()
289     /* 2nd case */     /* 2nd case */
290    
291     printf("\n2nd case - overspecified\n");     printf("\n2nd case - overspecified\n");
292       printf("Make coefficient matrix N\n");
293     make_coefficient_matrix_2();     make_coefficient_matrix_2();
294     print_matrix(N);     print_matrix(N, region);
295     printf("Computing NTN\n");     printf("Computing NTN\n");
296     NTN = ATA(N);     NTN = ATA(N);
297     printf("NTN\n\n");     print_matrix(NTN, region);
    print_matrix(NTN);  
298     printf("Compute rhs\n");     printf("Compute rhs\n");
299     ATv(N,D1,rhs);     ATv(N,D1,rhs);
300     range.low = 0;     range.low = 0;
301     range.high = 3;     range.high = 3;
302     print_vector(rhs,range);     print_vector(rhs,range);
303     sys = linsol_create();     sys = linsol_create();
304     linsol_set_coef_matrix(sys,NTN);     linsol_set_matrix(sys,NTN);
305     linsol_add_rhs(sys,rhs);     linsol_add_rhs(sys,rhs,FALSE);
306     linsol_reorder(sys);     printf("Reorder matrix\n");
307       linsol_reorder(sys, &region);
308       printf("Invert matrix\n");
309       linsol_invert(sys, &region);
310     printf("Solve with inverse rhs\n");     printf("Solve with inverse rhs\n");
311     linsol_solve(sys,rhs);     linsol_solve(sys,rhs);
312    
# Line 301  main() Line 322  main()
322     printf("Rank = %d\n",linsol_rank(sys));     printf("Rank = %d\n",linsol_rank(sys));
323     print_residuals(rhs,range);     print_residuals(rhs,range);
324     print_solution(rhs,range);     print_solution(rhs,range);
 }  
   
   
   
   
325    
326    #ifdef __WIN32__
327       Asc_PrintRemoveVTable(f_vtable_name);
328    #endif
329    
330      return 0;
331    }
Legend:
Removed from v.60  
changed lines
  Added in v.61
john.pye@anu.edu.au
 ViewVC Help
Powered by ViewVC 1.1.22