solver/example/example2.c

#include "utilities/ascConfig.h"
#include "solver/mtx.h"
#include "solver/linsol.h"

#ifdef __WIN32__
#include "utilities/ascPrintType.h"
#include "utilities/ascPrint.h"
#define f_vtable_name "example_vtable"
static struct Asc_PrintVTable f_vtable = {f_vtable_name, vfprintf, fflush, NULL};
#endif

static linsol_system_t sys;
static mtx_matrix_t N, NTN;
static real64 D1[10], D2[10];

static make_coefficient_matrix_1()
{
   mtx_coord_t coord;
 
   printf("Create N\n");
   N = mtx_create();
   mtx_set_order(N,10);

   /* eqn 0 */

   coord.row = 0;
   coord.col = 0;
   mtx_set_value(N,&coord,1.);

   /* eqn 1 */

   coord.row = 1;
   coord.col = 0;
   mtx_set_value(N,&coord,0.007);
   coord.col = 1;
   mtx_set_value(N,&coord,0.571);
   coord.col = 2;
   mtx_set_value(N,&coord,0.422);

   /* eqn 2 */

   coord.row = 2;
   coord.col = 1;
   mtx_set_value(N,&coord,0.125);
   coord.col = 2;
   mtx_set_value(N,&coord,0.75);
   coord.col = 3;
   mtx_set_value(N,&coord,0.125);

   /* eqn 3 */

   coord.row = 3;
   coord.col = 2;
   mtx_set_value(N,&coord,0.422);
   coord.col = 3;
   mtx_set_value(N,&coord,0.571);
   coord.col = 4;
   mtx_set_value(N,&coord,0.007);

   /* eqn 4 */

   coord.row = 4;
   coord.col = 4;
   mtx_set_value(N,&coord,1.);
}

static make_coefficient_matrix_2()
{
   mtx_coord_t coord;

   printf("Create N\n");
   N = mtx_create();
   mtx_set_order(N,10);

   /* eqn 0 */

   coord.row = 0;
   coord.col = 0;
   mtx_set_value(N,&coord,1.0);

   /* eqn 1 */

   coord.row = 1;
   coord.col = 0;
   mtx_set_value(N,&coord,0.15);
   coord.col = 1;
   mtx_set_value(N,&coord,0.662);
   coord.col = 2;
   mtx_set_value(N,&coord,0.188);

   /* eqn 2 */

   coord.row = 2;
   coord.col = 1;
   mtx_set_value(N,&coord,0.5);
   coord.col = 2;
   mtx_set_value(N,&coord,0.5);

   /* eqn 3 */

   coord.row = 3;
   coord.col = 1;
   mtx_set_value(N,&coord,0.188);
   coord.col = 2;
   mtx_set_value(N,&coord,0.662);
   coord.col = 3;
   mtx_set_value(N,&coord,0.15);

   /* eqn 4 */

   coord.row = 4;
   coord.col = 3;
   mtx_set_value(N,&coord,1.0);
}

static make_rhs()
{
   printf("\tD1 = { 0.0 }   D2 = { 0.0 }\n"
          "\t     { 1.5 }        { 2.0 }\n"
          "\t     { 3.0 }        { 2.5 }\n"
          "\t     { 4.5 }        { 2.0 }\n"
          "\t     { 6.0 }        { 0.0 }\n");
   D1[0] = 0.0;
   D1[1] = 1.5;
   D1[2] = 3.0;
   D1[3] = 4.5;
   D1[4] = 6.0;

   linsol_add_rhs(sys,D1,FALSE);

   D2[0] = 0.0;
   D2[1] = 2.0;
   D2[2] = 2.5;
   D2[3] = 2.0;
   D2[4] = 0.0;

   linsol_add_rhs(sys,D2,FALSE);
}

mtx_matrix_t ATA(mtx_matrix_t N)
{
   mtx_matrix_t NTN;
   int32 i;
   int32 v1,v2,v3;
   mtx_coord_t c1,c2,c3,c4;

   NTN = mtx_create();
   mtx_set_order(NTN,mtx_order(N));

   for (i = 0; i < mtx_order(N); i++) {
      c1.row = i;
      c1.col = mtx_FIRST;
      c2.row = c1.row;
      while (mtx_next_in_row(N,&c1,mtx_ALL_COLS)) {
         v1 = mtx_value(N,&c1);
         c2.col = mtx_FIRST;
         while (mtx_next_in_row(N,&c2,mtx_ALL_COLS)) {
            if (c2.col >= c1.col) {
               v2 = v1*mtx_value(N,&c2);
               c3.row = c2.col;
               c3.col = c1.col;
               v3 = v2 + mtx_value(NTN,&c3);
               mtx_set_value(NTN,&c3,v3);
               if (c3.row != c3.col) {
                  c4.row = c3.col;
                  c4.col = c3.row;
                  v3 = v2 + mtx_value(NTN,&c4);
                  mtx_set_value(NTN,&c4,v3);
               }
            }
         }
      }
   }
   return(NTN);
}

