1 /*********************************************************/
2 /* TAUCS */
3 /* Author: Sivan Toledo */
4 /*********************************************************/
5
6 #include <stdio.h>
7 #include <stdlib.h>
8 #include <assert.h>
9 #include <math.h>
10 #include <string.h>
11 #include "taucs.h"
12
13 #define RNDM ((double)rand()/(double)RAND_MAX);
14
15 /*ifndef added omer*/
16 #ifndef max
17 #define max(x,y) ( ((x) > (y)) ? (x) : (y) )
18 #endif
19 #ifndef mod
20 #define mod(x,n) ((x) % (n))
21 #endif
22 /*********************************************************/
23 /* */
24 /*********************************************************/
25
26 #ifdef TAUCS_CORE_DOUBLE
27
taucs_ccs_generate_mesh2d_negative(int n)28 taucs_ccs_matrix* taucs_ccs_generate_mesh2d_negative(int n)
29 {
30 taucs_ccs_matrix* m;
31 int N;
32 int nnz;
33 int x,y,i,j,ip;
34
35 taucs_printf("generate_mesh2d_negative: starting\n");
36
37 m = (taucs_ccs_matrix*) taucs_malloc(sizeof(taucs_ccs_matrix));
38 if (!m) {
39 taucs_printf("generate_mesh2d_negative: out of memory (1)\n");
40 return NULL;
41 }
42
43 N = n*n;
44 nnz = 4*N;
45
46 m->n = N;
47 m->flags = TAUCS_SYMMETRIC | TAUCS_LOWER | TAUCS_DOUBLE;
48 m->colptr = (int*) taucs_malloc((N+1) * sizeof(int));
49 m->rowind = (int*) taucs_malloc(nnz * sizeof(int));
50 m->values.d/*taucs_values*/ = (double*) taucs_malloc(nnz * sizeof(double));
51
52 if (!(m->colptr) || !(m->rowind) || !(m->rowind)) {
53 taucs_printf("generate_mesh2d_negative: out of memory (4): ncols=%d nnz=%d\n",N,nnz);
54 taucs_free(m->colptr); taucs_free(m->rowind); taucs_free(m->values.d/*taucs_values*/);
55 return NULL;
56 }
57
58 ip = 0;
59 for (y=0; y<n; y++) {
60 for (x=0; x<n; x++) {
61 j = x + y*n; /* convert mesh (x,y) location to index in vector */
62 /*printf("column %d xy %d,%d starts at %d\n",j,x,y,ip);*/
63 (m->colptr)[j] = ip;
64
65 i=mod(x+1,n) + (y )*n ; if (i>j) { (m->rowind)[ip]=i; (m->values.d/*taucs_values*/)[ip]=-1.0; ip++; }
66 i=(x ) + mod(y+1,n)*n ; if (i>j) { (m->rowind)[ip]=i; (m->values.d/*taucs_values*/)[ip]=+100.0; ip++; }
67 i=mod(x+n-1,n) + (y )*n ; if (i>j) { (m->rowind)[ip]=i; (m->values.d/*taucs_values*/)[ip]=-1.0; ip++; }
68 i=(x ) + mod(y+n-1,n)*n ; if (i>j) { (m->rowind)[ip]=i; (m->values.d/*taucs_values*/)[ip]=+100.0; ip++; }
69
70 /* i=mod(x+1,n) + mod(y+1,n)*n ; (m->rowind)[ip]=i; (m->taucs_values)[ip]=+1.0; ip++; */
71 /* i=mod(x+2,n) + (y )*n ; (m->rowind)[ip]=i; (m->taucs_values)[ip]=.0625; ip++; */
72 /* i=(x ) + mod(y+2,n)*n ; (m->rowind)[ip]=i; (m->taucs_values)[ip]=.0625; ip++; */
73
74 i=(x )+(y )*n; (m->rowind)[ip]=i;
75 /* (m->taucs_values)[ip]= 4.25; if (x==0 && y==0) (m->taucs_values)[ip] += 1; to make it nonsingular */
76 (m->values.d/*taucs_values*/)[ip]= 202.0; if (x==0 && y==0) (m->values.d/*taucs_values*/)[ip] += 1; /* to make it nonsingular */
77 ip++;
78 }
79 }
80
81 (m->colptr)[N] = ip;
82
83 taucs_printf("generate_mesh2d_negative: done: ncols=%d nnz=%d\n",N,ip);
84
85 return m;
86 }
87
88
89 taucs_ccs_matrix*
taucs_ccs_generate_mesh2d(int n,char * which)90 taucs_ccs_generate_mesh2d(int n,char *which)
91 {
92 taucs_ccs_matrix* m;
93 int N;
94 int nnz;
95 int x,y,i,j,ip;
96 double jump = 100;
97
98 taucs_printf("taucs_ccs_generate_mesh2d: starting\n");
99
100 m = (taucs_ccs_matrix*) taucs_malloc(sizeof(taucs_ccs_matrix));
101 if (!m) {
102 taucs_printf("generate_mesh2d: out of memory (1)\n");
103 return NULL;
104 }
105
106 N = n*n;
107 nnz = 3*N;
108
109 m->n = N;
110 m->m = N;
111 m->flags = TAUCS_SYMMETRIC | TAUCS_LOWER | TAUCS_DOUBLE;
112 m->colptr = (int*) taucs_malloc((N+1) * sizeof(int));
113 m->rowind = (int*) taucs_malloc(nnz * sizeof(int));
114 m->values.d/*taucs_values*/ = (double*) taucs_malloc(nnz * sizeof(double));
115
116 if (!(m->colptr) || !(m->rowind) || !(m->rowind)) {
117 taucs_printf("taucs_ccs_generate_mesh2d: out of memory: ncols=%d nnz=%d\n",N,nnz);
118 taucs_free(m->colptr); taucs_free(m->rowind); taucs_free(m->values.d/*taucs_values*/);
119 return NULL;
120 }
121
122 ip = 0;
123 for (y=0; y<n; y++) {
124 for (x=0; x<n; x++) {
125 j = x + y*n; /* convert mesh (x,y) location to index in vector */
126 /*printf("column %d xy %d,%d starts at %d\n",j,x,y,ip);*/
127 (m->colptr)[j] = ip;
128 /* if (x < n-1) { i=(x+1)+(y )*n; (m->rowind)[ip]=i; (m->taucs_values)[ip]=-1.0; ip++; } */
129
130 if (!strcmp(which,"anisotropic_y")) {
131 if (y < n-1) { i=(x )+(y+1)*n; (m->rowind)[ip]=i; (m->values.d/*taucs_values*/)[ip]=-jump; ip++; }
132 } else
133 if (y < n-1) { i=(x )+(y+1)*n; (m->rowind)[ip]=i; (m->values.d/*taucs_values*/)[ip]=-1.0; ip++; }
134
135 if (!strcmp(which,"anisotropic_x")) {
136 if (x < n-1) { i=(x+1)+(y )*n; (m->rowind)[ip]=i; (m->values.d/*taucs_values*/)[ip]=-jump; ip++; }
137 } else
138 if (x < n-1) { i=(x+1)+(y )*n; (m->rowind)[ip]=i; (m->values.d/*taucs_values*/)[ip]=-1.0; ip++; }
139
140 if (!strcmp(which,"anisotropic_y"))
141 {
142 i=(x )+(y )*n; (m->rowind)[ip]=i;
143 (m->values.d/*taucs_values*/)[ip]= 0.0;
144 if (x > 0) (m->values.d/*taucs_values*/)[ip] += 1.0;
145 if (y > 0) (m->values.d/*taucs_values*/)[ip] += jump;
146 if (x < n-1) (m->values.d/*taucs_values*/)[ip] += 1.0;
147 if (y < n-1) (m->values.d/*taucs_values*/)[ip] += jump;
148 if (x==0 && y==0) (m->values.d/*taucs_values*/)[ip] += 1.0; /* to make it nonsingular */
149 ip++;
150 }
151 else if (!strcmp(which,"anisotropic_x"))
152 {
153 i=(x )+(y )*n; (m->rowind)[ip]=i;
154 (m->values.d/*taucs_values*/)[ip]= 0.0;
155 if (x > 0) (m->values.d/*taucs_values*/)[ip] += jump;
156 if (y > 0) (m->values.d/*taucs_values*/)[ip] += 1.0;
157 if (x < n-1) (m->values.d/*taucs_values*/)[ip] += jump;
158 if (y < n-1) (m->values.d/*taucs_values*/)[ip] += 1.0;
159 if (x==0 && y==0) (m->values.d/*taucs_values*/)[ip] += 1.0; /* to make it nonsingular */
160 ip++;
161 }
162 else if (!strcmp(which,"dirichlet"))
163 {
164 i=(x )+(y )*n; (m->rowind)[ip]=i; (m->values.d/*taucs_values*/)[ip]= 4.0;
165 ip++;
166 }
167 else /* neumann */
168 {
169 i=(x )+(y )*n; (m->rowind)[ip]=i;
170 (m->values.d/*taucs_values*/)[ip]= 0.0;
171 if (x > 0) (m->values.d/*taucs_values*/)[ip] += 1.0;
172 if (y > 0) (m->values.d/*taucs_values*/)[ip] += 1.0;
173 if (x < n-1) (m->values.d/*taucs_values*/)[ip] += 1.0;
174 if (y < n-1) (m->values.d/*taucs_values*/)[ip] += 1.0;
175 if (x==0 && y==0) (m->values.d/*taucs_values*/)[ip] += 1.0; /* to make it nonsingular */
176 ip++;
177 }
178
179 }
180 }
181 (m->colptr)[N] = ip;
182
183 taucs_printf("taucs_ccs_generate_mesh2d: done, ncols=%d nnz=%d\n",N,ip);
184
185 /*
186 for (j=0; j<N; j++) {
187 for (ip=(m->colptr)[j]; ip < (m->colptr)[j+1]; ip++) {
188 i = (m->rowind)[ip];
189 taucs_printf("<%d %d %lg>\n",i,j,m->taucs_values[ip]);
190 }
191 }
192 */
193
194 return m;
195 }
196
197 taucs_ccs_matrix*
taucs_ccs_generate_mesh3d(int X,int Y,int Z)198 taucs_ccs_generate_mesh3d(int X, int Y, int Z)
199 {
200 taucs_ccs_matrix* m;
201 int N;
202 int nnz;
203 int x,y,z,i,j,ip;
204
205 taucs_printf("taucs_ccs_generate_mesh3d: starting\n");
206
207 m = (taucs_ccs_matrix*) taucs_malloc(sizeof(taucs_ccs_matrix));
208 if (!m) {
209 taucs_printf("taucs_ccs_generate_mesh3d: out of memory\n");
210 return NULL;
211 }
212
213 N = X*Y*Z;
214 nnz = 4*N;
215
216 m->n = N;
217 m->m = N;
218 m->flags = TAUCS_SYMMETRIC | TAUCS_LOWER | TAUCS_DOUBLE;
219 /*m->indshift = 0;*/
220 m->colptr = (int*) taucs_malloc((N+1) * sizeof(int));
221 m->rowind = (int*) taucs_malloc(nnz * sizeof(int));
222 m->values.d/*taucs_values*/ = (double*) taucs_malloc(nnz * sizeof(double));
223
224 if (!(m->colptr) || !(m->rowind) || !(m->rowind)) {
225 taucs_printf("taucs_ccs_generate_mesh3d: out of memory: ncols=%d nnz=%d\n",N,nnz);
226 taucs_free(m->colptr); taucs_free(m->rowind); taucs_free(m->values.d/*taucs_values*/);
227 return NULL;
228 }
229
230 ip = 0;
231 for (z=0; z<Z; z++) {
232 for (y=0; y<Y; y++) {
233 for (x=0; x<X; x++) {
234 j = z*X*Y + y*X + x;
235 /*printf("column %d xy %d,%d starts at %d\n",j,x,y,ip);*/
236 (m->colptr)[j] = ip;
237 if (x < X-1) { i=(z )*X*Y+(y )*X+(x+1); (m->rowind)[ip]=i; (m->values.d/*taucs_values*/)[ip]=-1.0; ip++; }
238 if (y < Y-1) { i=(z )*X*Y+(y+1)*X+(x ); (m->rowind)[ip]=i; (m->values.d/*taucs_values*/)[ip]=-1.0; ip++; }
239 if (z < Z-1) { i=(z+1)*X*Y+(y )*X+(x ); (m->rowind)[ip]=i; (m->values.d/*taucs_values*/)[ip]=-1.0; ip++; }
240 {
241 i=(z )*X*Y+(y )*X+(x ); (m->rowind)[ip]=i;
242 (m->values.d/*taucs_values*/)[ip]= 0.0;
243 if (x < X-1) (m->values.d/*taucs_values*/)[ip] += 1.0;
244 if (y < Y-1) (m->values.d/*taucs_values*/)[ip] += 1.0;
245 if (z < Z-1) (m->values.d/*taucs_values*/)[ip] += 1.0;
246 if (x > 0 ) (m->values.d/*taucs_values*/)[ip] += 1.0;
247 if (y > 0 ) (m->values.d/*taucs_values*/)[ip] += 1.0;
248 if (z > 0 ) (m->values.d/*taucs_values*/)[ip] += 1.0;
249 if (x==0 && y==0 && z==0) (m->values.d/*taucs_values*/)[ip] += 1.0;
250 ip++;
251 }
252 /* { i=(z )*X*Y+(y )*X+(x ); (m->rowind)[ip]=i; (m->taucs_values)[ip]= 6.0; ip++; } */
253 }
254 }
255 }
256 (m->colptr)[N] = ip;
257
258 taucs_printf("taucs_ccs_generate_mesh3d: done, ncols=%d nnz=%d\n",N,ip);
259
260 return m;
261 }
262
263 taucs_ccs_matrix*
taucs_ccs_generate_dense(int M,int N,int flags)264 taucs_ccs_generate_dense(int M, int N, int flags)
265 {
266 taucs_ccs_matrix* m;
267 int nnz;
268 int i,j,ip;/* x,y omer*/
269
270 taucs_printf("taucs_ccs_generate_dense: starting\n");
271
272 m = (taucs_ccs_matrix*) taucs_malloc(sizeof(taucs_ccs_matrix));
273 if (!m) {
274 taucs_printf("taucs_ccs_generate_dense: out of memory\n");
275 return NULL;
276 }
277
278 m->m = N;
279 m->n = N;
280 if (flags & TAUCS_SYMMETRIC) {
281 nnz = N*(N+1)/2;
282 m->flags = TAUCS_SYMMETRIC | TAUCS_LOWER | TAUCS_DOUBLE;
283 } else {
284 nnz = N*N;
285 m->flags = TAUCS_DOUBLE;
286 }
287
288 m->colptr = (int*) taucs_malloc((N+1) * sizeof(int));
289 m->rowind = (int*) taucs_malloc(nnz * sizeof(int));
290 m->values.d/*taucs_values*/ = (double*) taucs_malloc(nnz * sizeof(double));
291
292 if (!(m->colptr) || !(m->rowind) || !(m->rowind)) {
293 taucs_printf("taucs_ccs_generate_dense: out of memory: nrows=%d ncols=%d nnz=%d\n",M,N,nnz);
294 taucs_free(m->colptr); taucs_free(m->rowind); taucs_free(m->values.d/*taucs_values*/);
295 return NULL;
296 }
297
298 ip = 0;
299 for (j=0; j<N; j++) {
300 (m->colptr)[j] = ip;
301 if (flags & TAUCS_SYMMETRIC) {
302 for (i=j; i<N; i++) {
303 (m->rowind)[ip]=i;
304 (m->values.d/*taucs_values*/)[ip]=RNDM;
305 ip++;
306 }
307 } else {
308 for (i=0; i<M; i++) {
309 (m->rowind)[ip]=i;
310 (m->values.d/*taucs_values*/)[ip]=RNDM;
311 ip++;
312 }
313 }
314 }
315 (m->colptr)[N] = ip;
316
317 taucs_printf("taucs_ccs_generate_dense: done, nrows=%d ncols=%d nnz=%d\n",M,N,ip);
318
319 return m;
320 }
321
322 /* random resistor networks */
323
recursive_visit(int i,int * neighbors[],int degree[],int visited[])324 int recursive_visit(int i,
325 int* neighbors[],
326 int degree[],
327 int visited[])
328 {
329 int j,jp,count;
330 visited[i] = 1;
331 count = 1;
332 for (jp=0; jp<degree[i]; jp++) {
333 j = neighbors[i][jp];
334 if (! visited[j] ) count += recursive_visit(j,neighbors,degree,visited);
335 }
336 return count;
337 }
338
339 taucs_ccs_matrix*
taucs_ccs_generate_rrn(int X,int Y,int Z,double drop_probability,double rmin)340 taucs_ccs_generate_rrn(int X, int Y, int Z, double drop_probability, double rmin)
341 {
342 taucs_ccs_matrix* m;
343 taucs_ccs_matrix* l;
344 int N;
345 int nnz;
346 int x,y,z,i,j,k,ip,jp;
347 double* D; /* contributions to future diagonal elements */
348
349 int** neighbors;
350 int* degree;
351 int* visited;
352 int* reps;
353 int ncomponents;
354
355 int largest;
356 int largest_rep = 0; /* warning */
357
358 taucs_printf("taucs_ccs_generate_rrn: starting (%d %d %d %.4e %.4e)\n",
359 X,Y,Z,drop_probability,rmin);
360
361 if (drop_probability > 1.0 || drop_probability < 0.0) {
362 taucs_printf("taucs_ccs_generate_rrn: drop probability (%lg) must be in [0,1], setting to 0\n",
363 drop_probability);
364 drop_probability = 0.0;
365 }
366
367 if (rmin > 1.0 || rmin <= 0.0) {
368 taucs_printf("taucs_ccs_generate_rrn: rmin (%lg) must be in (0,1], setting to 1\n",
369 rmin);
370 rmin = 1.0;
371 }
372
373 m = (taucs_ccs_matrix*) taucs_malloc(sizeof(taucs_ccs_matrix));
374 if (!m) {
375 taucs_printf("taucs_ccs_generate_rrn: out of memory\n");
376 return NULL;
377 }
378
379 N = X*Y*Z;
380 nnz = 4*N; /* this is an upper bound */
381
382 m->n = N;
383 m->m = N;
384 m->flags = TAUCS_SYMMETRIC | TAUCS_LOWER | TAUCS_DOUBLE;
385 /*m->indshift = 0;*/
386 m->colptr = (int*) taucs_malloc((N+1) * sizeof(int));
387 m->rowind = (int*) taucs_malloc(nnz * sizeof(int));
388 m->values.d/*taucs_values*/ = (double*) taucs_malloc(nnz * sizeof(double));
389
390 D = (double*) taucs_malloc(N * sizeof(double));
391
392 if (!(m->colptr) || !(m->rowind) || !(m->rowind) || !D) {
393 taucs_printf("taucs_ccs_generate_rrn: out of memory: ncols=%d nnz=%d\n",N,nnz);
394 taucs_free(m->colptr); taucs_free(m->rowind); taucs_free(m->values.d/*taucs_values*/); taucs_free(D);
395 return NULL;
396 }
397
398 for (i=0; i<N; i++) D[i] = 0.0;
399
400 ip = 0;
401 for (z=0; z<Z; z++) {
402 for (y=0; y<Y; y++) {
403 for (x=0; x<X; x++) {
404 int j, je, jw, js, jn, ju, jd; /* indices for up, down, east, west, south, north */
405 int jp; /* pointer to the diagonal value */
406 double v;
407
408 j = z*X*Y + y*X + x;
409 jw = (x > 0 ) ? (z )*X*Y + (y )*X + (x-1) : (z )*X*Y + (y )*X + (X-1) ;
410 je = (x < X-1) ? (z )*X*Y + (y )*X + (x+1) : (z )*X*Y + (y )*X + (0 ) ;
411 js = (y > 0 ) ? (z )*X*Y + (y-1)*X + (x ) : (z )*X*Y + (Y-1)*X + (x ) ;
412 jn = (y < Y-1) ? (z )*X*Y + (y+1)*X + (x ) : (z )*X*Y + (0 )*X + (x ) ;
413 jd = (z > 0 ) ? (z-1)*X*Y + (y )*X + (x ) : (Z-1)*X*Y + (y )*X + (x ) ;
414 ju = (z < Z-1) ? (z+1)*X*Y + (y )*X + (x ) : (0 )*X*Y + (y )*X + (x ) ;
415
416 jw = (x > 0 ) ? (z )*X*Y + (y )*X + (x-1) : -1;
417 je = (x < X-1) ? (z )*X*Y + (y )*X + (x+1) : -1;
418 js = (y > 0 ) ? (z )*X*Y + (y-1)*X + (x ) : -1;
419 jn = (y < Y-1) ? (z )*X*Y + (y+1)*X + (x ) : -1;
420 jd = (z > 0 ) ? (z-1)*X*Y + (y )*X + (x ) : -1;
421 ju = (z < Z-1) ? (z+1)*X*Y + (y )*X + (x ) : -1;
422
423 if ( ((double)rand() / (double)RAND_MAX) < drop_probability) jw = -1;
424 if ( ((double)rand() / (double)RAND_MAX) < drop_probability) je = -1;
425 if ( ((double)rand() / (double)RAND_MAX) < drop_probability) js = -1;
426 if ( ((double)rand() / (double)RAND_MAX) < drop_probability) jn = -1;
427 if ( ((double)rand() / (double)RAND_MAX) < drop_probability) ju = -1;
428 if ( ((double)rand() / (double)RAND_MAX) < drop_probability) jd = -1;
429
430 /*printf("xyz=%d %d %d j's=%d %d %d %d %d %d %d\n",x,y,z,j,jw,je,js,jn,jd,ju);*/
431 /*printf("column %d xy %d,%d starts at %d\n",j,x,y,ip);*/
432 (m->colptr)[j] = ip;
433 jp = ip;
434
435 /*printf("j=%d D[]=%lf\n",j,D[j]);*/
436
437 (m->rowind)[ip]= j;
438 /*
439 if (x==0 && y==0 && z==0)
440 (m->taucs_values)[ip]= 1.0;
441 else
442 */
443 (m->values.d/*taucs_values*/)[ip]= D[j];
444 ip++;
445
446 if (jw != j != -1) {
447 if (jw > j) {
448 v = -1.0;
449 v = ((double)rand()/(double)RAND_MAX) > 0.99 ? -1.0 : -rmin;
450 v = -( rmin + (((double)rand()/(double)RAND_MAX) * (1.0-rmin)) );
451 /*printf(">> %g\n",v);*/
452 (m->rowind)[ip] = jw;
453 (m->values.d/*taucs_values*/)[ip] = v;
454 ip++;
455 (m->values.d/*taucs_values*/)[jp] -= v;
456 D[jw] -= v;
457 }
458 }
459
460 if (je != j && je != jw && je != -1) {
461 if (je > j) {
462 v = -1.0;
463 v = ((double)rand()/(double)RAND_MAX) > 0.99 ? -1.0 : -rmin;
464 v = -( rmin + (((double)rand()/(double)RAND_MAX) * (1.0-rmin)) );
465 /*printf(">> %g\n",v);*/
466 (m->rowind)[ip] = je;
467 (m->values.d/*taucs_values*/)[ip] = v;
468 ip++;
469 (m->values.d/*taucs_values*/)[jp] -= v;
470 D[je] -= v;
471 }
472 }
473
474 if (js != j && js != -1) {
475 if (js > j) {
476 v = -1.0;
477 v = ((double)rand()/(double)RAND_MAX) > 0.99 ? -1.0 : -rmin;
478 v = -( rmin + (((double)rand()/(double)RAND_MAX) * (1.0-rmin)) );
479 /*printf(">> %g\n",v);*/
480 (m->rowind)[ip] = js;
481 (m->values.d/*taucs_values*/)[ip] = v;
482 ip++;
483 (m->values.d/*taucs_values*/)[jp] -= v;
484 D[js] -= v;
485 }
486 }
487
488 if (jn != j && jn != js && jn != -1) {
489 if (jn > j) {
490 v = -1.0;
491 v = ((double)rand()/(double)RAND_MAX) > 0.99 ? -1.0 : -rmin;
492 v = -( rmin + (((double)rand()/(double)RAND_MAX) * (1.0-rmin)) );
493 /*printf(">> %g\n",v);*/
494 (m->rowind)[ip] = jn;
495 (m->values.d/*taucs_values*/)[ip] = v;
496 ip++;
497 (m->values.d/*taucs_values*/)[jp] -= v;
498 D[jn] -= v;
499 }
500 }
501
502 if (ju != j && ju != -1) {
503 if (ju > j) {
504 v = -1.0;
505 v = ((double)rand()/(double)RAND_MAX) > 0.99 ? -1.0 : -rmin;
506 v = -( rmin + (((double)rand()/(double)RAND_MAX) * (1.0-rmin)) );
507 /*printf(">> %g\n",v);*/
508 (m->rowind)[ip] = ju;
509 (m->values.d/*taucs_values*/)[ip] = v;
510 ip++;
511 (m->values.d/*taucs_values*/)[jp] -= v;
512 D[ju] -= v;
513 }
514 }
515
516 if (jd != j && jd != ju && jd != -1) {
517 if (jd > j) {
518 v = -1.0;
519 v = ((double)rand()/(double)RAND_MAX) > 0.99 ? -1.0 : -rmin;
520 v = -( rmin + (((double)rand()/(double)RAND_MAX) * (1.0-rmin)) );
521 /*printf(">> %g\n",v);*/
522 (m->rowind)[ip] = jd;
523 (m->values.d/*taucs_values*/)[ip] = v;
524 ip++;
525 (m->values.d/*taucs_values*/)[jp] -= v;
526 D[jd] -= v;
527 }
528 }
529
530 }
531 }
532 }
533 taucs_free(D);
534 (m->colptr)[N] = ip;
535
536 taucs_printf("taucs_ccs_generate_rrn: done, ncols=%d allocated nnz=%d real nnz=%d\n",
537 N,nnz,ip);
538
539
540 neighbors = (int**) taucs_malloc(N * sizeof(int*));
541 degree = (int*) taucs_malloc(N * sizeof(int));
542 visited = (int*) taucs_malloc(N * sizeof(int));
543 reps = (int*) taucs_malloc(N * sizeof(int));
544
545 for (i=0; i<N; i++) degree[i] = 0;
546
547 for (j=0; j<N; j++) {
548 for (ip=(m->colptr)[j]; ip<(m->colptr)[j+1]; ip++) {
549 i = (m->rowind)[ ip ];
550 if (i != j) {
551 degree[i]++;
552 degree[j]++;
553 }
554 }
555 }
556
557
558 for (i=0; i<N; i++) {
559 neighbors[i] = (int*) taucs_malloc(degree[i] * sizeof(int));
560 visited[i] = 0;
561 }
562
563 for (j=0; j<N; j++) {
564 for (ip=(m->colptr)[j]; ip<(m->colptr)[j+1]; ip++) {
565 i = (m->rowind)[ ip ];
566 if (i != j) {
567 neighbors[i][visited[i]] = j;
568 neighbors[j][visited[j]] = i;
569 assert(visited[i] < degree[i]);
570 assert(visited[j] < degree[j]);
571 visited[i]++;
572 visited[j]++;
573 }
574 }
575 }
576
577 for (i=0; i<N; i++) visited[i] = 0;
578 ncomponents = 0;
579 largest = -1;
580 for (i=0; i<N; i++) {
581 if (visited[i] == 0) {
582 int count;
583 reps[ncomponents] = i;
584 ncomponents++;
585 count = recursive_visit(i,neighbors,degree,visited);
586 if (count > largest) {
587 largest = count;
588 largest_rep = i;
589 }
590 /*printf("new connected component vertex %d, size=%d\n",i,count);*/
591 }
592 }
593 for (i=0; i<ncomponents; i++) {
594 j = reps[i];
595 /*printf("rep[%d] = %d\n",i,j);*/
596 (m->values.d/*taucs_values*/)[ (m->colptr)[j] ] += 1.0;
597 }
598 printf("found %d components, largest is %d, rep is %d\n",ncomponents,largest,largest_rep);
599 printf("found %d components\n",ncomponents);
600
601 for (i=0; i<N; i++) visited[i] = 0;
602 (void) recursive_visit(largest_rep,neighbors,degree,visited);
603
604 /* we now reuse the degree and reps vectors */
605
606 for (i=0; i<N; i++) degree[i] = reps[i] = -1;
607 j = 0;
608 for (i=0; i<N; i++) {
609 if (visited[i]) {
610 degree[i] = j;
611 reps[j] = i;
612 j++;
613 }
614 }
615
616 l = (taucs_ccs_matrix*) taucs_malloc(sizeof(taucs_ccs_matrix));
617 if (!l) {
618 taucs_printf("taucs_ccs_generate_rrn: out of memory\n");
619 return NULL;
620 }
621
622 nnz = (m->colptr)[N]; /* this is an upper bound */
623
624 l->n = largest;
625 l->m = largest;
626 l->flags = TAUCS_SYMMETRIC | TAUCS_LOWER | TAUCS_DOUBLE;
627 l->colptr = (int*) taucs_malloc((largest+1) * sizeof(int));
628 l->rowind = (int*) taucs_malloc(nnz * sizeof(int));
629 l->values.d/*taucs_values*/ = (double*) taucs_malloc(nnz * sizeof(double));
630
631 k = 0;
632 for (jp=0; jp<N; jp++) {
633 int iip;
634 j = degree[jp];
635 if (j == -1) continue;
636 assert(j < largest);
637 (l->colptr)[j] = k;
638 for (iip=(m->colptr)[jp]; iip<(m->colptr)[jp+1]; iip++) {
639 double v;
640 ip = (m->rowind)[iip];
641 v = (m->values.d/*taucs_values*/)[iip];
642 i = degree[ip];
643 assert(i >= j);
644 (l->rowind)[k] = i;
645 (l->values.d/*taucs_values*/)[k] = v;
646 k++;
647 }
648 }
649 (l->colptr)[largest] = k;
650
651 for (i=0; i<N; i++) taucs_free(neighbors[i]);
652 taucs_free(visited);
653 taucs_free(reps);
654 taucs_free(degree);
655 taucs_free(neighbors);
656
657 taucs_ccs_free(m);
658
659 return l;
660 }
661
taucs_vec_generate_continuous(int X,int Y,int Z,char * which)662 double* taucs_vec_generate_continuous(int X, int Y, int Z, char* which)
663 {
664 int x,y,z,j;/* i,k omer*/
665 double* V;
666 double dx,dy,dz;
667
668 V = (double*) taucs_malloc( X*Y*Z * sizeof(double));
669 if (!V) {
670 taucs_printf("taucs_vec_generate_continuous: out of memory\n");
671 return V;
672 }
673
674 for (z=0; z<Z; z++) {
675 for (y=0; y<Y; y++) {
676 for (x=0; x<X; x++) {
677 double v;
678
679 j = z*X*Y + y*X + x;
680
681 dx = (double) (x+1) / (double) X;
682 dy = (double) (y+1) / (double) Y;
683 dz = (double) (z+1) / (double) Z;
684
685 v = (dx*dy*dz*(1.0-dx)*(1.0-dy)*(1.0-dz));
686 v = v*v;
687 v = v*exp(dx*dx*dy*dz);
688
689 V[j] = v;
690 }
691 }
692 }
693
694 return V;
695 }
696
697 taucs_ccs_matrix*
taucs_ccs_generate_discontinuous(int X,int Y,int Z,double jump)698 taucs_ccs_generate_discontinuous(int X, int Y, int Z, double jump)
699 {
700 taucs_ccs_matrix* m;
701 /*taucs_ccs_matrix* l; omer*/
702 int N;
703 int nnz;
704 int x,y,z,i,ip;/*j,k,jp omer*/
705 double* D; /* contributions to future diagonal elements */
706
707 taucs_printf("taucs_ccs_generate_discontinuous: starting (%d %d %d %e)\n",
708 X,Y,Z,jump);
709
710
711 m = (taucs_ccs_matrix*) taucs_malloc(sizeof(taucs_ccs_matrix));
712 if (!m) {
713 taucs_printf("taucs_ccs_generate_discontinuous: out of memory\n");
714 return NULL;
715 }
716
717 N = X*Y*Z;
718 nnz = 4*N; /* this is an upper bound */
719
720 m->n = N;
721 m->m = N;
722 m->flags = TAUCS_SYMMETRIC | TAUCS_LOWER | TAUCS_DOUBLE;
723 /*m->indshift = 0;*/
724 m->colptr = (int*) taucs_malloc((N+1) * sizeof(int));
725 m->rowind = (int*) taucs_malloc(nnz * sizeof(int));
726 m->values.d/*taucs_values*/ = (double*) taucs_malloc(nnz * sizeof(double));
727
728 D = (double*) taucs_malloc(N * sizeof(double));
729
730 if (!(m->colptr) || !(m->rowind) || !(m->rowind) || !D) {
731 taucs_printf("taucs_ccs_generate_discontinuous: out of memory: ncols=%d nnz=%d\n",N,nnz);
732 taucs_free(m->colptr); taucs_free(m->rowind); taucs_free(m->values.d/*taucs_values*/); taucs_free(D);
733 return NULL;
734 }
735
736 for (i=0; i<N; i++) D[i] = 0.0;
737
738 ip = 0;
739 for (z=0; z<Z; z++) {
740 for (y=0; y<Y; y++) {
741 for (x=0; x<X; x++) {
742 int j, je, jw, js, jn, ju, jd; /* indices for up, down, east, west, south, north */
743 int jp; /* pointer to the diagonal value */
744 double v;
745 int cj, cjw, cje, cjs, cjn, cjd, cju; /* which region? */
746
747 j = z*X*Y + y*X + x;
748 jw = (x > 0 ) ? (z )*X*Y + (y )*X + (x-1) : (z )*X*Y + (y )*X + (X-1) ;
749 je = (x < X-1) ? (z )*X*Y + (y )*X + (x+1) : (z )*X*Y + (y )*X + (0 ) ;
750 js = (y > 0 ) ? (z )*X*Y + (y-1)*X + (x ) : (z )*X*Y + (Y-1)*X + (x ) ;
751 jn = (y < Y-1) ? (z )*X*Y + (y+1)*X + (x ) : (z )*X*Y + (0 )*X + (x ) ;
752 jd = (z > 0 ) ? (z-1)*X*Y + (y )*X + (x ) : (Z-1)*X*Y + (y )*X + (x ) ;
753 ju = (z < Z-1) ? (z+1)*X*Y + (y )*X + (x ) : (0 )*X*Y + (y )*X + (x ) ;
754
755 jw = (x > 0 ) ? (z )*X*Y + (y )*X + (x-1) : -1;
756 je = (x < X-1) ? (z )*X*Y + (y )*X + (x+1) : -1;
757 js = (y > 0 ) ? (z )*X*Y + (y-1)*X + (x ) : -1;
758 jn = (y < Y-1) ? (z )*X*Y + (y+1)*X + (x ) : -1;
759 jd = (z > 0 ) ? (z-1)*X*Y + (y )*X + (x ) : -1;
760 ju = (z < Z-1) ? (z+1)*X*Y + (y )*X + (x ) : -1;
761
762 /*printf("xyz=%d %d %d j's=%d %d %d %d %d %d %d\n",x,y,z,j,jw,je,js,jn,jd,ju);*/
763 /*printf("column %d xy %d,%d starts at %d\n",j,x,y,ip);*/
764 (m->colptr)[j] = ip;
765 jp = ip;
766
767 /*printf("j=%d D[]=%lf\n",j,D[j]);*/
768
769 (m->rowind)[ip]= j;
770 /* Nonsingular Neumann */
771 if (x==0 && y==0 && z==0)
772 (m->values.d/*taucs_values*/)[ip]= D[j] + 1.0;
773 else
774 (m->values.d/*taucs_values*/)[ip]= D[j];
775
776 /* Singular Neumann */
777 /*
778 (m->taucs_values)[ip] = D[j];
779 */
780
781 /* Dirichlet */
782 /*
783 (m->taucs_values)[ip] = D[j];
784 if (x==0 || x==X-1) (m->taucs_values)[ip] += 1.0;
785 if (y==0 || y==Y-1) (m->taucs_values)[ip] += 1.0;
786 if (z==0 || z==Z-1) (m->taucs_values)[ip] += 1.0;
787 */
788
789 ip++;
790
791 cj = ((x ) >= X/8 && (x ) < 7*X/8)
792 && ((y ) >= Y/8 && (y ) < 7*Y/8)
793 && ((z ) >= Z/8 && (z ) < 7*Z/8);
794 /*
795 cj = cj && !( ((x ) >= 2*X/8 && (x ) < 6*X/8)
796 && ((y ) >= 2*Y/8 && (y ) < 6*Y/8)
797 && ((z ) >= 2*Z/8 && (z ) < 6*Z/8));
798 */
799 cjw = ((x-1) >= X/8 && (x-1) < 7*X/8)
800 && ((y ) >= Y/8 && (y ) < 7*Y/8)
801 && ((z ) >= Z/8 && (z ) < 7*Z/8);
802 /*
803 cjw = cjw && !( ((x-1) >= 2*X/8 && (x-1) < 6*X/8)
804 && ((y ) >= 2*Y/8 && (y ) < 6*Y/8)
805 && ((z ) >= 2*Z/8 && (z ) < 6*Z/8));
806 */
807 cje = ((x+1) >= X/8 && (x+1) < 7*X/8)
808 && ((y ) >= Y/8 && (y ) < 7*Y/8)
809 && ((z ) >= Z/8 && (z ) < 7*Z/8);
810 /*
811 cje = cje && !( ((x+1) >= 2*X/8 && (x+1) < 6*X/8)
812 && ((y ) >= 2*Y/8 && (y ) < 6*Y/8)
813 && ((z ) >= 2*Z/8 && (z ) < 6*Z/8));
814 */
815 cjs = ((x ) >= X/8 && (x ) < 7*X/8)
816 && ((y-1) >= Y/8 && (y-1) < 7*Y/8)
817 && ((z ) >= Z/8 && (z ) < 7*Z/8);
818 /*
819 cjs = cjs && !( ((x ) >= 2*X/8 && (x ) < 6*X/8)
820 && ((y-1) >= 2*Y/8 && (y-1) < 6*Y/8)
821 && ((z ) >= 2*Z/8 && (z ) < 6*Z/8));
822 */
823 cjn = ((x ) >= X/8 && (x ) < 7*X/8)
824 && ((y+1) >= Y/8 && (y+1) < 7*Y/8)
825 && ((z ) >= Z/8 && (z ) < 7*Z/8);
826 /*
827 cjn = cjn && !( ((x ) >= 2*X/8 && (x ) < 6*X/8)
828 && ((y+1) >= 2*Y/8 && (y+1) < 6*Y/8)
829 && ((z ) >= 2*Z/8 && (z ) < 6*Z/8));
830 */
831 cjd = ((x ) >= X/8 && (x ) < 7*X/8)
832 && ((y ) >= Y/8 && (y ) < 7*Y/8)
833 && ((z-1) >= Z/8 && (z-1) < 7*Z/8);
834 /*
835 cjd = cjd && !( ((x ) >= 2*X/8 && (x ) < 6*X/8)
836 && ((y ) >= 2*Y/8 && (y ) < 6*Y/8)
837 && ((z-1) >= 2*Z/8 && (z-1) < 6*Z/8));
838 */
839 cju = ((x ) >= X/8 && (x ) < 7*X/8)
840 && ((y ) >= Y/8 && (y ) < 7*Y/8)
841 && ((z+1) >= Z/8 && (z+1) < 7*Z/8);
842 /*
843 cju = cju && !( ((x ) >= 2*X/8 && (x ) < 6*X/8)
844 && ((y ) >= 2*Y/8 && (y ) < 6*Y/8)
845 && ((z+1) >= 2*Z/8 && (z+1) < 6*Z/8));
846 */
847
848 if (jw != j && jw != -1) {
849 if (jw > j) {
850 v = -jump;
851 v = (x < X/8 || y < Y/8) ? -jump : -1.0;
852 v = (cj && cjw) ? -jump : -1.0;
853 /*printf(">> %g\n",v);*/
854 (m->rowind)[ip] = jw;
855 (m->values.d/*taucs_values*/)[ip] = v;
856 ip++;
857 (m->values.d/*taucs_values*/)[jp] -= v;
858 D[jw] -= v;
859 }
860 }
861
862 if (je != j && je != jw && je != -1) {
863 if (je > j) {
864 v = -jump;
865 v = ((x-1) < X/8 || y < Y/8) ? -jump : -1.0;
866 v = (cj && cje) ? -jump : -1.0;
867 /*printf(">> %g\n",v);*/
868 (m->rowind)[ip] = je;
869 (m->values.d/*taucs_values*/)[ip] = v;
870 ip++;
871 (m->values.d/*taucs_values*/)[jp] -= v;
872 D[je] -= v;
873 }
874 }
875
876 if (js != j && js != -1) {
877 if (js > j) {
878 v = -jump;
879 v = (y < Y/8 || x < X/8) ? -jump : -1.0;
880 v = (cj && cjs) ? -jump : -1.0;
881 /*printf(">> %g\n",v);*/
882 (m->rowind)[ip] = js;
883 (m->values.d/*taucs_values*/)[ip] = v;
884 ip++;
885 (m->values.d/*taucs_values*/)[jp] -= v;
886 D[js] -= v;
887 }
888 }
889
890 if (jn != j && jn != js && jn != -1) {
891 if (jn > j) {
892 v = -jump;
893 v = ((y-1) < Y/8 || x < X/8) ? -jump : -1.0;
894 v = (cj && cjn) ? -jump : -1.0;
895 /*printf(">> %g\n",v);*/
896 (m->rowind)[ip] = jn;
897 (m->values.d/*taucs_values*/)[ip] = v;
898 ip++;
899 (m->values.d/*taucs_values*/)[jp] -= v;
900 D[jn] -= v;
901 }
902 }
903
904 if (ju != j && ju != -1) {
905 if (ju > j) {
906 v = -1.0;
907 v = (cj && cju) ? -jump : -1.0;
908 /*printf(">> %g\n",v);*/
909 (m->rowind)[ip] = ju;
910 (m->values.d/*taucs_values*/)[ip] = v;
911 ip++;
912 (m->values.d/*taucs_values*/)[jp] -= v;
913 D[ju] -= v;
914 }
915 }
916
917 if (jd != j && jd != ju && jd != -1) {
918 if (jd > j) {
919 v = -1.0;
920 v = (cj && cjd) ? -jump : -1.0;
921 /*printf(">> %g\n",v);*/
922 (m->rowind)[ip] = jd;
923 (m->values.d/*taucs_values*/)[ip] = v;
924 ip++;
925 (m->values.d/*taucs_values*/)[jp] -= v;
926 D[jd] -= v;
927 }
928 }
929
930 }
931 }
932 }
933 taucs_free(D);
934 (m->colptr)[N] = ip;
935
936 taucs_printf("taucs_ccs_generate_discontinuous: done, ncols=%d allocated nnz=%d real nnz=%d\n",
937 N,nnz,ip);
938
939 /*taucs_ccs_write_ijv(m,"X.ijv");*/
940
941 return m;
942 }
943
944 #endif /* TAUCS_CORE_DOUBLE */
945
946 /*********************************************************/
947 /* */
948 /*********************************************************/
949