1 /*********************************************************/
2 /* TAUCS */
3 /* Author: Sivan Toledo */
4 /*********************************************************/
5 /*
6
7 TAUCS_CONFIG_BEGIN
8 TAUCS_CONFIG_DEFAULT OFF
9 TAUCS_CONFIG BASE
10 xxxTAUCS_CONFIG CILK
11 TAUCS_CONFIG DREAL
12 TAUCS_CONFIG FACTOR
13 TAUCS_CONFIG OOC_LLT
14 TAUCS_CONFIG INCOMPLETE_CHOL
15 TAUCS_CONFIG VAIDYA
16 TAUCS_CONFIG METIS
17 TAUCS_CONFIG GENMMD
18 TAUCS_CONFIG COLAMD
19 TAUCS_CONFIG AMD
20 TAUCS_CONFIG MALLOC_STUBS
21 TAUCS_CONFIG MATRIX_GENERATORS
22 TAUCS_CONFIG AD_HOC_TEST
23 TAUCS_CONFIG_END
24
25 */
26
27 #include <stdio.h>
28 #include <stdlib.h>
29 #include "taucs.h"
30
rnorm(taucs_ccs_matrix * A,void * x,void * b,void * aux)31 int rnorm(taucs_ccs_matrix* A, void* x, void* b, void* aux)
32 {
33 double relerr;
34
35 taucs_ccs_times_vec(A,x,aux);
36
37 taucs_vec_axpby(A->n,A->flags,1.0,aux,-1.0,b,aux);
38
39 relerr = taucs_vec_norm2(A->n,A->flags,aux)
40 / taucs_vec_norm2(A->n,A->flags,b);
41
42
43 taucs_printf("relative 2-norm of the residual %.2e \n",relerr);
44
45 if (relerr > 1e-8) {
46 taucs_printf("test failed, error seems too high (but\n");
47 taucs_printf("maybe the matrix is too ill-conditioned)\n");
48 return 1;
49 }
50 return 0;
51 }
52
test_spd_orderings(taucs_ccs_matrix * A,double * x,double * y,double * b,double * z)53 int test_spd_orderings(taucs_ccs_matrix* A,
54 double* x, double* y, double* b, double* z)
55 {
56 int rc;
57 char* metis[] = {"taucs.factor.LLT=true", "taucs.factor.ordering=metis", NULL};
58 char* genmmd[] = {"taucs.factor.LLT=true", "taucs.factor.ordering=genmmd", NULL};
59 char* amd[] = {"taucs.factor.LLT=true", "taucs.factor.ordering=amd", NULL};
60 char* colamd[] = {"taucs.factor.LLT=true", "taucs.factor.ordering=colamd", NULL};
61 void* opt_arg[] = { NULL };
62 int test = 200;
63
64 printf("TEST %d\n",test++);
65 rc = taucs_linsolve(A,NULL,1, y,b,metis,opt_arg);
66 if (rc != TAUCS_SUCCESS) return rc;
67 if (rnorm(A,y,b,z)) return TAUCS_ERROR;
68
69 printf("TEST %d\n",test++);
70 rc = taucs_linsolve(A,NULL,1, y,b,genmmd,opt_arg);
71 if (rc != TAUCS_SUCCESS) return rc;
72 if (rnorm(A,y,b,z)) return TAUCS_ERROR;
73
74 printf("TEST %d\n",test++);
75 rc = taucs_linsolve(A,NULL,1, y,b,amd,opt_arg);
76 if (rc != TAUCS_SUCCESS) return rc;
77 if (rnorm(A,y,b,z)) return TAUCS_ERROR;
78
79 /* colamd should fail on symmetric matrices */
80 printf("TEST %d\n",test++);
81 rc = taucs_linsolve(A,NULL,1, y,b,colamd,opt_arg);
82 if (rc == TAUCS_SUCCESS) return TAUCS_ERROR;
83
84 printf("TESING SYMMETRIC ORDERINGS SUCCEEDED\n");
85
86 return TAUCS_SUCCESS;
87 }
88
test_spd_factorizations(taucs_ccs_matrix * A,double * x,double * y,double * b,double * z)89 int test_spd_factorizations(taucs_ccs_matrix* A,
90 double* x, double* y, double* b, double* z)
91 {
92 int rc;
93 char* mf[] = {"taucs.factor.LLT=true", "taucs.factor.mf=true", NULL};
94 char* ll[] = {"taucs.factor.LLT=true", "taucs.factor.ll=true", NULL};
95 char* mfmd[] = {"taucs.factor.LLT=true", "taucs.factor.mf=true", "taucs.maxdepth=5", NULL};
96 char* llmd[] = {"taucs.factor.LLT=true", "taucs.factor.ll=true", "taucs.maxdepth=5", NULL};
97 char* ooc[] = {"taucs.factor.LLT=true", "taucs.ooc=true", "taucs.ooc.basename=taucs-test", NULL};
98 void* opt_arg[] = { NULL };
99
100 rc = taucs_linsolve(A,NULL,1, y,b,ooc,opt_arg);
101 if (rc != TAUCS_SUCCESS) return rc;
102 if (rnorm(A,y,b,z)) return TAUCS_ERROR;
103
104 rc = taucs_linsolve(A,NULL,1, y,b,mf,opt_arg);
105 if (rc != TAUCS_SUCCESS) return rc;
106 if (rnorm(A,y,b,z)) return TAUCS_ERROR;
107
108 rc = taucs_linsolve(A,NULL,1, y,b,ll,opt_arg);
109 if (rc != TAUCS_SUCCESS) return rc;
110 if (rnorm(A,y,b,z)) return TAUCS_ERROR;
111
112 /* low depth should fail */
113 rc = taucs_linsolve(A,NULL,1, y,b,mfmd,opt_arg);
114 if (rc == TAUCS_SUCCESS) return TAUCS_ERROR;
115
116 rc = taucs_linsolve(A,NULL,1, y,b,llmd,opt_arg);
117 if (rc == TAUCS_SUCCESS) return TAUCS_ERROR;
118
119 printf("TESING SPD FACTORIZATIONS SUCCEEDED\n");
120
121 return TAUCS_SUCCESS;
122 }
123
test_spd_factorsolve(taucs_ccs_matrix * A,double * x,double * y,double * b,double * z)124 int test_spd_factorsolve(taucs_ccs_matrix* A,
125 double* x, double* y, double* b, double* z)
126 {
127 int rc;
128 void* F = NULL;
129 char* factor[] = {"taucs.factor.LLT=true", NULL};
130 char* solve [] = {"taucs.factor=false", NULL};
131 void* opt_arg[] = { NULL };
132 int test = 100;
133
134 printf("TESING SPD FACTORSOLVE\n");
135
136 /* solve without a factorization should fail */
137 printf("TEST %d\n",test++);
138 rc = taucs_linsolve(A,NULL,1, y,b,solve,opt_arg);
139 if (rc == TAUCS_SUCCESS) return TAUCS_ERROR;
140
141 /* solve without a factorization should fail */
142 printf("TEST %d\n",test++);
143 rc = taucs_linsolve(A,&F,1, y,b,solve,opt_arg);
144 if (rc == TAUCS_SUCCESS) return TAUCS_ERROR;
145
146 /* this should work, factor, solve, free */
147 printf("TEST %d\n",test++);
148 rc = taucs_linsolve(A,&F,1, y,b,factor,opt_arg);
149 if (rc != TAUCS_SUCCESS) return rc;
150 if (rnorm(A,y,b,z)) return TAUCS_ERROR;
151 printf("TEST %d\n",test++);
152 rc = taucs_linsolve(NULL,&F,0, NULL,NULL,factor,opt_arg);
153 if (rc != TAUCS_SUCCESS) return rc;
154
155 /* now first factor, then solve, then free */
156 printf("TEST %d\n",test++);
157 rc = taucs_linsolve(A,&F,0,NULL,NULL,factor,opt_arg);
158 if (rc != TAUCS_SUCCESS) return rc;
159 printf("TEST %d\n",test++);
160 rc = taucs_linsolve(A,&F,1,y,b,solve,opt_arg);
161 if (rc != TAUCS_SUCCESS) return rc;
162 printf("TEST %d\n",test++);
163 rc = taucs_linsolve(NULL,&F,0, NULL,NULL,factor,opt_arg);
164 if (rc != TAUCS_SUCCESS) return rc;
165 if (rnorm(A,y,b,z)) return TAUCS_ERROR;
166
167 /* factor+solve for 4 rhs */
168 printf("TEST %d\n",test++);
169 rc = taucs_linsolve(A,&F,4, y,b,factor,opt_arg);
170 if (rc != TAUCS_SUCCESS) return rc;
171 if (rnorm(A,y,b,z)) return TAUCS_ERROR;
172 if (rnorm(A,y+1*(A->n),b+1*(A->n),z+1*(A->n))) return TAUCS_ERROR;
173 if (rnorm(A,y+2*(A->n),b+2*(A->n),z+2*(A->n))) return TAUCS_ERROR;
174 if (rnorm(A,y+3*(A->n),b+3*(A->n),z+3*(A->n))) return TAUCS_ERROR;
175 printf("TEST %d\n",test++);
176 rc = taucs_linsolve(NULL,&F,0, NULL,NULL,factor,opt_arg);
177 if (rc != TAUCS_SUCCESS) return rc;
178
179 printf("TESING SPD FACTORSOLVE SUCCEEDED\n");
180
181 return TAUCS_SUCCESS;
182 }
183
main()184 int main()
185 {
186 int n;
187 int xyz = 20;
188 int i;
189
190 taucs_ccs_matrix* A;
191
192 double* X;
193 double* B;
194 double* Y;
195 double* Z;
196
197 taucs_logfile("stdout");
198
199 A = taucs_ccs_generate_mesh3d(xyz,xyz,xyz);
200 if (!A) {
201 taucs_printf("Matrix generation failed\n");
202 return 1;
203 }
204
205 n = A->n;
206
207 X =(double*)malloc(4*n*sizeof(double));
208 B =(double*)malloc(4*n*sizeof(double));
209 Y =(double*)malloc(4*n*sizeof(double));
210 Z =(double*)malloc(4*n*sizeof(double));
211 if (!X || !B || !Y || !Z) {
212 taucs_printf("Vector allocation failed\n");
213 return 1;
214 }
215
216 for(i=0; i<4*n; i++) X[i]=((double)rand()/RAND_MAX);
217 taucs_ccs_times_vec(A,X ,B);
218 taucs_ccs_times_vec(A,X+1*n,B+1*n);
219 taucs_ccs_times_vec(A,X+2*n,B+2*n);
220 taucs_ccs_times_vec(A,X+3*n,B+3*n);
221
222 if (test_spd_factorsolve(A,X,Y,B,Z)) {
223 printf("SPD FACTOR-SOLVE FAILED\n");
224 return 1;
225 }
226
227 if (test_spd_orderings(A,X,Y,B,Z)) {
228 printf("SPD ORDERING FAILED\n");
229 return 1;
230 }
231
232 if (test_spd_factorizations(A,X,Y,B,Z)) {
233 printf("SPD FACTORIZATION FAILED\n");
234 return 1;
235 }
236
237 taucs_printf("test succeeded\n");
238 return 0;
239 }
240