static ATv(mtx_matrix_t A, real64 v[], real64 r[])
{
   mtx_coord_t c;
   int32 i;

   for (i = 0; i < mtx_order(A); i++) r[i] = 0.0;

   for (i = 0; i < mtx_order(A); i++) {
      c.row = i;
      c.col = mtx_FIRST;
      while (mtx_next_in_row(A,&c,mtx_ALL_COLS)) {
         r[c.col] = r[c.col] + mtx_value(A,&c)*v[i];
      }
   }
}

/* print out non-zero elements of a matrix in text by row */
static print_matrix(mtx_matrix_t N, mtx_region_t region)
{
   mtx_coord_t c;
   int32 i;
   int has_row;

   for (i = region.col.low; i <= region.col.high; i++) {
      c.row = i;
      c.col = mtx_FIRST;
      has_row = FALSE;
      while (mtx_next_in_row(N,&c,mtx_ALL_COLS)) {
         has_row = TRUE;
         printf("\tmtx[%d][%d] = %f\n",c.row,c.col,mtx_value(N,&c));
      }
      if (has_row) printf("\n");
   }
}

static print_vector(real64 *v, mtx_range_t range)
{
    int32 i;

    for (i = range.low; i <= range.high; i++) {
      printf("\tvec[%d] = %f\n",i,v[i]);
    }
}

static print_solution(real64 *rhs, mtx_range_t range)
{
    int32 i;

    printf("Solution:\n");
    for (i = range.low; i <= range.high; i++) {
      printf("\ts[%d] = %f\n",i,linsol_var_value(sys,rhs,i));
    }
}

static print_residuals(real64 *rhs, mtx_range_t range)
{
    int32 i;

    printf("Residuals:\n");
    for (i = range.low; i <= range.high; i++) {
      printf("\tr[%d] = %f\n",i,linsol_eqn_residual(sys,rhs,i));
    }
}

int main()
{
   real64 rhs[10];
   mtx_range_t range;
   mtx_region_t region = {0,4,0,4};

#ifdef __WIN32__
   Asc_PrintPushVTable(&f_vtable);
#endif

   sys = linsol_create();

   printf("\n1st Case - Square\n");

   printf("Make coefficient matrix N\n");
   make_coefficient_matrix_1();
   print_matrix(N, region);
   printf("Make right hand side(s) D1 & D2\n");
   make_rhs();

   printf("Number of right hand sides = %d\n"
         ,linsol_number_of_rhs(sys));

   printf("Set coefficient matrix\n");
   linsol_set_matrix(sys,N);
   printf("Reorder equations\n");
   linsol_reorder(sys,&region);

   printf("Invert matrix\n");
   linsol_invert(sys,&region);

   printf("Solve D1\n");
   linsol_solve(sys,D1);
   printf("Rank = %d\n",linsol_rank(sys));
   range.low = 0;
   range.high = 4;
   print_residuals(D1,range);
   print_solution(D1,range);

   printf("Solve D2\n");
   linsol_solve(sys,D2);
   printf("Rank = %d\n",linsol_rank(sys));
   print_residuals(D2,range);
   print_solution(D2,range);

   mtx_destroy(N);
   linsol_destroy(sys);

   /* 2nd case */

   printf("\n2nd case - overspecified\n");
   printf("Make coefficient matrix N\n");
   make_coefficient_matrix_2();
   print_matrix(N, region);
   printf("Computing NTN\n");
   NTN = ATA(N);
   print_matrix(NTN, region);
   printf("Compute rhs\n");
   ATv(N,D1,rhs);
   range.low = 0;
   range.high = 3;
   print_vector(rhs,range);
   sys = linsol_create();
   linsol_set_matrix(sys,NTN);
   linsol_add_rhs(sys,rhs,FALSE);
   printf("Reorder matrix\n");
   linsol_reorder(sys, &region);
   printf("Invert matrix\n");
   linsol_invert(sys, &region);
   printf("Solve with inverse rhs\n");
   linsol_solve(sys,rhs);

   printf("Rank = %d\n",linsol_rank(sys));
   print_residuals(rhs,range);
   print_solution(rhs,range);

   printf("Solve with inverse rhs\n");
   ATv(N,D2,rhs);
   print_vector(rhs,range);
   linsol_rhs_was_changed(sys,rhs);
   linsol_solve(sys,rhs);
   printf("Rank = %d\n",linsol_rank(sys));
   print_residuals(rhs,range);
   print_solution(rhs,range);

#ifdef __WIN32__
   Asc_PrintRemoveVTable(f_vtable_name);
#endif

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