solver/example/example2.c

#include "base.h"
#include "linsol.h"

static linsol_system_t sys;
static mtx_matrix_t N,NTN;
static mtx_number_t 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()
{
   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);

   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);
}

mtx_matrix_t ATA(N)
mtx_matrix_t N;
{
    mtx_matrix_t NTN;
    mtx_index_t i;
    mtx_number_t 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(A,v,r)
mtx_matrix_t A;
mtx_number_t v[],r[];
{
   mtx_coord_t c;
   mtx_index_t 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];
      }
   }
}

static print_matrix(N)
mtx_matrix_t N;
{
   mtx_coord_t c;
   mtx_index_t i;

   for (i = 0; i < mtx_order(N); i++) {
      c.row = i;
      c.col = mtx_FIRST;
      while (mtx_next_in_row(N,&c,mtx_ALL_COLS)) {
         printf("x[%d][%d] = %f\n",c.row,c.col,mtx_value(N,&c));
      }
      printf("\n");
   }
}

static print_vector(v, range)
mtx_number_t *v;
mtx_range_t range;

{
    mtx_index_t i;

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

static print_solution(rhs, range)
mtx_number_t *rhs;
mtx_range_t range;
{
    mtx_index_t i;

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

static print_residuals(rhs, range)
mtx_number_t *rhs;
mtx_range_t range;
{
    mtx_index_t i;

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

main()
{
   mtx_number_t rhs[10];
   mtx_range_t range;

   sys = linsol_create();

   printf("Make coefficient matrix N\n");
   make_coefficient_matrix_1();
   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_coef_matrix(sys,N);
   printf("Reorder equations\n");
   linsol_reorder(sys);

   printf("Solve with inverse 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 with inverse 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");
   make_coefficient_matrix_2();
   print_matrix(N);
   printf("Computing NTN\n");
   NTN = ATA(N);
   printf("NTN\n\n");
   print_matrix(NTN);
   printf("Compute rhs\n");
   ATv(N,D1,rhs);
   range.low = 0;
   range.high = 3;
   print_vector(rhs,range);
   sys = linsol_create();
   linsol_set_coef_matrix(sys,NTN);
   linsol_add_rhs(sys,rhs);
   linsol_reorder(sys);
   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);
}


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