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 METIS
15 TAUCS_CONFIG GENMMD
16 TAUCS_CONFIG COLAMD
17 TAUCS_CONFIG AMD
18 TAUCS_CONFIG GENERIC_COMPLEX
19 TAUCS_CONFIG MALLOC_STUBS
20 TAUCS_CONFIG MATRIX_GENERATORS
21 TAUCS_CONFIG AD_HOC_TEST
22 TAUCS_CONFIG_END
23
24 */
25
26 #include <stdio.h>
27 #include <stdlib.h>
28 #include "taucs.h"
29
30 #ifdef TAUCS_BLAS_UNDERSCORE
31 #define my_dnrm2 dnrm2_
32 #else
33 #define my_dnrm2 dnrm2
34 #endif
35
36 double my_dnrm2();
37
rnorm(taucs_ccs_matrix * A,double * x,double * y,double * aux)38 int rnorm(taucs_ccs_matrix* A, double* x, double* y, double* aux)
39 {
40 int i;
41 int one = 1;
42 double relerr;
43
44 for(i=0; i<A->n; i++) aux[i] = y[i]-x[i];
45 relerr = my_dnrm2(&(A->n),aux,&one)/my_dnrm2(&(A->n),x,&one);
46 taucs_printf("2-norm relative forward error %.2e \n",relerr);
47
48 if (relerr > 1e-8) {
49 taucs_printf("test failed, error seems too high (but\n");
50 taucs_printf("maybe the matrix is too ill-conditioned)\n");
51 return 1;
52 }
53 return 0;
54 }
55
test_spd_orderings(taucs_ccs_matrix * A,double * x,double * y,double * b,double * z)56 int test_spd_orderings(taucs_ccs_matrix* A,
57 double* x, double* y, double* b, double* z)
58 {
59 int rc;
60 char* metis[] = {"taucs.factor.ordering=metis", NULL};
61 char* genmmd[] = {"taucs.factor.ordering=genmmd", NULL};
62 char* colamd[] = {"taucs.factor.ordering=colamd", NULL};
63 char* amd[] = {"taucs.factor.ordering=amd", NULL};
64 void* opt_arg[] = { NULL };
65
66 rc = taucs_factor_solve(A,NULL,1, y,b,metis,opt_arg);
67 if (rc != TAUCS_SUCCESS) return rc;
68 if (rnorm(A,x,y,z)) return TAUCS_ERROR;
69
70 rc = taucs_factor_solve(A,NULL,1, y,b,genmmd,opt_arg);
71 if (rc != TAUCS_SUCCESS) return rc;
72 if (rnorm(A,x,y,z)) return TAUCS_ERROR;
73
74 rc = taucs_factor_solve(A,NULL,1, y,b,amd,opt_arg);
75 if (rc != TAUCS_SUCCESS) return rc;
76 if (rnorm(A,x,y,z)) return TAUCS_ERROR;
77
78 /* colamd should fail on symmetric matrices */
79 rc = taucs_factor_solve(A,NULL,1, y,b,colamd,opt_arg);
80 if (rc == TAUCS_SUCCESS) return TAUCS_ERROR;
81
82 printf("TESING SYMMETRIC ORDERINGS SUCCEDDED\n");
83
84 return TAUCS_SUCCESS;
85 }
86
test_spd_factorizations(taucs_ccs_matrix * A,double * x,double * y,double * b,double * z)87 int test_spd_factorizations(taucs_ccs_matrix* A,
88 double* x, double* y, double* b, double* z)
89 {
90 int rc;
91 char* mf[] = {"taucs.factor.mf=true", NULL};
92 char* ll[] = {"taucs.factor.ll=true", NULL};
93 char* mfmd[] = {"taucs.factor.mf=true", "taucs.maxdepth=5", NULL};
94 char* llmd[] = {"taucs.factor.ll=true", "taucs.maxdepth=5", NULL};
95 char* ooc[] = {"taucs.ooc=true", "taucs.ooc.basename=/tmp/taucs-test", NULL};
96 void* opt_arg[] = { NULL };
97
98 rc = taucs_factor_solve(A,NULL,1, y,b,ooc,opt_arg);
99 if (rc != TAUCS_SUCCESS) return rc;
100 if (rnorm(A,x,y,z)) return TAUCS_ERROR;
101
102 rc = taucs_factor_solve(A,NULL,1, y,b,mf,opt_arg);
103 if (rc != TAUCS_SUCCESS) return rc;
104 if (rnorm(A,x,y,z)) return TAUCS_ERROR;
105
106 rc = taucs_factor_solve(A,NULL,1, y,b,ll,opt_arg);
107 if (rc != TAUCS_SUCCESS) return rc;
108 if (rnorm(A,x,y,z)) return TAUCS_ERROR;
109
110 /* low depth should fail */
111 rc = taucs_factor_solve(A,NULL,1, y,b,mfmd,opt_arg);
112 if (rc == TAUCS_SUCCESS) return TAUCS_ERROR;
113
114 rc = taucs_factor_solve(A,NULL,1, y,b,llmd,opt_arg);
115 if (rc == TAUCS_SUCCESS) return TAUCS_ERROR;
116
117 printf("TESING SPD FACTORIZATIONS SUCCEDDED\n");
118
119 return TAUCS_SUCCESS;
120 }
121
test_spd_factorsolve(taucs_ccs_matrix * A,double * x,double * y,double * b,double * z)122 int test_spd_factorsolve(taucs_ccs_matrix* A,
123 double* x, double* y, double* b, double* z)
124 {
125 int rc;
126 void* F = NULL;
127 char* factor[] = {"taucs.solve=false", NULL};
128 char* solve [] = {"taucs.factor=false", NULL};
129 void* opt_arg[] = { NULL };
130
131 /* solve without a factorization should fail */
132 rc = taucs_factor_solve(A,NULL,1, y,b,solve,opt_arg);
133 if (rc == TAUCS_SUCCESS) return TAUCS_ERROR;
134
135 /* solve without a factorization should fail */
136 rc = taucs_factor_solve(A,&F,1, y,b,solve,opt_arg);
137 if (rc == TAUCS_SUCCESS) return TAUCS_ERROR;
138
139 rc = taucs_factor_solve(A,&F,1, y,b,factor,opt_arg);
140 if (rc != TAUCS_SUCCESS) return rc;
141
142 rc = taucs_factor_solve(A,&F,1, y,b,solve,opt_arg);
143 if (rc != TAUCS_SUCCESS) return rc;
144 if (rnorm(A,x,y,z)) return TAUCS_ERROR;
145
146 printf("TESING SPD FACTORSOLVE SUCCEDDED\n");
147
148 return TAUCS_SUCCESS;
149 }
150
main()151 int main()
152 {
153 int xyz = 20;
154 int i;
155
156 taucs_ccs_matrix* A;
157 void* f;
158 void* fc;
159
160 double* X;
161 double* B;
162 double* Y;
163 double* Z;
164
165 taucs_logfile("stdout");
166
167 A = taucs_ccs_generate_mesh3d(xyz,xyz,xyz);
168 if (!A) {
169 taucs_printf("Matrix generation failed\n");
170 return 1;
171 }
172
173 X =(double*)malloc((A->n)*sizeof(double));
174 B =(double*)malloc((A->n)*sizeof(double));
175 Y =(double*)malloc((A->n)*sizeof(double));
176 Z =(double*)malloc((A->n)*sizeof(double));
177 if (!X || !B || !Y || !Z) {
178 taucs_printf("Vector allocation failed\n");
179 return 1;
180 }
181
182 for(i=0; i<A->n; i++) X[i]=((double)random()/RAND_MAX);
183 taucs_ccs_times_vec(A,X,B);
184
185 if (test_spd_factorsolve(A,X,Y,B,Z)) {
186 printf("SPD FACTOR-SOLVE FAILED\n");
187 return 1;
188 }
189
190 if (test_spd_orderings(A,X,Y,B,Z)) {
191 printf("SPD ORDERING FAILED\n");
192 return 1;
193 }
194
195 if (test_spd_factorizations(A,X,Y,B,Z)) {
196 printf("SPD FACTORIZATION FAILED\n");
197 return 1;
198 }
199
200
201 taucs_printf("test succeeded\n");
202 return 0;
203 }
204