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: Note:
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()`
23: M*/
25: /*
26: spbas_memory_requirement:
27: Calculate the number of bytes needed to store tha matrix
28: */
29: size_t spbas_memory_requirement(spbas_matrix matrix)
30: {
31: size_t memreq = 6 * sizeof(PetscInt) + /* nrows, ncols, nnz, n_alloc_icol, n_alloc_val, col_idx_type */
32: sizeof(PetscBool) + /* block_data */
33: sizeof(PetscScalar **) + /* values */
34: sizeof(PetscScalar *) + /* alloc_val */
35: 2 * sizeof(PetscInt **) + /* icols, icols0 */
36: 2 * sizeof(PetscInt *) + /* row_nnz, alloc_icol */
37: matrix.nrows * sizeof(PetscInt) + /* row_nnz[*] */
38: matrix.nrows * sizeof(PetscInt *); /* icols[*] */
40: /* icol0[*] */
41: if (matrix.col_idx_type == SPBAS_OFFSET_ARRAY) memreq += matrix.nrows * sizeof(PetscInt);
43: /* icols[*][*] */
44: if (matrix.block_data) memreq += matrix.n_alloc_icol * sizeof(PetscInt);
45: else memreq += matrix.nnz * sizeof(PetscInt);
47: if (matrix.values) {
48: memreq += matrix.nrows * sizeof(PetscScalar *); /* values[*] */
49: /* values[*][*] */
50: if (matrix.block_data) memreq += matrix.n_alloc_val * sizeof(PetscScalar);
51: else memreq += matrix.nnz * sizeof(PetscScalar);
52: }
53: return memreq;
54: }
56: /*
57: spbas_allocate_pattern:
58: allocate the pattern arrays row_nnz, icols and optionally values
59: */
60: PetscErrorCode spbas_allocate_pattern(spbas_matrix *result, PetscBool do_values)
61: {
62: PetscInt nrows = result->nrows;
63: PetscInt col_idx_type = result->col_idx_type;
65: /* Allocate sparseness pattern */
66: PetscMalloc1(nrows, &result->row_nnz);
67: PetscMalloc1(nrows, &result->icols);
69: /* If offsets are given wrt an array, create array */
70: if (col_idx_type == SPBAS_OFFSET_ARRAY) {
71: PetscMalloc1(nrows, &result->icol0);
72: } else {
73: result->icol0 = NULL;
74: }
76: /* If values are given, allocate values array */
77: if (do_values) {
78: PetscMalloc1(nrows, &result->values);
79: } else {
80: result->values = NULL;
81: }
82: return 0;
83: }
85: /*
86: spbas_allocate_data:
87: in case of block_data:
88: Allocate the data arrays alloc_icol and optionally alloc_val,
89: set appropriate pointers from icols and values;
90: in case of !block_data:
91: Allocate the arrays icols[i] and optionally values[i]
92: */
93: PetscErrorCode spbas_allocate_data(spbas_matrix *result)
94: {
95: PetscInt i;
96: PetscInt nnz = result->nnz;
97: PetscInt nrows = result->nrows;
98: PetscInt r_nnz;
99: PetscBool do_values = (result->values) ? PETSC_TRUE : PETSC_FALSE;
100: PetscBool block_data = result->block_data;
102: if (block_data) {
103: /* Allocate the column number array and point to it */
104: result->n_alloc_icol = nnz;
106: PetscMalloc1(nnz, &result->alloc_icol);
108: result->icols[0] = result->alloc_icol;
109: for (i = 1; i < nrows; i++) result->icols[i] = result->icols[i - 1] + result->row_nnz[i - 1];
111: /* Allocate the value array and point to it */
112: if (do_values) {
113: result->n_alloc_val = nnz;
115: PetscMalloc1(nnz, &result->alloc_val);
117: result->values[0] = result->alloc_val;
118: for (i = 1; i < nrows; i++) result->values[i] = result->values[i - 1] + result->row_nnz[i - 1];
119: }
120: } else {
121: for (i = 0; i < nrows; i++) {
122: r_nnz = result->row_nnz[i];
123: PetscMalloc1(r_nnz, &result->icols[i]);
124: }
125: if (do_values) {
126: for (i = 0; i < nrows; i++) {
127: r_nnz = result->row_nnz[i];
128: PetscMalloc1(r_nnz, &result->values[i]);
129: }
130: }
131: }
132: return 0;
133: }
135: /*
136: spbas_row_order_icol
137: determine if row i1 should come
138: + before row i2 in the sorted rows (return -1),
139: + after (return 1)
140: + is identical (return 0).
141: */
142: int spbas_row_order_icol(PetscInt i1, PetscInt i2, PetscInt *irow_in, PetscInt *icol_in, PetscInt col_idx_type)
143: {
144: PetscInt j;
145: PetscInt nnz1 = irow_in[i1 + 1] - irow_in[i1];
146: PetscInt nnz2 = irow_in[i2 + 1] - irow_in[i2];
147: PetscInt *icol1 = &icol_in[irow_in[i1]];
148: PetscInt *icol2 = &icol_in[irow_in[i2]];
150: if (nnz1 < nnz2) return -1;
151: if (nnz1 > nnz2) return 1;
153: if (col_idx_type == SPBAS_COLUMN_NUMBERS) {
154: for (j = 0; j < nnz1; j++) {
155: if (icol1[j] < icol2[j]) return -1;
156: if (icol1[j] > icol2[j]) return 1;
157: }
158: } else if (col_idx_type == SPBAS_DIAGONAL_OFFSETS) {
159: for (j = 0; j < nnz1; j++) {
160: if (icol1[j] - i1 < icol2[j] - i2) return -1;
161: if (icol1[j] - i1 > icol2[j] - i2) return 1;
162: }
163: } else if (col_idx_type == SPBAS_OFFSET_ARRAY) {
164: for (j = 1; j < nnz1; j++) {
165: if (icol1[j] - icol1[0] < icol2[j] - icol2[0]) return -1;
166: if (icol1[j] - icol1[0] > icol2[j] - icol2[0]) return 1;
167: }
168: }
169: return 0;
170: }
172: /*
173: spbas_mergesort_icols:
174: return a sorting of the rows in which identical sparseness patterns are
175: next to each other
176: */
177: PetscErrorCode spbas_mergesort_icols(PetscInt nrows, PetscInt *irow_in, PetscInt *icol_in, PetscInt col_idx_type, PetscInt *isort)
178: {
179: PetscInt istep; /* Chunk-sizes of already sorted parts of arrays */
180: PetscInt i, i1, i2; /* Loop counters for (partly) sorted arrays */
181: PetscInt istart, i1end, i2end; /* start of newly sorted array part, end of both parts */
182: PetscInt *ialloc; /* Allocated arrays */
183: PetscInt *iswap; /* auxiliary pointers for swapping */
184: PetscInt *ihlp1; /* Pointers to new version of arrays, */
185: PetscInt *ihlp2; /* Pointers to previous version of arrays, */
187: PetscMalloc1(nrows, &ialloc);
189: ihlp1 = ialloc;
190: ihlp2 = isort;
192: /* Sorted array chunks are first 1 long, and increase until they are the complete array */
193: for (istep = 1; istep < nrows; istep *= 2) {
194: /*
195: Combine sorted parts
196: istart:istart+istep-1 and istart+istep-1:istart+2*istep-1
197: of ihlp2 and vhlp2
199: into one sorted part
200: istart:istart+2*istep-1
201: of ihlp1 and vhlp1
202: */
203: for (istart = 0; istart < nrows; istart += 2 * istep) {
204: /* Set counters and bound array part endings */
205: i1 = istart;
206: i1end = i1 + istep;
207: if (i1end > nrows) i1end = nrows;
208: i2 = istart + istep;
209: i2end = i2 + istep;
210: if (i2end > nrows) i2end = nrows;
212: /* Merge the two array parts */
213: for (i = istart; i < i2end; i++) {
214: if (i1 < i1end && i2 < i2end && spbas_row_order_icol(ihlp2[i1], ihlp2[i2], irow_in, icol_in, col_idx_type) < 0) {
215: ihlp1[i] = ihlp2[i1];
216: i1++;
217: } else if (i2 < i2end) {
218: ihlp1[i] = ihlp2[i2];
219: i2++;
220: } else {
221: ihlp1[i] = ihlp2[i1];
222: i1++;
223: }
224: }
225: }
227: /* Swap the two array sets */
228: iswap = ihlp2;
229: ihlp2 = ihlp1;
230: ihlp1 = iswap;
231: }
233: /* Copy one more time in case the sorted arrays are the temporary ones */
234: if (ihlp2 != isort) {
235: for (i = 0; i < nrows; i++) isort[i] = ihlp2[i];
236: }
237: PetscFree(ialloc);
238: return 0;
239: }
241: /*
242: spbas_compress_pattern:
243: calculate a compressed sparseness pattern for a sparseness pattern
244: given in compressed row storage. The compressed sparseness pattern may
245: require (much) less memory.
246: */
247: PetscErrorCode spbas_compress_pattern(PetscInt *irow_in, PetscInt *icol_in, PetscInt nrows, PetscInt ncols, PetscInt col_idx_type, spbas_matrix *B, PetscReal *mem_reduction)
248: {
249: PetscInt nnz = irow_in[nrows];
250: size_t mem_orig = (nrows + nnz) * sizeof(PetscInt);
251: size_t mem_compressed;
252: PetscInt *isort;
253: PetscInt *icols;
254: PetscInt row_nnz;
255: PetscInt *ipoint;
256: PetscBool *used;
257: PetscInt ptr;
258: PetscInt i, j;
259: const PetscBool no_values = PETSC_FALSE;
261: /* Allocate the structure of the new matrix */
262: B->nrows = nrows;
263: B->ncols = ncols;
264: B->nnz = nnz;
265: B->col_idx_type = col_idx_type;
266: B->block_data = PETSC_TRUE;
268: spbas_allocate_pattern(B, no_values);
270: /* When using an offset array, set it */
271: if (col_idx_type == SPBAS_OFFSET_ARRAY) {
272: for (i = 0; i < nrows; i++) B->icol0[i] = icol_in[irow_in[i]];
273: }
275: /* Allocate the ordering for the rows */
276: PetscMalloc1(nrows, &isort);
277: PetscMalloc1(nrows, &ipoint);
278: PetscCalloc1(nrows, &used);
280: for (i = 0; i < nrows; i++) {
281: B->row_nnz[i] = irow_in[i + 1] - irow_in[i];
282: isort[i] = i;
283: ipoint[i] = i;
284: }
286: /* Sort the rows so that identical columns will be next to each other */
287: spbas_mergesort_icols(nrows, irow_in, icol_in, col_idx_type, isort);
288: PetscInfo(NULL, "Rows have been sorted for patterns\n");
290: /* Replace identical rows with the first one in the list */
291: for (i = 1; i < nrows; i++) {
292: if (spbas_row_order_icol(isort[i - 1], isort[i], irow_in, icol_in, col_idx_type) == 0) ipoint[isort[i]] = ipoint[isort[i - 1]];
293: }
295: /* Collect the rows which are used*/
296: for (i = 0; i < nrows; i++) used[ipoint[i]] = PETSC_TRUE;
298: /* Calculate needed memory */
299: B->n_alloc_icol = 0;
300: for (i = 0; i < nrows; i++) {
301: if (used[i]) B->n_alloc_icol += B->row_nnz[i];
302: }
303: PetscMalloc1(B->n_alloc_icol, &B->alloc_icol);
305: /* Fill in the diagonal offsets for the rows which store their own data */
306: ptr = 0;
307: for (i = 0; i < B->nrows; i++) {
308: if (used[i]) {
309: B->icols[i] = &B->alloc_icol[ptr];
310: icols = &icol_in[irow_in[i]];
311: row_nnz = B->row_nnz[i];
312: if (col_idx_type == SPBAS_COLUMN_NUMBERS) {
313: for (j = 0; j < row_nnz; j++) B->icols[i][j] = icols[j];
314: } else if (col_idx_type == SPBAS_DIAGONAL_OFFSETS) {
315: for (j = 0; j < row_nnz; j++) B->icols[i][j] = icols[j] - i;
316: } else if (col_idx_type == SPBAS_OFFSET_ARRAY) {
317: for (j = 0; j < row_nnz; j++) B->icols[i][j] = icols[j] - icols[0];
318: }
319: ptr += B->row_nnz[i];
320: }
321: }
323: /* Point to the right places for all data */
324: for (i = 0; i < nrows; i++) B->icols[i] = B->icols[ipoint[i]];
325: PetscInfo(NULL, "Row patterns have been compressed\n");
326: PetscInfo(NULL, " (%g nonzeros per row)\n", (double)((PetscReal)nnz / (PetscReal)nrows));
328: PetscFree(isort);
329: PetscFree(used);
330: PetscFree(ipoint);
332: mem_compressed = spbas_memory_requirement(*B);
333: *mem_reduction = 100.0 * (PetscReal)(mem_orig - mem_compressed) / (PetscReal)mem_orig;
334: return 0;
335: }
337: /*
338: spbas_incomplete_cholesky
339: Incomplete Cholesky decomposition
340: */
341: #include <../src/mat/impls/aij/seq/bas/spbas_cholesky.h>
343: /*
344: spbas_delete : de-allocate the arrays owned by this matrix
345: */
346: PetscErrorCode spbas_delete(spbas_matrix matrix)
347: {
348: PetscInt i;
350: if (matrix.block_data) {
351: PetscFree(matrix.alloc_icol);
352: if (matrix.values) PetscFree(matrix.alloc_val);
353: } else {
354: for (i = 0; i < matrix.nrows; i++) PetscFree(matrix.icols[i]);
355: PetscFree(matrix.icols);
356: if (matrix.values) {
357: for (i = 0; i < matrix.nrows; i++) PetscFree(matrix.values[i]);
358: }
359: }
361: PetscFree(matrix.row_nnz);
362: PetscFree(matrix.icols);
363: if (matrix.col_idx_type == SPBAS_OFFSET_ARRAY) PetscFree(matrix.icol0);
364: PetscFree(matrix.values);
365: return 0;
366: }
368: /*
369: spbas_matrix_to_crs:
370: Convert an spbas_matrix to compessed row storage
371: */
372: PetscErrorCode spbas_matrix_to_crs(spbas_matrix matrix_A, MatScalar **val_out, PetscInt **irow_out, PetscInt **icol_out)
373: {
374: PetscInt nrows = matrix_A.nrows;
375: PetscInt nnz = matrix_A.nnz;
376: PetscInt i, j, r_nnz, i0;
377: PetscInt *irow;
378: PetscInt *icol;
379: PetscInt *icol_A;
380: MatScalar *val;
381: PetscScalar *val_A;
382: PetscInt col_idx_type = matrix_A.col_idx_type;
383: PetscBool do_values = matrix_A.values ? PETSC_TRUE : PETSC_FALSE;
385: PetscMalloc1(nrows + 1, &irow);
386: PetscMalloc1(nnz, &icol);
387: *icol_out = icol;
388: *irow_out = irow;
389: if (do_values) {
390: PetscMalloc1(nnz, &val);
391: *val_out = val;
392: *icol_out = icol;
393: *irow_out = irow;
394: }
396: irow[0] = 0;
397: for (i = 0; i < nrows; i++) {
398: r_nnz = matrix_A.row_nnz[i];
399: i0 = irow[i];
400: irow[i + 1] = i0 + r_nnz;
401: icol_A = matrix_A.icols[i];
403: if (do_values) {
404: val_A = matrix_A.values[i];
405: for (j = 0; j < r_nnz; j++) {
406: icol[i0 + j] = icol_A[j];
407: val[i0 + j] = val_A[j];
408: }
409: } else {
410: for (j = 0; j < r_nnz; j++) icol[i0 + j] = icol_A[j];
411: }
413: if (col_idx_type == SPBAS_DIAGONAL_OFFSETS) {
414: for (j = 0; j < r_nnz; j++) icol[i0 + j] += i;
415: } else if (col_idx_type == SPBAS_OFFSET_ARRAY) {
416: i0 = matrix_A.icol0[i];
417: for (j = 0; j < r_nnz; j++) icol[i0 + j] += i0;
418: }
419: }
420: return 0;
421: }
423: /*
424: spbas_transpose
425: return the transpose of a matrix
426: */
427: PetscErrorCode spbas_transpose(spbas_matrix in_matrix, spbas_matrix *result)
428: {
429: PetscInt col_idx_type = in_matrix.col_idx_type;
430: PetscInt nnz = in_matrix.nnz;
431: PetscInt ncols = in_matrix.nrows;
432: PetscInt nrows = in_matrix.ncols;
433: PetscInt i, j, k;
434: PetscInt r_nnz;
435: PetscInt *irow;
436: PetscInt icol0 = 0;
437: PetscScalar *val;
439: /* Copy input values */
440: result->nrows = nrows;
441: result->ncols = ncols;
442: result->nnz = nnz;
443: result->col_idx_type = SPBAS_COLUMN_NUMBERS;
444: result->block_data = PETSC_TRUE;
446: /* Allocate sparseness pattern */
447: spbas_allocate_pattern(result, in_matrix.values ? PETSC_TRUE : PETSC_FALSE);
449: /* Count the number of nonzeros in each row */
450: for (i = 0; i < nrows; i++) result->row_nnz[i] = 0;
452: for (i = 0; i < ncols; i++) {
453: r_nnz = in_matrix.row_nnz[i];
454: irow = in_matrix.icols[i];
455: if (col_idx_type == SPBAS_COLUMN_NUMBERS) {
456: for (j = 0; j < r_nnz; j++) result->row_nnz[irow[j]]++;
457: } else if (col_idx_type == SPBAS_DIAGONAL_OFFSETS) {
458: for (j = 0; j < r_nnz; j++) result->row_nnz[i + irow[j]]++;
459: } else if (col_idx_type == SPBAS_OFFSET_ARRAY) {
460: icol0 = in_matrix.icol0[i];
461: for (j = 0; j < r_nnz; j++) result->row_nnz[icol0 + irow[j]]++;
462: }
463: }
465: /* Set the pointers to the data */
466: spbas_allocate_data(result);
468: /* Reset the number of nonzeros in each row */
469: for (i = 0; i < nrows; i++) result->row_nnz[i] = 0;
471: /* Fill the data arrays */
472: if (in_matrix.values) {
473: for (i = 0; i < ncols; i++) {
474: r_nnz = in_matrix.row_nnz[i];
475: irow = in_matrix.icols[i];
476: val = in_matrix.values[i];
478: if (col_idx_type == SPBAS_COLUMN_NUMBERS) icol0 = 0;
479: else if (col_idx_type == SPBAS_DIAGONAL_OFFSETS) icol0 = i;
480: else if (col_idx_type == SPBAS_OFFSET_ARRAY) icol0 = in_matrix.icol0[i];
481: for (j = 0; j < r_nnz; j++) {
482: k = icol0 + irow[j];
483: result->icols[k][result->row_nnz[k]] = i;
484: result->values[k][result->row_nnz[k]] = val[j];
485: result->row_nnz[k]++;
486: }
487: }
488: } else {
489: for (i = 0; i < ncols; i++) {
490: r_nnz = in_matrix.row_nnz[i];
491: irow = in_matrix.icols[i];
493: if (col_idx_type == SPBAS_COLUMN_NUMBERS) icol0 = 0;
494: else if (col_idx_type == SPBAS_DIAGONAL_OFFSETS) icol0 = i;
495: else if (col_idx_type == SPBAS_OFFSET_ARRAY) icol0 = in_matrix.icol0[i];
497: for (j = 0; j < r_nnz; j++) {
498: k = icol0 + irow[j];
499: result->icols[k][result->row_nnz[k]] = i;
500: result->row_nnz[k]++;
501: }
502: }
503: }
504: return 0;
505: }
507: /*
508: spbas_mergesort
510: mergesort for an array of integers and an array of associated
511: reals
513: on output, icol[0..nnz-1] is increasing;
514: val[0..nnz-1] has undergone the same permutation as icol
516: NB: val may be NULL: in that case, only the integers are sorted
518: */
519: PetscErrorCode spbas_mergesort(PetscInt nnz, PetscInt *icol, PetscScalar *val)
520: {
521: PetscInt istep; /* Chunk-sizes of already sorted parts of arrays */
522: PetscInt i, i1, i2; /* Loop counters for (partly) sorted arrays */
523: PetscInt istart, i1end, i2end; /* start of newly sorted array part, end of both parts */
524: PetscInt *ialloc; /* Allocated arrays */
525: PetscScalar *valloc = NULL;
526: PetscInt *iswap; /* auxiliary pointers for swapping */
527: PetscScalar *vswap;
528: PetscInt *ihlp1; /* Pointers to new version of arrays, */
529: PetscScalar *vhlp1 = NULL; /* (arrays under construction) */
530: PetscInt *ihlp2; /* Pointers to previous version of arrays, */
531: PetscScalar *vhlp2 = NULL;
533: PetscMalloc1(nnz, &ialloc);
534: ihlp1 = ialloc;
535: ihlp2 = icol;
537: if (val) {
538: PetscMalloc1(nnz, &valloc);
539: vhlp1 = valloc;
540: vhlp2 = val;
541: }
543: /* Sorted array chunks are first 1 long, and increase until they are the complete array */
544: for (istep = 1; istep < nnz; istep *= 2) {
545: /*
546: Combine sorted parts
547: istart:istart+istep-1 and istart+istep-1:istart+2*istep-1
548: of ihlp2 and vhlp2
550: into one sorted part
551: istart:istart+2*istep-1
552: of ihlp1 and vhlp1
553: */
554: for (istart = 0; istart < nnz; istart += 2 * istep) {
555: /* Set counters and bound array part endings */
556: i1 = istart;
557: i1end = i1 + istep;
558: if (i1end > nnz) i1end = nnz;
559: i2 = istart + istep;
560: i2end = i2 + istep;
561: if (i2end > nnz) i2end = nnz;
563: /* Merge the two array parts */
564: if (val) {
565: for (i = istart; i < i2end; i++) {
566: if (i1 < i1end && i2 < i2end && ihlp2[i1] < ihlp2[i2]) {
567: ihlp1[i] = ihlp2[i1];
568: vhlp1[i] = vhlp2[i1];
569: i1++;
570: } else if (i2 < i2end) {
571: ihlp1[i] = ihlp2[i2];
572: vhlp1[i] = vhlp2[i2];
573: i2++;
574: } else {
575: ihlp1[i] = ihlp2[i1];
576: vhlp1[i] = vhlp2[i1];
577: i1++;
578: }
579: }
580: } else {
581: for (i = istart; i < i2end; i++) {
582: if (i1 < i1end && i2 < i2end && ihlp2[i1] < ihlp2[i2]) {
583: ihlp1[i] = ihlp2[i1];
584: i1++;
585: } else if (i2 < i2end) {
586: ihlp1[i] = ihlp2[i2];
587: i2++;
588: } else {
589: ihlp1[i] = ihlp2[i1];
590: i1++;
591: }
592: }
593: }
594: }
596: /* Swap the two array sets */
597: iswap = ihlp2;
598: ihlp2 = ihlp1;
599: ihlp1 = iswap;
600: vswap = vhlp2;
601: vhlp2 = vhlp1;
602: vhlp1 = vswap;
603: }
605: /* Copy one more time in case the sorted arrays are the temporary ones */
606: if (ihlp2 != icol) {
607: for (i = 0; i < nnz; i++) icol[i] = ihlp2[i];
608: if (val) {
609: for (i = 0; i < nnz; i++) val[i] = vhlp2[i];
610: }
611: }
613: PetscFree(ialloc);
614: if (val) PetscFree(valloc);
615: return 0;
616: }
618: /*
619: spbas_apply_reordering_rows:
620: apply the given reordering to the rows: matrix_A = matrix_A(perm,:);
621: */
622: PetscErrorCode spbas_apply_reordering_rows(spbas_matrix *matrix_A, const PetscInt *permutation)
623: {
624: PetscInt i, j, ip;
625: PetscInt nrows = matrix_A->nrows;
626: PetscInt *row_nnz;
627: PetscInt **icols;
628: PetscBool do_values = matrix_A->values ? PETSC_TRUE : PETSC_FALSE;
629: PetscScalar **vals = NULL;
633: if (do_values) PetscMalloc1(nrows, &vals);
634: PetscMalloc1(nrows, &row_nnz);
635: PetscMalloc1(nrows, &icols);
637: for (i = 0; i < nrows; i++) {
638: ip = permutation[i];
639: if (do_values) vals[i] = matrix_A->values[ip];
640: icols[i] = matrix_A->icols[ip];
641: row_nnz[i] = matrix_A->row_nnz[ip];
642: for (j = 0; j < row_nnz[i]; j++) icols[i][j] += ip - i;
643: }
645: if (do_values) PetscFree(matrix_A->values);
646: PetscFree(matrix_A->icols);
647: PetscFree(matrix_A->row_nnz);
649: if (do_values) matrix_A->values = vals;
650: matrix_A->icols = icols;
651: matrix_A->row_nnz = row_nnz;
652: return 0;
653: }
655: /*
656: spbas_apply_reordering_cols:
657: apply the given reordering to the columns: matrix_A(:,perm) = matrix_A;
658: */
659: PetscErrorCode spbas_apply_reordering_cols(spbas_matrix *matrix_A, const PetscInt *permutation)
660: {
661: PetscInt i, j;
662: PetscInt nrows = matrix_A->nrows;
663: PetscInt row_nnz;
664: PetscInt *icols;
665: PetscBool do_values = matrix_A->values ? PETSC_TRUE : PETSC_FALSE;
666: PetscScalar *vals = NULL;
670: for (i = 0; i < nrows; i++) {
671: icols = matrix_A->icols[i];
672: row_nnz = matrix_A->row_nnz[i];
673: if (do_values) vals = matrix_A->values[i];
675: for (j = 0; j < row_nnz; j++) icols[j] = permutation[i + icols[j]] - i;
676: spbas_mergesort(row_nnz, icols, vals);
677: }
678: return 0;
679: }
681: /*
682: spbas_apply_reordering:
683: apply the given reordering: matrix_A(perm,perm) = matrix_A;
684: */
685: PetscErrorCode spbas_apply_reordering(spbas_matrix *matrix_A, const PetscInt *permutation, const PetscInt *inv_perm)
686: {
687: spbas_apply_reordering_rows(matrix_A, inv_perm);
688: spbas_apply_reordering_cols(matrix_A, permutation);
689: return 0;
690: }
692: PetscErrorCode spbas_pattern_only(PetscInt nrows, PetscInt ncols, PetscInt *ai, PetscInt *aj, spbas_matrix *result)
693: {
694: spbas_matrix retval;
695: PetscInt i, j, i0, r_nnz;
697: /* Copy input values */
698: retval.nrows = nrows;
699: retval.ncols = ncols;
700: retval.nnz = ai[nrows];
702: retval.block_data = PETSC_TRUE;
703: retval.col_idx_type = SPBAS_DIAGONAL_OFFSETS;
705: /* Allocate output matrix */
706: spbas_allocate_pattern(&retval, PETSC_FALSE);
707: for (i = 0; i < nrows; i++) retval.row_nnz[i] = ai[i + 1] - ai[i];
708: spbas_allocate_data(&retval);
709: /* Copy the structure */
710: for (i = 0; i < retval.nrows; i++) {
711: i0 = ai[i];
712: r_nnz = ai[i + 1] - i0;
714: for (j = 0; j < r_nnz; j++) retval.icols[i][j] = aj[i0 + j] - i;
715: }
716: *result = retval;
717: return 0;
718: }
720: /*
721: spbas_mark_row_power:
722: Mark the columns in row 'row' which are nonzero in
723: matrix^2log(marker).
724: */
725: PetscErrorCode spbas_mark_row_power(PetscInt *iwork, /* marker-vector */
726: PetscInt row, /* row for which the columns are marked */
727: spbas_matrix *in_matrix, /* matrix for which the power is being calculated */
728: PetscInt marker, /* marker-value: 2^power */
729: PetscInt minmrk, /* lower bound for marked points */
730: PetscInt maxmrk) /* upper bound for marked points */
731: {
732: PetscInt i, j, nnz;
734: nnz = in_matrix->row_nnz[row];
736: /* For higher powers, call this function recursively */
737: if (marker > 1) {
738: for (i = 0; i < nnz; i++) {
739: j = row + in_matrix->icols[row][i];
740: if (minmrk <= j && j < maxmrk && iwork[j] < marker) {
741: spbas_mark_row_power(iwork, row + in_matrix->icols[row][i], in_matrix, marker / 2, minmrk, maxmrk);
742: iwork[j] |= marker;
743: }
744: }
745: } else {
746: /* Mark the columns reached */
747: for (i = 0; i < nnz; i++) {
748: j = row + in_matrix->icols[row][i];
749: if (minmrk <= j && j < maxmrk) iwork[j] |= 1;
750: }
751: }
752: return 0;
753: }
755: /*
756: spbas_power
757: Calculate sparseness patterns for incomplete Cholesky decompositions
758: of a given order: (almost) all nonzeros of the matrix^(order+1) which
759: are inside the band width are found and stored in the output sparseness
760: pattern.
761: */
762: PetscErrorCode spbas_power(spbas_matrix in_matrix, PetscInt power, spbas_matrix *result)
763: {
764: spbas_matrix retval;
765: PetscInt nrows = in_matrix.nrows;
766: PetscInt ncols = in_matrix.ncols;
767: PetscInt i, j, kend;
768: PetscInt nnz, inz;
769: PetscInt *iwork;
770: PetscInt marker;
771: PetscInt maxmrk = 0;
778: /* Copy input values*/
779: retval.nrows = ncols;
780: retval.ncols = nrows;
781: retval.nnz = 0;
782: retval.col_idx_type = SPBAS_DIAGONAL_OFFSETS;
783: retval.block_data = PETSC_FALSE;
785: /* Allocate sparseness pattern */
786: spbas_allocate_pattern(&retval, in_matrix.values ? PETSC_TRUE : PETSC_FALSE);
788: /* Allocate marker array: note sure the max needed so use the max of the two */
789: PetscCalloc1(PetscMax(ncols, nrows), &iwork);
791: /* Calculate marker values */
792: marker = 1;
793: for (i = 1; i < power; i++) marker *= 2;
795: for (i = 0; i < nrows; i++) {
796: /* Calculate the pattern for each row */
798: nnz = in_matrix.row_nnz[i];
799: kend = i + in_matrix.icols[i][nnz - 1];
800: if (maxmrk <= kend) maxmrk = kend + 1;
801: spbas_mark_row_power(iwork, i, &in_matrix, marker, i, maxmrk);
803: /* Count the columns*/
804: nnz = 0;
805: for (j = i; j < maxmrk; j++) nnz += (iwork[j] != 0);
807: /* Allocate the column indices */
808: retval.row_nnz[i] = nnz;
809: PetscMalloc1(nnz, &retval.icols[i]);
811: /* Administrate the column indices */
812: inz = 0;
813: for (j = i; j < maxmrk; j++) {
814: if (iwork[j]) {
815: retval.icols[i][inz] = j - i;
816: inz++;
817: iwork[j] = 0;
818: }
819: }
820: retval.nnz += nnz;
821: };
822: PetscFree(iwork);
823: *result = retval;
824: return 0;
825: }
827: /*
828: spbas_keep_upper:
829: remove the lower part of the matrix: keep the upper part
830: */
831: PetscErrorCode spbas_keep_upper(spbas_matrix *inout_matrix)
832: {
833: PetscInt i, j;
834: PetscInt jstart;
837: for (i = 0; i < inout_matrix->nrows; i++) {
838: for (jstart = 0; (jstart < inout_matrix->row_nnz[i]) && (inout_matrix->icols[i][jstart] < 0); jstart++) { }
839: if (jstart > 0) {
840: for (j = 0; j < inout_matrix->row_nnz[i] - jstart; j++) inout_matrix->icols[i][j] = inout_matrix->icols[i][j + jstart];
842: if (inout_matrix->values) {
843: for (j = 0; j < inout_matrix->row_nnz[i] - jstart; j++) inout_matrix->values[i][j] = inout_matrix->values[i][j + jstart];
844: }
846: inout_matrix->row_nnz[i] -= jstart;
848: inout_matrix->icols[i] = (PetscInt *)realloc((void *)inout_matrix->icols[i], inout_matrix->row_nnz[i] * sizeof(PetscInt));
850: if (inout_matrix->values) inout_matrix->values[i] = (PetscScalar *)realloc((void *)inout_matrix->values[i], inout_matrix->row_nnz[i] * sizeof(PetscScalar));
851: inout_matrix->nnz -= jstart;
852: }
853: }
854: return 0;
855: }