Actual source code: spbas.c
petsc-3.9.4 2018-09-11
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: Since this currently hooked up to use drop tolerance it should produce the same factors and hence convergence as the PETSc ICC, for higher
18: 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
19: we recommend not using this functionality.
21: .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: {
63: PetscInt nrows = result->nrows;
64: PetscInt col_idx_type = result->col_idx_type;
67: /* Allocate sparseness pattern */
68: PetscMalloc1(nrows,&result->row_nnz);
69: PetscMalloc1(nrows,&result->icols);
71: /* If offsets are given wrt an array, create array */
72: if (col_idx_type == SPBAS_OFFSET_ARRAY) {
73: PetscMalloc1(nrows,&result->icol0);
74: } else {
75: result->icol0 = NULL;
76: }
78: /* If values are given, allocate values array */
79: if (do_values) {
80: PetscMalloc1(nrows,&result->values);
81: } else {
82: result->values = NULL;
83: }
84: return(0);
85: }
87: /*
88: spbas_allocate_data:
89: in case of block_data:
90: Allocate the data arrays alloc_icol and optionally alloc_val,
91: set appropriate pointers from icols and values;
92: in case of !block_data:
93: Allocate the arrays icols[i] and optionally values[i]
94: */
95: PetscErrorCode spbas_allocate_data(spbas_matrix * result)
96: {
97: PetscInt i;
98: PetscInt nnz = result->nnz;
99: PetscInt nrows = result->nrows;
100: PetscInt r_nnz;
102: PetscBool do_values = (result->values) ? PETSC_TRUE : PETSC_FALSE;
103: PetscBool block_data = result->block_data;
106: if (block_data) {
107: /* Allocate the column number array and point to it */
108: result->n_alloc_icol = nnz;
110: PetscMalloc1(nnz, &result->alloc_icol);
112: result->icols[0] = result->alloc_icol;
113: for (i=1; i<nrows; i++) {
114: result->icols[i] = result->icols[i-1] + result->row_nnz[i-1];
115: }
117: /* Allocate the value array and point to it */
118: if (do_values) {
119: result->n_alloc_val = nnz;
121: PetscMalloc1(nnz, &result->alloc_val);
123: result->values[0] = result->alloc_val;
124: for (i=1; i<nrows; i++) {
125: result->values[i] = result->values[i-1] + result->row_nnz[i-1];
126: }
127: }
128: } else {
129: for (i=0; i<nrows; i++) {
130: r_nnz = result->row_nnz[i];
131: PetscMalloc1(r_nnz, &result->icols[i]);
132: }
133: if (do_values) {
134: for (i=0; i<nrows; i++) {
135: r_nnz = result->row_nnz[i];
136: PetscMalloc1(r_nnz, &result->values[i]);
137: }
138: }
139: }
140: return(0);
141: }
143: /*
144: spbas_row_order_icol
145: determine if row i1 should come
146: + before row i2 in the sorted rows (return -1),
147: + after (return 1)
148: + is identical (return 0).
149: */
150: int spbas_row_order_icol(PetscInt i1, PetscInt i2, PetscInt *irow_in, PetscInt *icol_in,PetscInt col_idx_type)
151: {
152: PetscInt j;
153: PetscInt nnz1 = irow_in[i1+1] - irow_in[i1];
154: PetscInt nnz2 = irow_in[i2+1] - irow_in[i2];
155: PetscInt * icol1 = &icol_in[irow_in[i1]];
156: PetscInt * icol2 = &icol_in[irow_in[i2]];
158: if (nnz1<nnz2) return -1;
159: if (nnz1>nnz2) return 1;
161: if (col_idx_type == SPBAS_COLUMN_NUMBERS) {
162: for (j=0; j<nnz1; j++) {
163: if (icol1[j]< icol2[j]) return -1;
164: if (icol1[j]> icol2[j]) return 1;
165: }
166: } else if (col_idx_type == SPBAS_DIAGONAL_OFFSETS) {
167: for (j=0; j<nnz1; j++) {
168: if (icol1[j]-i1< icol2[j]-i2) return -1;
169: if (icol1[j]-i1> icol2[j]-i2) return 1;
170: }
171: } else if (col_idx_type == SPBAS_OFFSET_ARRAY) {
172: for (j=1; j<nnz1; j++) {
173: if (icol1[j]-icol1[0] < icol2[j]-icol2[0]) return -1;
174: if (icol1[j]-icol1[0] > icol2[j]-icol2[0]) return 1;
175: }
176: }
177: return 0;
178: }
180: /*
181: spbas_mergesort_icols:
182: return a sorting of the rows in which identical sparseness patterns are
183: next to each other
184: */
185: PetscErrorCode spbas_mergesort_icols(PetscInt nrows, PetscInt * irow_in, PetscInt * icol_in,PetscInt col_idx_type, PetscInt *isort)
186: {
188: PetscInt istep; /* Chunk-sizes of already sorted parts of arrays */
189: PetscInt i, i1, i2; /* Loop counters for (partly) sorted arrays */
190: PetscInt istart, i1end, i2end; /* start of newly sorted array part, end of both parts */
191: PetscInt *ialloc; /* Allocated arrays */
192: PetscInt *iswap; /* auxiliary pointers for swapping */
193: PetscInt *ihlp1; /* Pointers to new version of arrays, */
194: PetscInt *ihlp2; /* Pointers to previous version of arrays, */
197: PetscMalloc1(nrows,&ialloc);
199: ihlp1 = ialloc;
200: ihlp2 = isort;
202: /* Sorted array chunks are first 1 long, and increase until they are the complete array */
203: for (istep=1; istep<nrows; istep*=2) {
204: /*
205: Combine sorted parts
206: istart:istart+istep-1 and istart+istep-1:istart+2*istep-1
207: of ihlp2 and vhlp2
209: into one sorted part
210: istart:istart+2*istep-1
211: of ihlp1 and vhlp1
212: */
213: for (istart=0; istart<nrows; istart+=2*istep) {
214: /* Set counters and bound array part endings */
215: i1=istart; i1end = i1+istep; if (i1end>nrows) i1end=nrows;
216: i2=istart+istep; i2end = i2+istep; if (i2end>nrows) i2end=nrows;
218: /* Merge the two array parts */
219: for (i=istart; i<i2end; i++) {
220: if (i1<i1end && i2<i2end && spbas_row_order_icol(ihlp2[i1], ihlp2[i2], irow_in, icol_in, col_idx_type) < 0) {
221: ihlp1[i] = ihlp2[i1];
222: i1++;
223: } else if (i2<i2end) {
224: ihlp1[i] = ihlp2[i2];
225: i2++;
226: } else {
227: ihlp1[i] = ihlp2[i1];
228: i1++;
229: }
230: }
231: }
233: /* Swap the two array sets */
234: iswap = ihlp2; ihlp2 = ihlp1; ihlp1 = iswap;
235: }
237: /* Copy one more time in case the sorted arrays are the temporary ones */
238: if (ihlp2 != isort) {
239: for (i=0; i<nrows; i++) isort[i] = ihlp2[i];
240: }
241: PetscFree(ialloc);
242: return(0);
243: }
247: /*
248: spbas_compress_pattern:
249: calculate a compressed sparseness pattern for a sparseness pattern
250: given in compressed row storage. The compressed sparseness pattern may
251: require (much) less memory.
252: */
253: PetscErrorCode spbas_compress_pattern(PetscInt *irow_in, PetscInt *icol_in, PetscInt nrows, PetscInt ncols, PetscInt col_idx_type, spbas_matrix *B,PetscReal *mem_reduction)
254: {
255: PetscInt nnz = irow_in[nrows];
256: size_t mem_orig = (nrows + nnz) * sizeof(PetscInt);
257: size_t mem_compressed;
258: PetscErrorCode ierr;
259: PetscInt *isort;
260: PetscInt *icols;
261: PetscInt row_nnz;
262: PetscInt *ipoint;
263: PetscBool *used;
264: PetscInt ptr;
265: PetscInt i,j;
266: const PetscBool no_values = PETSC_FALSE;
269: /* Allocate the structure of the new matrix */
270: B->nrows = nrows;
271: B->ncols = ncols;
272: B->nnz = nnz;
273: B->col_idx_type = col_idx_type;
274: B->block_data = PETSC_TRUE;
276: spbas_allocate_pattern(B, no_values);
278: /* When using an offset array, set it */
279: if (col_idx_type==SPBAS_OFFSET_ARRAY) {
280: for (i=0; i<nrows; i++) B->icol0[i] = icol_in[irow_in[i]];
281: }
283: /* Allocate the ordering for the rows */
284: PetscMalloc1(nrows,&isort);
285: PetscMalloc1(nrows,&ipoint);
286: PetscMalloc1(nrows,&used);
288: /* Initialize the sorting */
289: PetscMemzero((void*) used, nrows*sizeof(PetscBool));
290: for (i = 0; i<nrows; i++) {
291: B->row_nnz[i] = irow_in[i+1]-irow_in[i];
292: isort[i] = i;
293: ipoint[i] = i;
294: }
296: /* Sort the rows so that identical columns will be next to each other */
297: spbas_mergesort_icols(nrows, irow_in, icol_in, col_idx_type, isort);
298: PetscInfo(NULL,"Rows have been sorted for patterns\n");
300: /* Replace identical rows with the first one in the list */
301: for (i=1; i<nrows; i++) {
302: if (spbas_row_order_icol(isort[i-1], isort[i], irow_in, icol_in, col_idx_type) == 0) {
303: ipoint[isort[i]] = ipoint[isort[i-1]];
304: }
305: }
307: /* Collect the rows which are used*/
308: for (i=0; i<nrows; i++) used[ipoint[i]] = PETSC_TRUE;
310: /* Calculate needed memory */
311: B->n_alloc_icol = 0;
312: for (i=0; i<nrows; i++) {
313: if (used[i]) B->n_alloc_icol += B->row_nnz[i];
314: }
315: PetscMalloc1(B->n_alloc_icol,&B->alloc_icol);
317: /* Fill in the diagonal offsets for the rows which store their own data */
318: ptr = 0;
319: for (i=0; i<B->nrows; i++) {
320: if (used[i]) {
321: B->icols[i] = &B->alloc_icol[ptr];
322: icols = &icol_in[irow_in[i]];
323: row_nnz = B->row_nnz[i];
324: if (col_idx_type == SPBAS_COLUMN_NUMBERS) {
325: for (j=0; j<row_nnz; j++) {
326: B->icols[i][j] = icols[j];
327: }
328: } else if (col_idx_type == SPBAS_DIAGONAL_OFFSETS) {
329: for (j=0; j<row_nnz; j++) {
330: B->icols[i][j] = icols[j]-i;
331: }
332: } else if (col_idx_type == SPBAS_OFFSET_ARRAY) {
333: for (j=0; j<row_nnz; j++) {
334: B->icols[i][j] = icols[j]-icols[0];
335: }
336: }
337: ptr += B->row_nnz[i];
338: }
339: }
341: /* Point to the right places for all data */
342: for (i=0; i<nrows; i++) {
343: B->icols[i] = B->icols[ipoint[i]];
344: }
345: PetscInfo(NULL,"Row patterns have been compressed\n");
346: PetscInfo1(NULL," (%g nonzeros per row)\n", (double) ((PetscReal) nnz / (PetscReal) nrows));
348: ierr=PetscFree(isort);
349: ierr=PetscFree(used);
350: ierr=PetscFree(ipoint);
352: mem_compressed = spbas_memory_requirement(*B);
353: *mem_reduction = 100.0 * (PetscReal)(mem_orig-mem_compressed)/ (PetscReal) mem_orig;
354: return(0);
355: }
357: /*
358: spbas_incomplete_cholesky
359: Incomplete Cholesky decomposition
360: */
361: #include <../src/mat/impls/aij/seq/bas/spbas_cholesky.h>
363: /*
364: spbas_delete : de-allocate the arrays owned by this matrix
365: */
366: PetscErrorCode spbas_delete(spbas_matrix matrix)
367: {
368: PetscInt i;
372: if (matrix.block_data) {
373: ierr=PetscFree(matrix.alloc_icol);
374: if (matrix.values) {ierr=PetscFree(matrix.alloc_val);}
375: } else {
376: for (i=0; i<matrix.nrows; i++) { ierr=PetscFree(matrix.icols[i]);}
377: PetscFree(matrix.icols);
378: if (matrix.values) {
379: for (i=0; i<matrix.nrows; i++) { ierr=PetscFree(matrix.values[i]);}
380: }
381: }
383: ierr=PetscFree(matrix.row_nnz);
384: ierr=PetscFree(matrix.icols);
385: if (matrix.col_idx_type == SPBAS_OFFSET_ARRAY) {ierr=PetscFree(matrix.icol0);}
386: ierr=PetscFree(matrix.values);
387: return(0);
388: }
390: /*
391: spbas_matrix_to_crs:
392: Convert an spbas_matrix to compessed row storage
393: */
394: PetscErrorCode spbas_matrix_to_crs(spbas_matrix matrix_A,MatScalar **val_out, PetscInt **irow_out, PetscInt **icol_out)
395: {
396: PetscInt nrows = matrix_A.nrows;
397: PetscInt nnz = matrix_A.nnz;
398: PetscInt i,j,r_nnz,i0;
399: PetscInt *irow;
400: PetscInt *icol;
401: PetscInt *icol_A;
402: MatScalar *val;
403: PetscScalar *val_A;
404: PetscInt col_idx_type = matrix_A.col_idx_type;
405: PetscBool do_values = matrix_A.values ? PETSC_TRUE : PETSC_FALSE;
409: PetscMalloc1(nrows+1, &irow);
410: PetscMalloc1(nnz, &icol);
411: *icol_out = icol;
412: *irow_out = irow;
413: if (do_values) {
414: PetscMalloc1(nnz, &val);
415: *val_out = val; *icol_out = icol; *irow_out=irow;
416: }
418: irow[0]=0;
419: for (i=0; i<nrows; i++) {
420: r_nnz = matrix_A.row_nnz[i];
421: i0 = irow[i];
422: irow[i+1] = i0 + r_nnz;
423: icol_A = matrix_A.icols[i];
425: if (do_values) {
426: val_A = matrix_A.values[i];
427: for (j=0; j<r_nnz; j++) {
428: icol[i0+j] = icol_A[j];
429: val[i0+j] = val_A[j];
430: }
431: } else {
432: for (j=0; j<r_nnz; j++) icol[i0+j] = icol_A[j];
433: }
435: if (col_idx_type == SPBAS_DIAGONAL_OFFSETS) {
436: for (j=0; j<r_nnz; j++) icol[i0+j] += i;
437: } else if (col_idx_type == SPBAS_OFFSET_ARRAY) {
438: i0 = matrix_A.icol0[i];
439: for (j=0; j<r_nnz; j++) icol[i0+j] += i0;
440: }
441: }
442: return(0);
443: }
446: /*
447: spbas_transpose
448: return the transpose of a matrix
449: */
450: PetscErrorCode spbas_transpose(spbas_matrix in_matrix, spbas_matrix * result)
451: {
452: PetscInt col_idx_type = in_matrix.col_idx_type;
453: PetscInt nnz = in_matrix.nnz;
454: PetscInt ncols = in_matrix.nrows;
455: PetscInt nrows = in_matrix.ncols;
456: PetscInt i,j,k;
457: PetscInt r_nnz;
458: PetscInt *irow;
459: PetscInt icol0 = 0;
460: PetscScalar * val;
464: /* Copy input values */
465: result->nrows = nrows;
466: result->ncols = ncols;
467: result->nnz = nnz;
468: result->col_idx_type = SPBAS_COLUMN_NUMBERS;
469: result->block_data = PETSC_TRUE;
471: /* Allocate sparseness pattern */
472: spbas_allocate_pattern(result, in_matrix.values ? PETSC_TRUE : PETSC_FALSE);
474: /* Count the number of nonzeros in each row */
475: for (i = 0; i<nrows; i++) result->row_nnz[i] = 0;
477: for (i=0; i<ncols; i++) {
478: r_nnz = in_matrix.row_nnz[i];
479: irow = in_matrix.icols[i];
480: if (col_idx_type == SPBAS_COLUMN_NUMBERS) {
481: for (j=0; j<r_nnz; j++) result->row_nnz[irow[j]]++;
482: } else if (col_idx_type == SPBAS_DIAGONAL_OFFSETS) {
483: for (j=0; j<r_nnz; j++) result->row_nnz[i+irow[j]]++;
484: } else if (col_idx_type == SPBAS_OFFSET_ARRAY) {
485: icol0=in_matrix.icol0[i];
486: for (j=0; j<r_nnz; j++) result->row_nnz[icol0+irow[j]]++;
487: }
488: }
490: /* Set the pointers to the data */
491: spbas_allocate_data(result);
493: /* Reset the number of nonzeros in each row */
494: for (i = 0; i<nrows; i++) result->row_nnz[i] = 0;
496: /* Fill the data arrays */
497: if (in_matrix.values) {
498: for (i=0; i<ncols; i++) {
499: r_nnz = in_matrix.row_nnz[i];
500: irow = in_matrix.icols[i];
501: val = in_matrix.values[i];
503: if (col_idx_type == SPBAS_COLUMN_NUMBERS) icol0 = 0;
504: else if (col_idx_type == SPBAS_DIAGONAL_OFFSETS) icol0 = i;
505: else if (col_idx_type == SPBAS_OFFSET_ARRAY) icol0 = in_matrix.icol0[i];
506: for (j=0; j<r_nnz; j++) {
507: k = icol0 + irow[j];
508: result->icols[k][result->row_nnz[k]] = i;
509: result->values[k][result->row_nnz[k]] = val[j];
510: result->row_nnz[k]++;
511: }
512: }
513: } else {
514: for (i=0; i<ncols; i++) {
515: r_nnz = in_matrix.row_nnz[i];
516: irow = in_matrix.icols[i];
518: if (col_idx_type == SPBAS_COLUMN_NUMBERS) icol0=0;
519: else if (col_idx_type == SPBAS_DIAGONAL_OFFSETS) icol0=i;
520: else if (col_idx_type == SPBAS_OFFSET_ARRAY) icol0=in_matrix.icol0[i];
522: for (j=0; j<r_nnz; j++) {
523: k = icol0 + irow[j];
524: result->icols[k][result->row_nnz[k]] = i;
525: result->row_nnz[k]++;
526: }
527: }
528: }
529: return(0);
530: }
532: /*
533: spbas_mergesort
535: mergesort for an array of integers and an array of associated
536: reals
538: on output, icol[0..nnz-1] is increasing;
539: val[0..nnz-1] has undergone the same permutation as icol
541: NB: val may be NULL: in that case, only the integers are sorted
543: */
544: PetscErrorCode spbas_mergesort(PetscInt nnz, PetscInt *icol, PetscScalar *val)
545: {
546: PetscInt istep; /* Chunk-sizes of already sorted parts of arrays */
547: PetscInt i, i1, i2; /* Loop counters for (partly) sorted arrays */
548: PetscInt istart, i1end, i2end; /* start of newly sorted array part, end of both parts */
549: PetscInt *ialloc; /* Allocated arrays */
550: PetscScalar *valloc=NULL;
551: PetscInt *iswap; /* auxiliary pointers for swapping */
552: PetscScalar *vswap;
553: PetscInt *ihlp1; /* Pointers to new version of arrays, */
554: PetscScalar *vhlp1=NULL; /* (arrays under construction) */
555: PetscInt *ihlp2; /* Pointers to previous version of arrays, */
556: PetscScalar *vhlp2=NULL;
559: PetscMalloc1(nnz,&ialloc);
560: ihlp1 = ialloc;
561: ihlp2 = icol;
563: if (val) {
564: PetscMalloc1(nnz,&valloc);
565: vhlp1 = valloc;
566: vhlp2 = val;
567: }
570: /* Sorted array chunks are first 1 long, and increase until they are the complete array */
571: for (istep=1; istep<nnz; istep*=2) {
572: /*
573: Combine sorted parts
574: istart:istart+istep-1 and istart+istep-1:istart+2*istep-1
575: of ihlp2 and vhlp2
577: into one sorted part
578: istart:istart+2*istep-1
579: of ihlp1 and vhlp1
580: */
581: for (istart=0; istart<nnz; istart+=2*istep) {
582: /* Set counters and bound array part endings */
583: i1=istart; i1end = i1+istep; if (i1end>nnz) i1end=nnz;
584: i2=istart+istep; i2end = i2+istep; if (i2end>nnz) i2end=nnz;
586: /* Merge the two array parts */
587: if (val) {
588: for (i=istart; i<i2end; i++) {
589: if (i1<i1end && i2<i2end && ihlp2[i1] < ihlp2[i2]) {
590: ihlp1[i] = ihlp2[i1];
591: vhlp1[i] = vhlp2[i1];
592: i1++;
593: } else if (i2<i2end) {
594: ihlp1[i] = ihlp2[i2];
595: vhlp1[i] = vhlp2[i2];
596: i2++;
597: } else {
598: ihlp1[i] = ihlp2[i1];
599: vhlp1[i] = vhlp2[i1];
600: i1++;
601: }
602: }
603: } else {
604: for (i=istart; i<i2end; i++) {
605: if (i1<i1end && i2<i2end && ihlp2[i1] < ihlp2[i2]) {
606: ihlp1[i] = ihlp2[i1];
607: i1++;
608: } else if (i2<i2end) {
609: ihlp1[i] = ihlp2[i2];
610: i2++;
611: } else {
612: ihlp1[i] = ihlp2[i1];
613: i1++;
614: }
615: }
616: }
617: }
619: /* Swap the two array sets */
620: iswap = ihlp2; ihlp2 = ihlp1; ihlp1 = iswap;
621: vswap = vhlp2; vhlp2 = vhlp1; vhlp1 = vswap;
622: }
624: /* Copy one more time in case the sorted arrays are the temporary ones */
625: if (ihlp2 != icol) {
626: for (i=0; i<nnz; i++) icol[i] = ihlp2[i];
627: if (val) {
628: for (i=0; i<nnz; i++) val[i] = vhlp2[i];
629: }
630: }
632: PetscFree(ialloc);
633: if (val) {PetscFree(valloc);}
634: return(0);
635: }
637: /*
638: spbas_apply_reordering_rows:
639: apply the given reordering to the rows: matrix_A = matrix_A(perm,:);
640: */
641: PetscErrorCode spbas_apply_reordering_rows(spbas_matrix *matrix_A, const PetscInt *permutation)
642: {
643: PetscInt i,j,ip;
644: PetscInt nrows=matrix_A->nrows;
645: PetscInt * row_nnz;
646: PetscInt **icols;
647: PetscBool do_values = matrix_A->values ? PETSC_TRUE : PETSC_FALSE;
648: PetscScalar **vals = NULL;
652: if (matrix_A->col_idx_type != SPBAS_DIAGONAL_OFFSETS) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP_SYS,"must have diagonal offsets in pattern\n");
654: if (do_values) {
655: PetscMalloc1(nrows, &vals);
656: }
657: PetscMalloc1(nrows, &row_nnz);
658: PetscMalloc1(nrows, &icols);
660: for (i=0; i<nrows; i++) {
661: ip = permutation[i];
662: if (do_values) vals[i] = matrix_A->values[ip];
663: icols[i] = matrix_A->icols[ip];
664: row_nnz[i] = matrix_A->row_nnz[ip];
665: for (j=0; j<row_nnz[i]; j++) icols[i][j] += ip-i;
666: }
668: if (do_values) { PetscFree(matrix_A->values);}
669: PetscFree(matrix_A->icols);
670: PetscFree(matrix_A->row_nnz);
672: if (do_values) matrix_A->values = vals;
673: matrix_A->icols = icols;
674: matrix_A->row_nnz = row_nnz;
675: return(0);
676: }
679: /*
680: spbas_apply_reordering_cols:
681: apply the given reordering to the columns: matrix_A(:,perm) = matrix_A;
682: */
683: PetscErrorCode spbas_apply_reordering_cols(spbas_matrix *matrix_A,const PetscInt *permutation)
684: {
685: PetscInt i,j;
686: PetscInt nrows=matrix_A->nrows;
687: PetscInt row_nnz;
688: PetscInt *icols;
689: PetscBool do_values = matrix_A->values ? PETSC_TRUE : PETSC_FALSE;
690: PetscScalar *vals = NULL;
694: if (matrix_A->col_idx_type != SPBAS_DIAGONAL_OFFSETS) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP_SYS, "must have diagonal offsets in pattern\n");
696: for (i=0; i<nrows; i++) {
697: icols = matrix_A->icols[i];
698: row_nnz = matrix_A->row_nnz[i];
699: if (do_values) vals = matrix_A->values[i];
701: for (j=0; j<row_nnz; j++) {
702: icols[j] = permutation[i+icols[j]]-i;
703: }
704: spbas_mergesort(row_nnz, icols, vals);
705: }
706: return(0);
707: }
709: /*
710: spbas_apply_reordering:
711: apply the given reordering: matrix_A(perm,perm) = matrix_A;
712: */
713: PetscErrorCode spbas_apply_reordering(spbas_matrix *matrix_A, const PetscInt *permutation, const PetscInt * inv_perm)
714: {
718: spbas_apply_reordering_rows(matrix_A, inv_perm);
719: spbas_apply_reordering_cols(matrix_A, permutation);
720: return(0);
721: }
723: PetscErrorCode spbas_pattern_only(PetscInt nrows, PetscInt ncols, PetscInt *ai, PetscInt *aj, spbas_matrix * result)
724: {
725: spbas_matrix retval;
726: PetscInt i, j, i0, r_nnz;
730: /* Copy input values */
731: retval.nrows = nrows;
732: retval.ncols = ncols;
733: retval.nnz = ai[nrows];
735: retval.block_data = PETSC_TRUE;
736: retval.col_idx_type = SPBAS_DIAGONAL_OFFSETS;
738: /* Allocate output matrix */
739: spbas_allocate_pattern(&retval, PETSC_FALSE);
740: for (i=0; i<nrows; i++) retval.row_nnz[i] = ai[i+1]-ai[i];
741: spbas_allocate_data(&retval);
742: /* Copy the structure */
743: for (i = 0; i<retval.nrows; i++) {
744: i0 = ai[i];
745: r_nnz = ai[i+1]-i0;
747: for (j=0; j<r_nnz; j++) {
748: retval.icols[i][j] = aj[i0+j]-i;
749: }
750: }
751: *result = retval;
752: return(0);
753: }
756: /*
757: spbas_mark_row_power:
758: Mark the columns in row 'row' which are nonzero in
759: matrix^2log(marker).
760: */
761: PetscErrorCode spbas_mark_row_power(PetscInt *iwork, /* marker-vector */
762: PetscInt row, /* row for which the columns are marked */
763: spbas_matrix * in_matrix, /* matrix for which the power is being calculated */
764: PetscInt marker, /* marker-value: 2^power */
765: PetscInt minmrk, /* lower bound for marked points */
766: PetscInt maxmrk) /* upper bound for marked points */
767: {
769: PetscInt i,j, nnz;
772: nnz = in_matrix->row_nnz[row];
774: /* For higher powers, call this function recursively */
775: if (marker>1) {
776: for (i=0; i<nnz; i++) {
777: j = row + in_matrix->icols[row][i];
778: if (minmrk<=j && j<maxmrk && iwork[j] < marker) {
779: spbas_mark_row_power(iwork, row + in_matrix->icols[row][i],in_matrix, marker/2,minmrk,maxmrk);
780: iwork[j] |= marker;
781: }
782: }
783: } else {
784: /* Mark the columns reached */
785: for (i=0; i<nnz; i++) {
786: j = row + in_matrix->icols[row][i];
787: if (minmrk<=j && j<maxmrk) iwork[j] |= 1;
788: }
789: }
790: return(0);
791: }
794: /*
795: spbas_power
796: Calculate sparseness patterns for incomplete Cholesky decompositions
797: of a given order: (almost) all nonzeros of the matrix^(order+1) which
798: are inside the band width are found and stored in the output sparseness
799: pattern.
800: */
801: PetscErrorCode spbas_power(spbas_matrix in_matrix,PetscInt power, spbas_matrix * result)
802: {
803: spbas_matrix retval;
804: PetscInt nrows = in_matrix.nrows;
805: PetscInt ncols = in_matrix.ncols;
806: PetscInt i, j, kend;
807: PetscInt nnz, inz;
808: PetscInt *iwork;
809: PetscInt marker;
810: PetscInt maxmrk=0;
814: if (in_matrix.col_idx_type != SPBAS_DIAGONAL_OFFSETS) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP_SYS,"must have diagonal offsets in pattern\n");
815: if (ncols != nrows) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "Dimension error\n");
816: if (in_matrix.values) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "Input array must be sparseness pattern (no values)");
817: if (power<=0) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP_SYS, "Power must be 1 or up");
819: /* Copy input values*/
820: retval.nrows = ncols;
821: retval.ncols = nrows;
822: retval.nnz = 0;
823: retval.col_idx_type = SPBAS_DIAGONAL_OFFSETS;
824: retval.block_data = PETSC_FALSE;
826: /* Allocate sparseness pattern */
827: spbas_allocate_pattern(&retval, in_matrix.values ? PETSC_TRUE : PETSC_FALSE);
829: /* Allocate marker array */
830: PetscMalloc1(nrows, &iwork);
832: /* Erase the pattern for this row */
833: PetscMemzero((void*) iwork, retval.nrows*sizeof(PetscInt));
835: /* Calculate marker values */
836: marker = 1; for (i=1; i<power; i++) marker*=2;
838: for (i=0; i<nrows; i++) {
839: /* Calculate the pattern for each row */
841: nnz = in_matrix.row_nnz[i];
842: kend = i+in_matrix.icols[i][nnz-1];
843: if (maxmrk<=kend) maxmrk=kend+1;
844: spbas_mark_row_power(iwork, i, &in_matrix, marker, i, maxmrk);
846: /* Count the columns*/
847: nnz = 0;
848: for (j=i; j<maxmrk; j++) nnz+= (iwork[j]!=0);
850: /* Allocate the column indices */
851: retval.row_nnz[i] = nnz;
852: PetscMalloc1(nnz,&retval.icols[i]);
854: /* Administrate the column indices */
855: inz = 0;
856: for (j=i; j<maxmrk; j++) {
857: if (iwork[j]) {
858: retval.icols[i][inz] = j-i;
859: inz++;
860: iwork[j]=0;
861: }
862: }
863: retval.nnz += nnz;
864: };
865: PetscFree(iwork);
866: *result = retval;
867: return(0);
868: }
872: /*
873: spbas_keep_upper:
874: remove the lower part of the matrix: keep the upper part
875: */
876: PetscErrorCode spbas_keep_upper(spbas_matrix * inout_matrix)
877: {
878: PetscInt i, j;
879: PetscInt jstart;
882: if (inout_matrix->block_data) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP_SYS, "Not yet for block data matrices\n");
883: for (i=0; i<inout_matrix->nrows; i++) {
884: for (jstart=0; (jstart<inout_matrix->row_nnz[i]) && (inout_matrix->icols[i][jstart]<0); jstart++) {}
885: if (jstart>0) {
886: for (j=0; j<inout_matrix->row_nnz[i]-jstart; j++) {
887: inout_matrix->icols[i][j] = inout_matrix->icols[i][j+jstart];
888: }
890: if (inout_matrix->values) {
891: for (j=0; j<inout_matrix->row_nnz[i]-jstart; j++) {
892: inout_matrix->values[i][j] = inout_matrix->values[i][j+jstart];
893: }
894: }
896: inout_matrix->row_nnz[i] -= jstart;
898: inout_matrix->icols[i] = (PetscInt*) realloc((void*) inout_matrix->icols[i], inout_matrix->row_nnz[i]*sizeof(PetscInt));
900: if (inout_matrix->values) {
901: inout_matrix->values[i] = (PetscScalar*) realloc((void*) inout_matrix->values[i], inout_matrix->row_nnz[i]*sizeof(PetscScalar));
902: }
903: inout_matrix->nnz -= jstart;
904: }
905: }
906: return(0);
907: }