Actual source code: spbas.c
1: #include <../src/mat/impls/aij/seq/aij.h>
2: #include <../src/mat/impls/aij/seq/bas/spbas.h>
4: /*MC
5: MATSOLVERBAS - Provides ICC(k) with drop tolerance
7: Works with MATAIJ matrices
9: Options Database Keys:
10: + -pc_factor_levels <l> - number of levels of fill
11: - -pc_factor_drop_tolerance - is not currently hooked up to do anything
13: Level: intermediate
15: Contributed by: Bas van 't Hof
17: Notes:
18: Since this currently hooked up to use drop tolerance it should produce the same factors and hence convergence as the PETSc ICC, for higher
19: levels of fill it does not. This needs to be investigated. Unless you are interested in drop tolerance ICC and willing to work through the code
20: we recommend not using this functionality.
22: .seealso: PCFactorSetMatSolverType(), MatSolverType, PCFactorSetLevels(), PCFactorSetDropTolerance()
24: M*/
26: /*
27: spbas_memory_requirement:
28: Calculate the number of bytes needed to store tha matrix
29: */
30: size_t spbas_memory_requirement(spbas_matrix matrix)
31: {
32: size_t memreq = 6 * sizeof(PetscInt) + /* nrows, ncols, nnz, n_alloc_icol, n_alloc_val, col_idx_type */
33: sizeof(PetscBool) + /* block_data */
34: sizeof(PetscScalar**) + /* values */
35: sizeof(PetscScalar*) + /* alloc_val */
36: 2 * sizeof(PetscInt**) + /* icols, icols0 */
37: 2 * sizeof(PetscInt*) + /* row_nnz, alloc_icol */
38: matrix.nrows * sizeof(PetscInt) + /* row_nnz[*] */
39: matrix.nrows * sizeof(PetscInt*); /* icols[*] */
41: /* icol0[*] */
42: if (matrix.col_idx_type == SPBAS_OFFSET_ARRAY) memreq += matrix.nrows * sizeof(PetscInt);
44: /* icols[*][*] */
45: if (matrix.block_data) memreq += matrix.n_alloc_icol * sizeof(PetscInt);
46: else memreq += matrix.nnz * sizeof(PetscInt);
48: if (matrix.values) {
49: memreq += matrix.nrows * sizeof(PetscScalar*); /* values[*] */
50: /* values[*][*] */
51: if (matrix.block_data) memreq += matrix.n_alloc_val * sizeof(PetscScalar);
52: else memreq += matrix.nnz * sizeof(PetscScalar);
53: }
54: return memreq;
55: }
57: /*
58: spbas_allocate_pattern:
59: allocate the pattern arrays row_nnz, icols and optionally values
60: */
61: PetscErrorCode spbas_allocate_pattern(spbas_matrix * result, PetscBool do_values)
62: {
63: PetscInt nrows = result->nrows;
64: PetscInt col_idx_type = result->col_idx_type;
66: /* Allocate sparseness pattern */
67: PetscMalloc1(nrows,&result->row_nnz);
68: PetscMalloc1(nrows,&result->icols);
70: /* If offsets are given wrt an array, create array */
71: if (col_idx_type == SPBAS_OFFSET_ARRAY) {
72: PetscMalloc1(nrows,&result->icol0);
73: } else {
74: result->icol0 = NULL;
75: }
77: /* If values are given, allocate values array */
78: if (do_values) {
79: PetscMalloc1(nrows,&result->values);
80: } else {
81: result->values = NULL;
82: }
83: return 0;
84: }
86: /*
87: spbas_allocate_data:
88: in case of block_data:
89: Allocate the data arrays alloc_icol and optionally alloc_val,
90: set appropriate pointers from icols and values;
91: in case of !block_data:
92: Allocate the arrays icols[i] and optionally values[i]
93: */
94: PetscErrorCode spbas_allocate_data(spbas_matrix * result)
95: {
96: PetscInt i;
97: PetscInt nnz = result->nnz;
98: PetscInt nrows = result->nrows;
99: PetscInt r_nnz;
100: PetscBool do_values = (result->values) ? PETSC_TRUE : PETSC_FALSE;
101: PetscBool block_data = result->block_data;
103: if (block_data) {
104: /* Allocate the column number array and point to it */
105: result->n_alloc_icol = nnz;
107: PetscMalloc1(nnz, &result->alloc_icol);
109: result->icols[0] = result->alloc_icol;
110: for (i=1; i<nrows; i++) {
111: result->icols[i] = result->icols[i-1] + result->row_nnz[i-1];
112: }
114: /* Allocate the value array and point to it */
115: if (do_values) {
116: result->n_alloc_val = nnz;
118: PetscMalloc1(nnz, &result->alloc_val);
120: result->values[0] = result->alloc_val;
121: for (i=1; i<nrows; i++) {
122: result->values[i] = result->values[i-1] + result->row_nnz[i-1];
123: }
124: }
125: } else {
126: for (i=0; i<nrows; i++) {
127: r_nnz = result->row_nnz[i];
128: PetscMalloc1(r_nnz, &result->icols[i]);
129: }
130: if (do_values) {
131: for (i=0; i<nrows; i++) {
132: r_nnz = result->row_nnz[i];
133: PetscMalloc1(r_nnz, &result->values[i]);
134: }
135: }
136: }
137: return 0;
138: }
140: /*
141: spbas_row_order_icol
142: determine if row i1 should come
143: + before row i2 in the sorted rows (return -1),
144: + after (return 1)
145: + is identical (return 0).
146: */
147: int spbas_row_order_icol(PetscInt i1, PetscInt i2, PetscInt *irow_in, PetscInt *icol_in,PetscInt col_idx_type)
148: {
149: PetscInt j;
150: PetscInt nnz1 = irow_in[i1+1] - irow_in[i1];
151: PetscInt nnz2 = irow_in[i2+1] - irow_in[i2];
152: PetscInt * icol1 = &icol_in[irow_in[i1]];
153: PetscInt * icol2 = &icol_in[irow_in[i2]];
155: if (nnz1<nnz2) return -1;
156: if (nnz1>nnz2) return 1;
158: if (col_idx_type == SPBAS_COLUMN_NUMBERS) {
159: for (j=0; j<nnz1; j++) {
160: if (icol1[j]< icol2[j]) return -1;
161: if (icol1[j]> icol2[j]) return 1;
162: }
163: } else if (col_idx_type == SPBAS_DIAGONAL_OFFSETS) {
164: for (j=0; j<nnz1; j++) {
165: if (icol1[j]-i1< icol2[j]-i2) return -1;
166: if (icol1[j]-i1> icol2[j]-i2) return 1;
167: }
168: } else if (col_idx_type == SPBAS_OFFSET_ARRAY) {
169: for (j=1; j<nnz1; j++) {
170: if (icol1[j]-icol1[0] < icol2[j]-icol2[0]) return -1;
171: if (icol1[j]-icol1[0] > icol2[j]-icol2[0]) return 1;
172: }
173: }
174: return 0;
175: }
177: /*
178: spbas_mergesort_icols:
179: return a sorting of the rows in which identical sparseness patterns are
180: next to each other
181: */
182: PetscErrorCode spbas_mergesort_icols(PetscInt nrows, PetscInt * irow_in, PetscInt * icol_in,PetscInt col_idx_type, PetscInt *isort)
183: {
184: PetscInt istep; /* Chunk-sizes of already sorted parts of arrays */
185: PetscInt i, i1, i2; /* Loop counters for (partly) sorted arrays */
186: PetscInt istart, i1end, i2end; /* start of newly sorted array part, end of both parts */
187: PetscInt *ialloc; /* Allocated arrays */
188: PetscInt *iswap; /* auxiliary pointers for swapping */
189: PetscInt *ihlp1; /* Pointers to new version of arrays, */
190: PetscInt *ihlp2; /* Pointers to previous version of arrays, */
192: PetscMalloc1(nrows,&ialloc);
194: ihlp1 = ialloc;
195: ihlp2 = isort;
197: /* Sorted array chunks are first 1 long, and increase until they are the complete array */
198: for (istep=1; istep<nrows; istep*=2) {
199: /*
200: Combine sorted parts
201: istart:istart+istep-1 and istart+istep-1:istart+2*istep-1
202: of ihlp2 and vhlp2
204: into one sorted part
205: istart:istart+2*istep-1
206: of ihlp1 and vhlp1
207: */
208: for (istart=0; istart<nrows; istart+=2*istep) {
209: /* Set counters and bound array part endings */
210: i1=istart; i1end = i1+istep; if (i1end>nrows) i1end=nrows;
211: i2=istart+istep; i2end = i2+istep; if (i2end>nrows) i2end=nrows;
213: /* Merge the two array parts */
214: for (i=istart; i<i2end; i++) {
215: if (i1<i1end && i2<i2end && spbas_row_order_icol(ihlp2[i1], ihlp2[i2], irow_in, icol_in, col_idx_type) < 0) {
216: ihlp1[i] = ihlp2[i1];
217: i1++;
218: } else if (i2<i2end) {
219: ihlp1[i] = ihlp2[i2];
220: i2++;
221: } else {
222: ihlp1[i] = ihlp2[i1];
223: i1++;
224: }
225: }
226: }
228: /* Swap the two array sets */
229: iswap = ihlp2; ihlp2 = ihlp1; ihlp1 = iswap;
230: }
232: /* Copy one more time in case the sorted arrays are the temporary ones */
233: if (ihlp2 != isort) {
234: for (i=0; i<nrows; i++) isort[i] = ihlp2[i];
235: }
236: PetscFree(ialloc);
237: return 0;
238: }
240: /*
241: spbas_compress_pattern:
242: calculate a compressed sparseness pattern for a sparseness pattern
243: given in compressed row storage. The compressed sparseness pattern may
244: require (much) less memory.
245: */
246: PetscErrorCode spbas_compress_pattern(PetscInt *irow_in, PetscInt *icol_in, PetscInt nrows, PetscInt ncols, PetscInt col_idx_type, spbas_matrix *B,PetscReal *mem_reduction)
247: {
248: PetscInt nnz = irow_in[nrows];
249: size_t mem_orig = (nrows + nnz) * sizeof(PetscInt);
250: size_t mem_compressed;
251: PetscInt *isort;
252: PetscInt *icols;
253: PetscInt row_nnz;
254: PetscInt *ipoint;
255: PetscBool *used;
256: PetscInt ptr;
257: PetscInt i,j;
258: const PetscBool no_values = PETSC_FALSE;
260: /* Allocate the structure of the new matrix */
261: B->nrows = nrows;
262: B->ncols = ncols;
263: B->nnz = nnz;
264: B->col_idx_type = col_idx_type;
265: B->block_data = PETSC_TRUE;
267: spbas_allocate_pattern(B, no_values);
269: /* When using an offset array, set it */
270: if (col_idx_type==SPBAS_OFFSET_ARRAY) {
271: for (i=0; i<nrows; i++) B->icol0[i] = icol_in[irow_in[i]];
272: }
274: /* Allocate the ordering for the rows */
275: PetscMalloc1(nrows,&isort);
276: PetscMalloc1(nrows,&ipoint);
277: PetscCalloc1(nrows,&used);
279: for (i = 0; i<nrows; i++) {
280: B->row_nnz[i] = irow_in[i+1]-irow_in[i];
281: isort[i] = i;
282: ipoint[i] = i;
283: }
285: /* Sort the rows so that identical columns will be next to each other */
286: spbas_mergesort_icols(nrows, irow_in, icol_in, col_idx_type, isort);
287: PetscInfo(NULL,"Rows have been sorted for patterns\n");
289: /* Replace identical rows with the first one in the list */
290: for (i=1; i<nrows; i++) {
291: if (spbas_row_order_icol(isort[i-1], isort[i], irow_in, icol_in, col_idx_type) == 0) {
292: ipoint[isort[i]] = ipoint[isort[i-1]];
293: }
294: }
296: /* Collect the rows which are used*/
297: for (i=0; i<nrows; i++) used[ipoint[i]] = PETSC_TRUE;
299: /* Calculate needed memory */
300: B->n_alloc_icol = 0;
301: for (i=0; i<nrows; i++) {
302: if (used[i]) B->n_alloc_icol += B->row_nnz[i];
303: }
304: PetscMalloc1(B->n_alloc_icol,&B->alloc_icol);
306: /* Fill in the diagonal offsets for the rows which store their own data */
307: ptr = 0;
308: for (i=0; i<B->nrows; i++) {
309: if (used[i]) {
310: B->icols[i] = &B->alloc_icol[ptr];
311: icols = &icol_in[irow_in[i]];
312: row_nnz = B->row_nnz[i];
313: if (col_idx_type == SPBAS_COLUMN_NUMBERS) {
314: for (j=0; j<row_nnz; j++) {
315: B->icols[i][j] = icols[j];
316: }
317: } else if (col_idx_type == SPBAS_DIAGONAL_OFFSETS) {
318: for (j=0; j<row_nnz; j++) {
319: B->icols[i][j] = icols[j]-i;
320: }
321: } else if (col_idx_type == SPBAS_OFFSET_ARRAY) {
322: for (j=0; j<row_nnz; j++) {
323: B->icols[i][j] = icols[j]-icols[0];
324: }
325: }
326: ptr += B->row_nnz[i];
327: }
328: }
330: /* Point to the right places for all data */
331: for (i=0; i<nrows; i++) {
332: B->icols[i] = B->icols[ipoint[i]];
333: }
334: PetscInfo(NULL,"Row patterns have been compressed\n");
335: PetscInfo(NULL," (%g nonzeros per row)\n", (double) ((PetscReal) nnz / (PetscReal) nrows));
337: PetscFree(isort);
338: PetscFree(used);
339: PetscFree(ipoint);
341: mem_compressed = spbas_memory_requirement(*B);
342: *mem_reduction = 100.0 * (PetscReal)(mem_orig-mem_compressed)/ (PetscReal) mem_orig;
343: return 0;
344: }
346: /*
347: spbas_incomplete_cholesky
348: Incomplete Cholesky decomposition
349: */
350: #include <../src/mat/impls/aij/seq/bas/spbas_cholesky.h>
352: /*
353: spbas_delete : de-allocate the arrays owned by this matrix
354: */
355: PetscErrorCode spbas_delete(spbas_matrix matrix)
356: {
357: PetscInt i;
359: if (matrix.block_data) {
360: PetscFree(matrix.alloc_icol);
361: if (matrix.values) PetscFree(matrix.alloc_val);
362: } else {
363: for (i=0; i<matrix.nrows; i++) PetscFree(matrix.icols[i]);
364: PetscFree(matrix.icols);
365: if (matrix.values) {
366: for (i=0; i<matrix.nrows; i++) PetscFree(matrix.values[i]);
367: }
368: }
370: PetscFree(matrix.row_nnz);
371: PetscFree(matrix.icols);
372: if (matrix.col_idx_type == SPBAS_OFFSET_ARRAY) PetscFree(matrix.icol0);
373: PetscFree(matrix.values);
374: return 0;
375: }
377: /*
378: spbas_matrix_to_crs:
379: Convert an spbas_matrix to compessed row storage
380: */
381: PetscErrorCode spbas_matrix_to_crs(spbas_matrix matrix_A,MatScalar **val_out, PetscInt **irow_out, PetscInt **icol_out)
382: {
383: PetscInt nrows = matrix_A.nrows;
384: PetscInt nnz = matrix_A.nnz;
385: PetscInt i,j,r_nnz,i0;
386: PetscInt *irow;
387: PetscInt *icol;
388: PetscInt *icol_A;
389: MatScalar *val;
390: PetscScalar *val_A;
391: PetscInt col_idx_type = matrix_A.col_idx_type;
392: PetscBool do_values = matrix_A.values ? PETSC_TRUE : PETSC_FALSE;
394: PetscMalloc1(nrows+1, &irow);
395: PetscMalloc1(nnz, &icol);
396: *icol_out = icol;
397: *irow_out = irow;
398: if (do_values) {
399: PetscMalloc1(nnz, &val);
400: *val_out = val; *icol_out = icol; *irow_out=irow;
401: }
403: irow[0]=0;
404: for (i=0; i<nrows; i++) {
405: r_nnz = matrix_A.row_nnz[i];
406: i0 = irow[i];
407: irow[i+1] = i0 + r_nnz;
408: icol_A = matrix_A.icols[i];
410: if (do_values) {
411: val_A = matrix_A.values[i];
412: for (j=0; j<r_nnz; j++) {
413: icol[i0+j] = icol_A[j];
414: val[i0+j] = val_A[j];
415: }
416: } else {
417: for (j=0; j<r_nnz; j++) icol[i0+j] = icol_A[j];
418: }
420: if (col_idx_type == SPBAS_DIAGONAL_OFFSETS) {
421: for (j=0; j<r_nnz; j++) icol[i0+j] += i;
422: } else if (col_idx_type == SPBAS_OFFSET_ARRAY) {
423: i0 = matrix_A.icol0[i];
424: for (j=0; j<r_nnz; j++) icol[i0+j] += i0;
425: }
426: }
427: return 0;
428: }
430: /*
431: spbas_transpose
432: return the transpose of a matrix
433: */
434: PetscErrorCode spbas_transpose(spbas_matrix in_matrix, spbas_matrix * result)
435: {
436: PetscInt col_idx_type = in_matrix.col_idx_type;
437: PetscInt nnz = in_matrix.nnz;
438: PetscInt ncols = in_matrix.nrows;
439: PetscInt nrows = in_matrix.ncols;
440: PetscInt i,j,k;
441: PetscInt r_nnz;
442: PetscInt *irow;
443: PetscInt icol0 = 0;
444: PetscScalar * val;
446: /* Copy input values */
447: result->nrows = nrows;
448: result->ncols = ncols;
449: result->nnz = nnz;
450: result->col_idx_type = SPBAS_COLUMN_NUMBERS;
451: result->block_data = PETSC_TRUE;
453: /* Allocate sparseness pattern */
454: spbas_allocate_pattern(result, in_matrix.values ? PETSC_TRUE : PETSC_FALSE);
456: /* Count the number of nonzeros in each row */
457: for (i = 0; i<nrows; i++) result->row_nnz[i] = 0;
459: for (i=0; i<ncols; i++) {
460: r_nnz = in_matrix.row_nnz[i];
461: irow = in_matrix.icols[i];
462: if (col_idx_type == SPBAS_COLUMN_NUMBERS) {
463: for (j=0; j<r_nnz; j++) result->row_nnz[irow[j]]++;
464: } else if (col_idx_type == SPBAS_DIAGONAL_OFFSETS) {
465: for (j=0; j<r_nnz; j++) result->row_nnz[i+irow[j]]++;
466: } else if (col_idx_type == SPBAS_OFFSET_ARRAY) {
467: icol0=in_matrix.icol0[i];
468: for (j=0; j<r_nnz; j++) result->row_nnz[icol0+irow[j]]++;
469: }
470: }
472: /* Set the pointers to the data */
473: spbas_allocate_data(result);
475: /* Reset the number of nonzeros in each row */
476: for (i = 0; i<nrows; i++) result->row_nnz[i] = 0;
478: /* Fill the data arrays */
479: if (in_matrix.values) {
480: for (i=0; i<ncols; i++) {
481: r_nnz = in_matrix.row_nnz[i];
482: irow = in_matrix.icols[i];
483: val = in_matrix.values[i];
485: if (col_idx_type == SPBAS_COLUMN_NUMBERS) icol0 = 0;
486: else if (col_idx_type == SPBAS_DIAGONAL_OFFSETS) icol0 = i;
487: else if (col_idx_type == SPBAS_OFFSET_ARRAY) icol0 = in_matrix.icol0[i];
488: for (j=0; j<r_nnz; j++) {
489: k = icol0 + irow[j];
490: result->icols[k][result->row_nnz[k]] = i;
491: result->values[k][result->row_nnz[k]] = val[j];
492: result->row_nnz[k]++;
493: }
494: }
495: } else {
496: for (i=0; i<ncols; i++) {
497: r_nnz = in_matrix.row_nnz[i];
498: irow = in_matrix.icols[i];
500: if (col_idx_type == SPBAS_COLUMN_NUMBERS) icol0=0;
501: else if (col_idx_type == SPBAS_DIAGONAL_OFFSETS) icol0=i;
502: else if (col_idx_type == SPBAS_OFFSET_ARRAY) icol0=in_matrix.icol0[i];
504: for (j=0; j<r_nnz; j++) {
505: k = icol0 + irow[j];
506: result->icols[k][result->row_nnz[k]] = i;
507: result->row_nnz[k]++;
508: }
509: }
510: }
511: return 0;
512: }
514: /*
515: spbas_mergesort
517: mergesort for an array of integers and an array of associated
518: reals
520: on output, icol[0..nnz-1] is increasing;
521: val[0..nnz-1] has undergone the same permutation as icol
523: NB: val may be NULL: in that case, only the integers are sorted
525: */
526: PetscErrorCode spbas_mergesort(PetscInt nnz, PetscInt *icol, PetscScalar *val)
527: {
528: PetscInt istep; /* Chunk-sizes of already sorted parts of arrays */
529: PetscInt i, i1, i2; /* Loop counters for (partly) sorted arrays */
530: PetscInt istart, i1end, i2end; /* start of newly sorted array part, end of both parts */
531: PetscInt *ialloc; /* Allocated arrays */
532: PetscScalar *valloc=NULL;
533: PetscInt *iswap; /* auxiliary pointers for swapping */
534: PetscScalar *vswap;
535: PetscInt *ihlp1; /* Pointers to new version of arrays, */
536: PetscScalar *vhlp1=NULL; /* (arrays under construction) */
537: PetscInt *ihlp2; /* Pointers to previous version of arrays, */
538: PetscScalar *vhlp2=NULL;
540: PetscMalloc1(nnz,&ialloc);
541: ihlp1 = ialloc;
542: ihlp2 = icol;
544: if (val) {
545: PetscMalloc1(nnz,&valloc);
546: vhlp1 = valloc;
547: vhlp2 = val;
548: }
550: /* Sorted array chunks are first 1 long, and increase until they are the complete array */
551: for (istep=1; istep<nnz; istep*=2) {
552: /*
553: Combine sorted parts
554: istart:istart+istep-1 and istart+istep-1:istart+2*istep-1
555: of ihlp2 and vhlp2
557: into one sorted part
558: istart:istart+2*istep-1
559: of ihlp1 and vhlp1
560: */
561: for (istart=0; istart<nnz; istart+=2*istep) {
562: /* Set counters and bound array part endings */
563: i1=istart; i1end = i1+istep; if (i1end>nnz) i1end=nnz;
564: i2=istart+istep; i2end = i2+istep; if (i2end>nnz) i2end=nnz;
566: /* Merge the two array parts */
567: if (val) {
568: for (i=istart; i<i2end; i++) {
569: if (i1<i1end && i2<i2end && ihlp2[i1] < ihlp2[i2]) {
570: ihlp1[i] = ihlp2[i1];
571: vhlp1[i] = vhlp2[i1];
572: i1++;
573: } else if (i2<i2end) {
574: ihlp1[i] = ihlp2[i2];
575: vhlp1[i] = vhlp2[i2];
576: i2++;
577: } else {
578: ihlp1[i] = ihlp2[i1];
579: vhlp1[i] = vhlp2[i1];
580: i1++;
581: }
582: }
583: } else {
584: for (i=istart; i<i2end; i++) {
585: if (i1<i1end && i2<i2end && ihlp2[i1] < ihlp2[i2]) {
586: ihlp1[i] = ihlp2[i1];
587: i1++;
588: } else if (i2<i2end) {
589: ihlp1[i] = ihlp2[i2];
590: i2++;
591: } else {
592: ihlp1[i] = ihlp2[i1];
593: i1++;
594: }
595: }
596: }
597: }
599: /* Swap the two array sets */
600: iswap = ihlp2; ihlp2 = ihlp1; ihlp1 = iswap;
601: vswap = vhlp2; vhlp2 = vhlp1; vhlp1 = vswap;
602: }
604: /* Copy one more time in case the sorted arrays are the temporary ones */
605: if (ihlp2 != icol) {
606: for (i=0; i<nnz; i++) icol[i] = ihlp2[i];
607: if (val) {
608: for (i=0; i<nnz; i++) val[i] = vhlp2[i];
609: }
610: }
612: PetscFree(ialloc);
613: if (val) PetscFree(valloc);
614: return 0;
615: }
617: /*
618: spbas_apply_reordering_rows:
619: apply the given reordering to the rows: matrix_A = matrix_A(perm,:);
620: */
621: PetscErrorCode spbas_apply_reordering_rows(spbas_matrix *matrix_A, const PetscInt *permutation)
622: {
623: PetscInt i,j,ip;
624: PetscInt nrows=matrix_A->nrows;
625: PetscInt * row_nnz;
626: PetscInt **icols;
627: PetscBool do_values = matrix_A->values ? PETSC_TRUE : PETSC_FALSE;
628: PetscScalar **vals = NULL;
632: if (do_values) {
633: PetscMalloc1(nrows, &vals);
634: }
635: PetscMalloc1(nrows, &row_nnz);
636: PetscMalloc1(nrows, &icols);
638: for (i=0; i<nrows; i++) {
639: ip = permutation[i];
640: if (do_values) vals[i] = matrix_A->values[ip];
641: icols[i] = matrix_A->icols[ip];
642: row_nnz[i] = matrix_A->row_nnz[ip];
643: for (j=0; j<row_nnz[i]; j++) icols[i][j] += ip-i;
644: }
646: if (do_values) PetscFree(matrix_A->values);
647: PetscFree(matrix_A->icols);
648: PetscFree(matrix_A->row_nnz);
650: if (do_values) matrix_A->values = vals;
651: matrix_A->icols = icols;
652: matrix_A->row_nnz = row_nnz;
653: return 0;
654: }
656: /*
657: spbas_apply_reordering_cols:
658: apply the given reordering to the columns: matrix_A(:,perm) = matrix_A;
659: */
660: PetscErrorCode spbas_apply_reordering_cols(spbas_matrix *matrix_A,const PetscInt *permutation)
661: {
662: PetscInt i,j;
663: PetscInt nrows=matrix_A->nrows;
664: PetscInt row_nnz;
665: PetscInt *icols;
666: PetscBool do_values = matrix_A->values ? PETSC_TRUE : PETSC_FALSE;
667: PetscScalar *vals = NULL;
671: for (i=0; i<nrows; i++) {
672: icols = matrix_A->icols[i];
673: row_nnz = matrix_A->row_nnz[i];
674: if (do_values) vals = matrix_A->values[i];
676: for (j=0; j<row_nnz; j++) {
677: icols[j] = permutation[i+icols[j]]-i;
678: }
679: spbas_mergesort(row_nnz, icols, vals);
680: }
681: return 0;
682: }
684: /*
685: spbas_apply_reordering:
686: apply the given reordering: matrix_A(perm,perm) = matrix_A;
687: */
688: PetscErrorCode spbas_apply_reordering(spbas_matrix *matrix_A, const PetscInt *permutation, const PetscInt * inv_perm)
689: {
690: spbas_apply_reordering_rows(matrix_A, inv_perm);
691: spbas_apply_reordering_cols(matrix_A, permutation);
692: return 0;
693: }
695: PetscErrorCode spbas_pattern_only(PetscInt nrows, PetscInt ncols, PetscInt *ai, PetscInt *aj, spbas_matrix * result)
696: {
697: spbas_matrix retval;
698: PetscInt i, j, i0, r_nnz;
700: /* Copy input values */
701: retval.nrows = nrows;
702: retval.ncols = ncols;
703: retval.nnz = ai[nrows];
705: retval.block_data = PETSC_TRUE;
706: retval.col_idx_type = SPBAS_DIAGONAL_OFFSETS;
708: /* Allocate output matrix */
709: spbas_allocate_pattern(&retval, PETSC_FALSE);
710: for (i=0; i<nrows; i++) retval.row_nnz[i] = ai[i+1]-ai[i];
711: spbas_allocate_data(&retval);
712: /* Copy the structure */
713: for (i = 0; i<retval.nrows; i++) {
714: i0 = ai[i];
715: r_nnz = ai[i+1]-i0;
717: for (j=0; j<r_nnz; j++) {
718: retval.icols[i][j] = aj[i0+j]-i;
719: }
720: }
721: *result = retval;
722: return 0;
723: }
725: /*
726: spbas_mark_row_power:
727: Mark the columns in row 'row' which are nonzero in
728: matrix^2log(marker).
729: */
730: PetscErrorCode spbas_mark_row_power(PetscInt *iwork, /* marker-vector */
731: PetscInt row, /* row for which the columns are marked */
732: spbas_matrix * in_matrix, /* matrix for which the power is being calculated */
733: PetscInt marker, /* marker-value: 2^power */
734: PetscInt minmrk, /* lower bound for marked points */
735: PetscInt maxmrk) /* upper bound for marked points */
736: {
737: PetscInt i,j, nnz;
739: nnz = in_matrix->row_nnz[row];
741: /* For higher powers, call this function recursively */
742: if (marker>1) {
743: for (i=0; i<nnz; i++) {
744: j = row + in_matrix->icols[row][i];
745: if (minmrk<=j && j<maxmrk && iwork[j] < marker) {
746: spbas_mark_row_power(iwork, row + in_matrix->icols[row][i],in_matrix, marker/2,minmrk,maxmrk);
747: iwork[j] |= marker;
748: }
749: }
750: } else {
751: /* Mark the columns reached */
752: for (i=0; i<nnz; i++) {
753: j = row + in_matrix->icols[row][i];
754: if (minmrk<=j && j<maxmrk) iwork[j] |= 1;
755: }
756: }
757: return 0;
758: }
760: /*
761: spbas_power
762: Calculate sparseness patterns for incomplete Cholesky decompositions
763: of a given order: (almost) all nonzeros of the matrix^(order+1) which
764: are inside the band width are found and stored in the output sparseness
765: pattern.
766: */
767: PetscErrorCode spbas_power(spbas_matrix in_matrix,PetscInt power, spbas_matrix * result)
768: {
769: spbas_matrix retval;
770: PetscInt nrows = in_matrix.nrows;
771: PetscInt ncols = in_matrix.ncols;
772: PetscInt i, j, kend;
773: PetscInt nnz, inz;
774: PetscInt *iwork;
775: PetscInt marker;
776: PetscInt maxmrk=0;
783: /* Copy input values*/
784: retval.nrows = ncols;
785: retval.ncols = nrows;
786: retval.nnz = 0;
787: retval.col_idx_type = SPBAS_DIAGONAL_OFFSETS;
788: retval.block_data = PETSC_FALSE;
790: /* Allocate sparseness pattern */
791: spbas_allocate_pattern(&retval, in_matrix.values ? PETSC_TRUE : PETSC_FALSE);
793: /* Allocate marker array: note sure the max needed so use the max of the two */
794: PetscCalloc1(PetscMax(ncols,nrows), &iwork);
796: /* Calculate marker values */
797: marker = 1; for (i=1; i<power; i++) marker*=2;
799: for (i=0; i<nrows; i++) {
800: /* Calculate the pattern for each row */
802: nnz = in_matrix.row_nnz[i];
803: kend = i+in_matrix.icols[i][nnz-1];
804: if (maxmrk<=kend) maxmrk=kend+1;
805: spbas_mark_row_power(iwork, i, &in_matrix, marker, i, maxmrk);
807: /* Count the columns*/
808: nnz = 0;
809: for (j=i; j<maxmrk; j++) nnz+= (iwork[j]!=0);
811: /* Allocate the column indices */
812: retval.row_nnz[i] = nnz;
813: PetscMalloc1(nnz,&retval.icols[i]);
815: /* Administrate the column indices */
816: inz = 0;
817: for (j=i; j<maxmrk; j++) {
818: if (iwork[j]) {
819: retval.icols[i][inz] = j-i;
820: inz++;
821: iwork[j] = 0;
822: }
823: }
824: retval.nnz += nnz;
825: };
826: PetscFree(iwork);
827: *result = retval;
828: return 0;
829: }
831: /*
832: spbas_keep_upper:
833: remove the lower part of the matrix: keep the upper part
834: */
835: PetscErrorCode spbas_keep_upper(spbas_matrix * inout_matrix)
836: {
837: PetscInt i, j;
838: PetscInt jstart;
841: for (i=0; i<inout_matrix->nrows; i++) {
842: for (jstart=0; (jstart<inout_matrix->row_nnz[i]) && (inout_matrix->icols[i][jstart]<0); jstart++) {}
843: if (jstart>0) {
844: for (j=0; j<inout_matrix->row_nnz[i]-jstart; j++) {
845: inout_matrix->icols[i][j] = inout_matrix->icols[i][j+jstart];
846: }
848: if (inout_matrix->values) {
849: for (j=0; j<inout_matrix->row_nnz[i]-jstart; j++) {
850: inout_matrix->values[i][j] = inout_matrix->values[i][j+jstart];
851: }
852: }
854: inout_matrix->row_nnz[i] -= jstart;
856: inout_matrix->icols[i] = (PetscInt*) realloc((void*) inout_matrix->icols[i], inout_matrix->row_nnz[i]*sizeof(PetscInt));
858: if (inout_matrix->values) {
859: inout_matrix->values[i] = (PetscScalar*) realloc((void*) inout_matrix->values[i], inout_matrix->row_nnz[i]*sizeof(PetscScalar));
860: }
861: inout_matrix->nnz -= jstart;
862: }
863: }
864: return 0;
865: }