Actual source code: spbas.c
petsc-3.3-p7 2013-05-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>
11: - -pc_factor_drop_tolerance
13: Level: beginner
15: Contributed by: Bas van 't Hof
17: .seealso: PCFactorSetMatSolverPackage(), MatSolverPackage, PCFactorSetLevels(), PCFactorSetDropTolerance()
19: M*/
21: /*
22: spbas_memory_requirement:
23: Calculate the number of bytes needed to store tha matrix
24: */
27: long int spbas_memory_requirement( spbas_matrix matrix)
28: {
29: long int memreq = 6 * sizeof(PetscInt) + /* nrows, ncols, nnz, n_alloc_icol,
30: n_alloc_val, col_idx_type */
31: sizeof(PetscBool) + /* block_data */
32: sizeof(PetscScalar**) + /* values */
33: sizeof(PetscScalar*) + /* alloc_val */
34: 2 * sizeof(int**) + /* icols, icols0 */
35: 2 * sizeof(PetscInt*) + /* row_nnz, alloc_icol */
36: matrix.nrows * sizeof(PetscInt) + /* row_nnz[*] */
37: matrix.nrows * sizeof(PetscInt*); /* icols[*] */
39: /* icol0[*] */
40: if (matrix.col_idx_type == SPBAS_OFFSET_ARRAY) { memreq += matrix.nrows * sizeof(PetscInt); }
42: /* icols[*][*] */
43: if (matrix.block_data) { memreq += matrix.n_alloc_icol * sizeof(PetscInt); }
44: else { memreq += matrix.nnz * sizeof(PetscInt); }
46: if (matrix.values) {
47: memreq += matrix.nrows * sizeof(PetscScalar*); /* values[*] */
48: /* values[*][*] */
49: if (matrix.block_data) { memreq += matrix.n_alloc_val * sizeof(PetscScalar); }
50: else { memreq += matrix.nnz * sizeof(PetscScalar); }
51: }
52: return memreq;
53: }
55: /*
56: spbas_allocate_pattern:
57: allocate the pattern arrays row_nnz, icols and optionally values
58: */
61: PetscErrorCode spbas_allocate_pattern( spbas_matrix * result, PetscBool do_values)
62: {
64: PetscInt nrows = result->nrows;
65: PetscInt col_idx_type = result->col_idx_type;
68: /* Allocate sparseness pattern */
69: PetscMalloc(nrows*sizeof(PetscInt),&result->row_nnz);
70: PetscMalloc(nrows*sizeof(PetscInt*),&result->icols);
72: /* If offsets are given wrt an array, create array */
73: if (col_idx_type == SPBAS_OFFSET_ARRAY) {
74: PetscMalloc(nrows*sizeof(PetscInt),&result->icol0);
75: } else {
76: result->icol0 = PETSC_NULL;
77: }
78:
79: /* If values are given, allocate values array */
80: if (do_values) {
81: PetscMalloc(nrows*sizeof(PetscScalar*),&result->values);
82: } else {
83: result->values = PETSC_NULL;
84: }
85: return(0);
86: }
88: /*
89: spbas_allocate_data:
90: in case of block_data:
91: Allocate the data arrays alloc_icol and optionally alloc_val,
92: set appropriate pointers from icols and values;
93: in case of !block_data:
94: Allocate the arrays icols[i] and optionally values[i]
95: */
98: PetscErrorCode spbas_allocate_data( spbas_matrix * result)
99: {
100: PetscInt i;
101: PetscInt nnz = result->nnz;
102: PetscInt nrows = result->nrows;
103: PetscInt r_nnz;
105: PetscBool do_values = (result->values != PETSC_NULL) ? PETSC_TRUE : PETSC_FALSE;
106: PetscBool block_data = result->block_data;
109: if (block_data) {
110: /* Allocate the column number array and point to it */
111: result->n_alloc_icol = nnz;
112: PetscMalloc(nnz*sizeof(PetscInt), &result->alloc_icol);
113: result->icols[0]= result->alloc_icol;
114: for (i=1; i<nrows; i++) {
115: result->icols[i] = result->icols[i-1] + result->row_nnz[i-1];
116: }
118: /* Allocate the value array and point to it */
119: if (do_values) {
120: result->n_alloc_val= nnz;
121: PetscMalloc(nnz*sizeof(PetscScalar), &result->alloc_val);
122: result->values[0]= result->alloc_val;
123: for (i=1; i<nrows; i++) {
124: result->values[i] = result->values[i-1] + result->row_nnz[i-1];
125: }
126: }
127: } else {
128: for (i=0; i<nrows; i++) {
129: r_nnz = result->row_nnz[i];
130: PetscMalloc(r_nnz*sizeof(PetscInt), &result->icols[i]);
131: }
132: if (do_values) {
133: for (i=0; i<nrows; i++) {
134: r_nnz = result->row_nnz[i];
135: PetscMalloc(r_nnz*sizeof(PetscScalar), &result->values[i]);
136: }
137: }
138: }
139: return(0);
140: }
142: /*
143: spbas_row_order_icol
144: determine if row i1 should come
145: + before row i2 in the sorted rows (return -1),
146: + after (return 1)
147: + is identical (return 0).
148: */
151: int spbas_row_order_icol(PetscInt i1, PetscInt i2, PetscInt *irow_in, PetscInt *icol_in,PetscInt col_idx_type)
152: {
153: PetscInt j;
154: PetscInt nnz1 = irow_in[i1+1] - irow_in[i1];
155: PetscInt nnz2 = irow_in[i2+1] - irow_in[i2];
156: PetscInt * icol1 = &icol_in[irow_in[i1]];
157: PetscInt * icol2 = &icol_in[irow_in[i2]];
159: if (nnz1<nnz2) {return -1;}
160: if (nnz1>nnz2) {return 1;}
162: if (col_idx_type == SPBAS_COLUMN_NUMBERS) {
163: for (j=0; j<nnz1; j++) {
164: if (icol1[j]< icol2[j]) {return -1;}
165: if (icol1[j]> icol2[j]) {return 1;}
166: }
167: } else if (col_idx_type == SPBAS_DIAGONAL_OFFSETS) {
168: for (j=0; j<nnz1; j++) {
169: if (icol1[j]-i1< icol2[j]-i2) {return -1;}
170: if (icol1[j]-i1> icol2[j]-i2) {return 1;}
171: }
172: } else if (col_idx_type == SPBAS_OFFSET_ARRAY) {
173: for (j=1; j<nnz1; j++) {
174: if (icol1[j]-icol1[0] < icol2[j]-icol2[0]) {return -1;}
175: if (icol1[j]-icol1[0] > icol2[j]-icol2[0]) {return 1;}
176: }
177: }
178: return 0;
179: }
181: /*
182: spbas_mergesort_icols:
183: return a sorting of the rows in which identical sparseness patterns are
184: next to each other
185: */
188: PetscErrorCode spbas_mergesort_icols( PetscInt nrows, PetscInt * irow_in, PetscInt * icol_in,PetscInt col_idx_type, PetscInt *isort)
189: {
191: PetscInt istep; /* Chunk-sizes of already sorted parts of arrays */
192: PetscInt i, i1, i2; /* Loop counters for (partly) sorted arrays */
193: PetscInt istart, i1end, i2end; /* start of newly sorted array part, end of both parts */
194: PetscInt *ialloc; /* Allocated arrays */
195: PetscInt *iswap; /* auxiliary pointers for swapping */
196: PetscInt *ihlp1; /* Pointers to new version of arrays, */
197: PetscInt *ihlp2; /* Pointers to previous version of arrays, */
201: PetscMalloc(nrows*sizeof(PetscInt),&ialloc);
203: ihlp1 = ialloc;
204: ihlp2 = isort;
206: /* Sorted array chunks are first 1 long, and increase until they are the complete array */
207: for (istep=1; istep<nrows; istep*=2) {
208: /*
209: Combine sorted parts
210: istart:istart+istep-1 and istart+istep-1:istart+2*istep-1
211: of ihlp2 and vhlp2
212:
213: into one sorted part
214: istart:istart+2*istep-1
215: of ihlp1 and vhlp1
216: */
217: for (istart=0; istart<nrows; istart+=2*istep) {
219: /* Set counters and bound array part endings */
220: i1=istart; i1end = i1+istep; if (i1end>nrows) {i1end=nrows;}
221: i2=istart+istep; i2end = i2+istep; if (i2end>nrows) {i2end=nrows;}
223: /* Merge the two array parts */
224: for (i=istart; i<i2end; i++) {
225: if (i1<i1end && i2<i2end && spbas_row_order_icol( ihlp2[i1], ihlp2[i2], irow_in, icol_in, col_idx_type) < 0) {
226: ihlp1[i] = ihlp2[i1];
227: i1++;
228: } else if (i2<i2end ) {
229: ihlp1[i] = ihlp2[i2];
230: i2++;
231: } else {
232: ihlp1[i] = ihlp2[i1];
233: i1++;
234: }
235: }
236: }
238: /* Swap the two array sets */
239: iswap = ihlp2; ihlp2 = ihlp1; ihlp1 = iswap;
240: }
242: /* Copy one more time in case the sorted arrays are the temporary ones */
243: if (ihlp2 != isort) {
244: for (i=0; i<nrows; i++) { isort[i] = ihlp2[i]; }
245: }
246: PetscFree(ialloc);
247: return(0);
248: }
252: /*
253: spbas_compress_pattern:
254: calculate a compressed sparseness pattern for a sparseness pattern
255: given in compressed row storage. The compressed sparseness pattern may
256: require (much) less memory.
257: */
260: PetscErrorCode spbas_compress_pattern(PetscInt *irow_in, PetscInt *icol_in, PetscInt nrows, PetscInt ncols, PetscInt col_idx_type, spbas_matrix *B,PetscReal *mem_reduction)
261: {
262: PetscInt nnz = irow_in[nrows];
263: long int mem_orig = (nrows + nnz) * sizeof(PetscInt);
264: long int mem_compressed;
265: PetscErrorCode ierr;
266: PetscInt *isort;
267: PetscInt *icols;
268: PetscInt row_nnz;
269: PetscInt *ipoint;
270: PetscBool *used;
271: PetscInt ptr;
272: PetscInt i,j;
273: const PetscBool no_values = PETSC_FALSE;
276: /* Allocate the structure of the new matrix */
277: B->nrows = nrows;
278: B->ncols = ncols;
279: B->nnz = nnz;
280: B->col_idx_type= col_idx_type;
281: B->block_data = PETSC_TRUE;
282: spbas_allocate_pattern( B, no_values);
284: /* When using an offset array, set it */
285: if (col_idx_type==SPBAS_OFFSET_ARRAY) {
286: for (i=0; i<nrows; i++) {B->icol0[i] = icol_in[irow_in[i]];}
287: }
289: /* Allocate the ordering for the rows */
290: PetscMalloc(nrows*sizeof(PetscInt),&isort);
291: PetscMalloc(nrows*sizeof(PetscInt),&ipoint);
292: PetscMalloc(nrows*sizeof(PetscBool),&used);
294: /* Initialize the sorting */
295: PetscMemzero((void*) used, nrows*sizeof(PetscBool));
296: for (i = 0; i<nrows; i++) {
297: B->row_nnz[i] = irow_in[i+1]-irow_in[i];
298: isort[i] = i;
299: ipoint[i]= i;
300: }
302: /* Sort the rows so that identical columns will be next to each other */
303: spbas_mergesort_icols( nrows, irow_in, icol_in, col_idx_type, isort);
304: PetscInfo(PETSC_NULL,"Rows have been sorted for patterns\n");
306: /* Replace identical rows with the first one in the list */
307: for (i=1; i<nrows; i++) {
308: if (spbas_row_order_icol(isort[i-1], isort[i], irow_in, icol_in, col_idx_type) == 0) {
309: ipoint[isort[i]] = ipoint[isort[i-1]];
310: }
311: }
313: /* Collect the rows which are used*/
314: for (i=0; i<nrows; i++) {used[ipoint[i]] = PETSC_TRUE;}
316: /* Calculate needed memory */
317: B->n_alloc_icol = 0;
318: for (i=0; i<nrows; i++) {
319: if (used[i]) {B->n_alloc_icol += B->row_nnz[i];}
320: }
321: PetscMalloc(B->n_alloc_icol*sizeof(PetscInt),&B->alloc_icol);
323: /* Fill in the diagonal offsets for the rows which store their own data */
324: ptr = 0;
325: for (i=0; i<B->nrows; i++) {
326: if (used[i]) {
327: B->icols[i] = &B->alloc_icol[ptr];
328: icols = &icol_in[irow_in[i]];
329: row_nnz = B->row_nnz[i];
330: if (col_idx_type == SPBAS_COLUMN_NUMBERS) {
331: for (j=0; j<row_nnz; j++) {
332: B->icols[i][j] = icols[j];
333: }
334: } else if (col_idx_type == SPBAS_DIAGONAL_OFFSETS) {
335: for (j=0; j<row_nnz; j++) {
336: B->icols[i][j] = icols[j]-i;
337: }
338: } else if (col_idx_type == SPBAS_OFFSET_ARRAY) {
339: for (j=0; j<row_nnz; j++) {
340: B->icols[i][j] = icols[j]-icols[0];
341: }
342: }
343: ptr += B->row_nnz[i];
344: }
345: }
347: /* Point to the right places for all data */
348: for (i=0; i<nrows; i++) {
349: B->icols[i] = B->icols[ipoint[i]];
350: }
351: PetscInfo(PETSC_NULL,"Row patterns have been compressed\n");
352: PetscInfo1(PETSC_NULL," (%G nonzeros per row)\n", (PetscReal) nnz / (PetscReal) nrows);
353:
354: ierr=PetscFree(isort);
355: ierr=PetscFree(used);
356: ierr=PetscFree(ipoint);
358: mem_compressed = spbas_memory_requirement( *B );
359: *mem_reduction = 100.0 * (PetscReal)(mem_orig-mem_compressed)/ (PetscReal) mem_orig;
360: return(0);
361: }
363: /*
364: spbas_incomplete_cholesky
365: Incomplete Cholesky decomposition
366: */
367: #include <../src/mat/impls/aij/seq/bas/spbas_cholesky.h>
369: /*
370: spbas_delete : de-allocate the arrays owned by this matrix
371: */
374: PetscErrorCode spbas_delete(spbas_matrix matrix)
375: {
376: PetscInt i;
380: if (matrix.block_data) {
381: ierr=PetscFree(matrix.alloc_icol);
382: if (matrix.values){ierr=PetscFree(matrix.alloc_val);}
383: } else {
384: for (i=0; i<matrix.nrows; i++) { ierr=PetscFree(matrix.icols[i]);}
385: PetscFree(matrix.icols);
386: if (matrix.values) {
387: for (i=0; i<matrix.nrows; i++) { ierr=PetscFree(matrix.values[i]);}
388: }
389: }
391: ierr=PetscFree(matrix.row_nnz);
392: ierr=PetscFree(matrix.icols);
393: if (matrix.col_idx_type == SPBAS_OFFSET_ARRAY) {ierr=PetscFree(matrix.icol0);}
394: ierr=PetscFree(matrix.values);
395: return(0);
396: }
398: /*
399: spbas_matrix_to_crs:
400: Convert an spbas_matrix to compessed row storage
401: */
404: PetscErrorCode spbas_matrix_to_crs(spbas_matrix matrix_A,MatScalar **val_out, PetscInt **irow_out, PetscInt **icol_out)
405: {
406: PetscInt nrows = matrix_A.nrows;
407: PetscInt nnz = matrix_A.nnz;
408: PetscInt i,j,r_nnz,i0;
409: PetscInt *irow;
410: PetscInt *icol;
411: PetscInt *icol_A;
412: MatScalar *val;
413: PetscScalar *val_A;
414: PetscInt col_idx_type = matrix_A.col_idx_type;
415: PetscBool do_values = matrix_A.values ? PETSC_TRUE : PETSC_FALSE;
419: PetscMalloc( sizeof(PetscInt) * (nrows+1), &irow);
420: PetscMalloc( sizeof(PetscInt) * nnz, &icol);
421: *icol_out = icol;
422: *irow_out=irow;
423: if (do_values) {
424: PetscMalloc( sizeof(MatScalar) * nnz, &val);
425: *val_out = val; *icol_out = icol; *irow_out=irow;
426: }
428: irow[0]=0;
429: for (i=0; i<nrows; i++) {
430: r_nnz = matrix_A.row_nnz[i];
431: i0 = irow[i];
432: irow[i+1] = i0 + r_nnz;
433: icol_A = matrix_A.icols[i];
435: if (do_values) {
436: val_A = matrix_A.values[i];
437: for (j=0; j<r_nnz; j++) {
438: icol[i0+j] = icol_A[j];
439: val[i0+j] = val_A[j];
440: }
441: } else {
442: for (j=0; j<r_nnz; j++) { icol[i0+j] = icol_A[j]; }
443: }
445: if (col_idx_type == SPBAS_DIAGONAL_OFFSETS) {
446: for (j=0; j<r_nnz; j++) { icol[i0+j] += i; }
447: }
448: else if (col_idx_type == SPBAS_OFFSET_ARRAY) {
449: i0 = matrix_A.icol0[i];
450: for (j=0; j<r_nnz; j++) { icol[i0+j] += i0; }
451: }
452: }
453: return(0);
454: }
457: /*
458: spbas_transpose
459: return the transpose of a matrix
460: */
463: PetscErrorCode spbas_transpose( spbas_matrix in_matrix, spbas_matrix * result)
464: {
465: PetscInt col_idx_type = in_matrix.col_idx_type;
466: PetscInt nnz = in_matrix.nnz;
467: PetscInt ncols = in_matrix.nrows;
468: PetscInt nrows = in_matrix.ncols;
469: PetscInt i,j,k;
470: PetscInt r_nnz;
471: PetscInt *irow;
472: PetscInt icol0 = 0;
473: PetscScalar * val;
477: /* Copy input values */
478: result->nrows = nrows;
479: result->ncols = ncols;
480: result->nnz = nnz;
481: result->col_idx_type = SPBAS_COLUMN_NUMBERS;
482: result->block_data = PETSC_TRUE;
484: /* Allocate sparseness pattern */
485: spbas_allocate_pattern(result, in_matrix.values ? PETSC_TRUE : PETSC_FALSE);
487: /* Count the number of nonzeros in each row */
488: for (i = 0; i<nrows; i++) { result->row_nnz[i] = 0; }
490: for (i=0; i<ncols; i++) {
491: r_nnz = in_matrix.row_nnz[i];
492: irow = in_matrix.icols[i];
493: if (col_idx_type == SPBAS_COLUMN_NUMBERS) {
494: for (j=0; j<r_nnz; j++) { result->row_nnz[irow[j]]++; }
495: } else if (col_idx_type == SPBAS_DIAGONAL_OFFSETS) {
496: for (j=0; j<r_nnz; j++) { result->row_nnz[i+irow[j]]++; }
497: } else if (col_idx_type == SPBAS_OFFSET_ARRAY) {
498: icol0=in_matrix.icol0[i];
499: for (j=0; j<r_nnz; j++) { result->row_nnz[icol0+irow[j]]++; }
500: }
501: }
503: /* Set the pointers to the data */
504: spbas_allocate_data(result);
506: /* Reset the number of nonzeros in each row */
507: for (i = 0; i<nrows; i++) { result->row_nnz[i] = 0; }
509: /* Fill the data arrays */
510: if (in_matrix.values) {
511: for (i=0; i<ncols; i++) {
512: r_nnz = in_matrix.row_nnz[i];
513: irow = in_matrix.icols[i];
514: val = in_matrix.values[i];
516: if (col_idx_type == SPBAS_COLUMN_NUMBERS) {icol0=0;}
517: else if (col_idx_type == SPBAS_DIAGONAL_OFFSETS) {icol0=i;}
518: else if (col_idx_type == SPBAS_OFFSET_ARRAY) {icol0=in_matrix.icol0[i];}
519: for (j=0; j<r_nnz; j++) {
520: k = icol0 + irow[j];
521: result->icols[k][result->row_nnz[k]] = i;
522: result->values[k][result->row_nnz[k]] = val[j];
523: result->row_nnz[k]++;
524: }
525: }
526: } else {
527: for (i=0; i<ncols; i++) {
528: r_nnz = in_matrix.row_nnz[i];
529: irow = in_matrix.icols[i];
530:
531: if (col_idx_type == SPBAS_COLUMN_NUMBERS) {icol0=0;}
532: else if (col_idx_type == SPBAS_DIAGONAL_OFFSETS) {icol0=i;}
533: else if (col_idx_type == SPBAS_OFFSET_ARRAY) {icol0=in_matrix.icol0[i];}
535: for (j=0; j<r_nnz; j++) {
536: k = icol0 + irow[j];
537: result->icols[k][result->row_nnz[k]] = i;
538: result->row_nnz[k]++;
539: }
540: }
541: }
542: return(0);
543: }
545: /*
546: spbas_mergesort
548: mergesort for an array of intergers and an array of associated
549: reals
551: on output, icol[0..nnz-1] is increasing;
552: val[0..nnz-1] has undergone the same permutation as icol
554: NB: val may be PETSC_NULL: in that case, only the integers are sorted
556: */
559: PetscErrorCode spbas_mergesort(PetscInt nnz, PetscInt *icol, PetscScalar *val)
560: {
561: PetscInt istep; /* Chunk-sizes of already sorted parts of arrays */
562: PetscInt i, i1, i2; /* Loop counters for (partly) sorted arrays */
563: PetscInt istart, i1end, i2end; /* start of newly sorted array part, end of both parts */
564: PetscInt *ialloc; /* Allocated arrays */
565: PetscScalar *valloc=PETSC_NULL;
566: PetscInt *iswap; /* auxiliary pointers for swapping */
567: PetscScalar *vswap;
568: PetscInt *ihlp1; /* Pointers to new version of arrays, */
569: PetscScalar *vhlp1=PETSC_NULL; /* (arrays under construction) */
570: PetscInt *ihlp2; /* Pointers to previous version of arrays, */
571: PetscScalar *vhlp2=PETSC_NULL;
574: PetscMalloc(nnz*sizeof(PetscInt),&ialloc);
575: ihlp1 = ialloc;
576: ihlp2 = icol;
578: if (val) {
579: PetscMalloc(nnz*sizeof(PetscScalar),&valloc);
580: vhlp1 = valloc;
581: vhlp2 = val;
582: }
585: /* Sorted array chunks are first 1 long, and increase until they are the complete array */
586: for (istep=1; istep<nnz; istep*=2) {
587: /*
588: Combine sorted parts
589: istart:istart+istep-1 and istart+istep-1:istart+2*istep-1
590: of ihlp2 and vhlp2
591:
592: into one sorted part
593: istart:istart+2*istep-1
594: of ihlp1 and vhlp1
595: */
596: for (istart=0; istart<nnz; istart+=2*istep) {
598: /* Set counters and bound array part endings */
599: i1=istart; i1end = i1+istep; if (i1end>nnz) {i1end=nnz;}
600: i2=istart+istep; i2end = i2+istep; if (i2end>nnz) {i2end=nnz;}
602: /* Merge the two array parts */
603: if (val) {
604: for (i=istart; i<i2end; i++) {
605: if (i1<i1end && i2<i2end && ihlp2[i1] < ihlp2[i2]) {
606: ihlp1[i] = ihlp2[i1];
607: vhlp1[i] = vhlp2[i1];
608: i1++;
609: } else if (i2<i2end ){
610: ihlp1[i] = ihlp2[i2];
611: vhlp1[i] = vhlp2[i2];
612: i2++;
613: } else {
614: ihlp1[i] = ihlp2[i1];
615: vhlp1[i] = vhlp2[i1];
616: i1++;
617: }
618: }
619: } else {
620: for (i=istart; i<i2end; i++) {
621: if (i1<i1end && i2<i2end && ihlp2[i1] < ihlp2[i2]) {
622: ihlp1[i] = ihlp2[i1];
623: i1++;
624: } else if (i2<i2end ) {
625: ihlp1[i] = ihlp2[i2];
626: i2++;
627: } else {
628: ihlp1[i] = ihlp2[i1];
629: i1++;
630: }
631: }
632: }
633: }
635: /* Swap the two array sets */
636: iswap = ihlp2; ihlp2 = ihlp1; ihlp1 = iswap;
637: vswap = vhlp2; vhlp2 = vhlp1; vhlp1 = vswap;
638: }
640: /* Copy one more time in case the sorted arrays are the temporary ones */
641: if (ihlp2 != icol) {
642: for (i=0; i<nnz; i++) { icol[i] = ihlp2[i]; }
643: if (val) {
644: for (i=0; i<nnz; i++) { val[i] = vhlp2[i]; }
645: }
646: }
648: PetscFree(ialloc);
649: if(val){PetscFree(valloc);}
650: return(0);
651: }
653: /*
654: spbas_apply_reordering_rows:
655: apply the given reordering to the rows: matrix_A = matrix_A(perm,:);
656: */
659: PetscErrorCode spbas_apply_reordering_rows(spbas_matrix *matrix_A, const PetscInt *permutation)
660: {
661: PetscInt i,j,ip;
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=PETSC_NULL;
670: if (matrix_A->col_idx_type != SPBAS_DIAGONAL_OFFSETS) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP_SYS,"must have diagonal offsets in pattern\n");
672: if (do_values) {
673: PetscMalloc( sizeof(PetscScalar*)*nrows, &vals);
674: }
675: PetscMalloc( sizeof(PetscInt)*nrows, &row_nnz);
676: PetscMalloc( sizeof(PetscInt*)*nrows, &icols);
678: for (i=0; i<nrows;i++) {
679: ip = permutation[i];
680: if (do_values) {vals[i] = matrix_A->values[ip];}
681: icols[i] = matrix_A->icols[ip];
682: row_nnz[i] = matrix_A->row_nnz[ip];
683: for (j=0; j<row_nnz[i]; j++) { icols[i][j] += ip-i; }
684: }
686: if (do_values){ PetscFree(matrix_A->values);}
687: PetscFree(matrix_A->icols);
688: PetscFree(matrix_A->row_nnz);
690: if (do_values) { matrix_A->values = vals; }
691: matrix_A->icols = icols;
692: matrix_A->row_nnz = row_nnz;
693:
694: return(0);
695: }
698: /*
699: spbas_apply_reordering_cols:
700: apply the given reordering to the columns: matrix_A(:,perm) = matrix_A;
701: */
704: PetscErrorCode spbas_apply_reordering_cols( spbas_matrix *matrix_A,const PetscInt *permutation)
705: {
706: PetscInt i,j;
707: PetscInt nrows=matrix_A->nrows;
708: PetscInt row_nnz;
709: PetscInt *icols;
710: PetscBool do_values = matrix_A->values ? PETSC_TRUE : PETSC_FALSE;
711: PetscScalar *vals=PETSC_NULL;
715: if (matrix_A->col_idx_type != SPBAS_DIAGONAL_OFFSETS) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP_SYS, "must have diagonal offsets in pattern\n");
717: for (i=0; i<nrows;i++) {
718: icols = matrix_A->icols[i];
719: row_nnz = matrix_A->row_nnz[i];
720: if (do_values) { vals = matrix_A->values[i]; }
722: for (j=0; j<row_nnz; j++) {
723: icols[j] = permutation[i+icols[j]]-i;
724: }
725: spbas_mergesort(row_nnz, icols, vals);
726: }
727: return(0);
728: }
730: /*
731: spbas_apply_reordering:
732: apply the given reordering: matrix_A(perm,perm) = matrix_A;
733: */
736: PetscErrorCode spbas_apply_reordering( spbas_matrix *matrix_A, const PetscInt *permutation, const PetscInt * inv_perm)
737: {
740: spbas_apply_reordering_rows( matrix_A, inv_perm);
741: spbas_apply_reordering_cols( matrix_A, permutation);
742: return(0);
743: }
747: PetscErrorCode spbas_pattern_only(PetscInt nrows, PetscInt ncols, PetscInt *ai, PetscInt *aj, spbas_matrix * result)
748: {
749: spbas_matrix retval;
750: PetscInt i, j, i0, r_nnz;
755: /* Copy input values */
756: retval.nrows = nrows;
757: retval.ncols = ncols;
758: retval.nnz = ai[nrows];
760: retval.block_data = PETSC_TRUE;
761: retval.col_idx_type = SPBAS_DIAGONAL_OFFSETS;
763: /* Allocate output matrix */
764: spbas_allocate_pattern(&retval, PETSC_FALSE);
765: for (i=0; i<nrows; i++) {retval.row_nnz[i] = ai[i+1]-ai[i];}
766: spbas_allocate_data(&retval);
767: /* Copy the structure */
768: for (i = 0; i<retval.nrows; i++) {
769: i0 = ai[i];
770: r_nnz = ai[i+1]-i0;
772: for (j=0; j<r_nnz; j++) {
773: retval.icols[i][j] = aj[i0+j]-i;
774: }
775: }
776: *result = retval;
777: return(0);
778: }
781: /*
782: spbas_mark_row_power:
783: Mark the columns in row 'row' which are nonzero in
784: matrix^2log(marker).
785: */
788: PetscErrorCode spbas_mark_row_power(
789: PetscInt *iwork, /* marker-vector */
790: PetscInt row, /* row for which the columns are marked */
791: spbas_matrix * in_matrix, /* matrix for which the power is being calculated */
792: PetscInt marker, /* marker-value: 2^power */
793: PetscInt minmrk, /* lower bound for marked points */
794: PetscInt maxmrk) /* upper bound for marked points */
795: {
797: PetscInt i,j, nnz;
800: nnz = in_matrix->row_nnz[row];
802: /* For higher powers, call this function recursively */
803: if (marker>1) {
804: for (i=0; i<nnz;i++) {
805: j = row + in_matrix->icols[row][i];
806: if (minmrk<=j && j<maxmrk && iwork[j] < marker ) {
807: spbas_mark_row_power( iwork, row + in_matrix->icols[row][i],in_matrix, marker/2,minmrk,maxmrk);
808: iwork[j] |= marker;
809: }
810: }
811: } else {
812: /* Mark the columns reached */
813: for (i=0; i<nnz;i++) {
814: j = row + in_matrix->icols[row][i];
815: if (minmrk<=j && j<maxmrk) { iwork[j] |= 1; }
816: }
817: }
818: return(0);
819: }
822: /*
823: spbas_power
824: Calculate sparseness patterns for incomplete Cholesky decompositions
825: of a given order: (almost) all nonzeros of the matrix^(order+1) which
826: are inside the band width are found and stored in the output sparseness
827: pattern.
828: */
831: PetscErrorCode spbas_power (spbas_matrix in_matrix,PetscInt power, spbas_matrix * result)
832: {
833: spbas_matrix retval;
834: PetscInt nrows = in_matrix.nrows;
835: PetscInt ncols = in_matrix.ncols;
836: PetscInt i, j, kend;
837: PetscInt nnz, inz;
838: PetscInt *iwork;
839: PetscInt marker;
840: PetscInt maxmrk=0;
844: if (in_matrix.col_idx_type != SPBAS_DIAGONAL_OFFSETS) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP_SYS,"must have diagonal offsets in pattern\n");
845: if (ncols != nrows) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "Dimension error\n");
846: if (in_matrix.values) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "Input array must be sparseness pattern (no values)");
847: if (power<=0) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP_SYS, "Power must be 1 or up");
849: /* Copy input values*/
850: retval.nrows = ncols;
851: retval.ncols = nrows;
852: retval.nnz = 0;
853: retval.col_idx_type = SPBAS_DIAGONAL_OFFSETS;
854: retval.block_data = PETSC_FALSE;
856: /* Allocate sparseness pattern */
857: spbas_allocate_pattern(&retval, in_matrix.values ? PETSC_TRUE : PETSC_FALSE);
859: /* Allocate marker array */
860: PetscMalloc(nrows * sizeof(PetscInt), &iwork);
862: /* Erase the pattern for this row */
863: PetscMemzero( (void *) iwork, retval.nrows*sizeof(PetscInt));
865: /* Calculate marker values */
866: marker = 1; for (i=1; i<power; i++) {marker*=2;}
868: for (i=0; i<nrows; i++) {
869: /* Calculate the pattern for each row */
871: nnz = in_matrix.row_nnz[i];
872: kend = i+in_matrix.icols[i][nnz-1];
873: if (maxmrk<=kend) {maxmrk=kend+1;}
874: spbas_mark_row_power( iwork, i, &in_matrix, marker, i, maxmrk);
876: /* Count the columns*/
877: nnz = 0;
878: for (j=i; j<maxmrk; j++) {nnz+= (iwork[j]!=0);}
880: /* Allocate the column indices */
881: retval.row_nnz[i] = nnz;
882: PetscMalloc(nnz*sizeof(PetscInt),&retval.icols[i]);
884: /* Administrate the column indices */
885: inz = 0;
886: for (j=i; j<maxmrk; j++) {
887: if (iwork[j]) {
888: retval.icols[i][inz] = j-i;
889: inz++;
890: iwork[j]=0;
891: }
892: }
893: retval.nnz += nnz;
894: };
895: PetscFree(iwork);
896: *result = retval;
897: return(0);
898: }
902: /*
903: spbas_keep_upper:
904: remove the lower part of the matrix: keep the upper part
905: */
908: PetscErrorCode spbas_keep_upper( spbas_matrix * inout_matrix)
909: {
910: PetscInt i, j;
911: PetscInt jstart;
914: if (inout_matrix->block_data) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP_SYS, "Not yet for block data matrices\n");
915: for (i=0; i<inout_matrix->nrows; i++) {
916: for (jstart=0; (jstart<inout_matrix->row_nnz[i]) && (inout_matrix->icols[i][jstart]<0); jstart++){}
917: if (jstart>0) {
918: for (j=0; j<inout_matrix->row_nnz[i]-jstart; j++) {
919: inout_matrix->icols[i][j] = inout_matrix->icols[i][j+jstart];
920: }
922: if (inout_matrix->values) {
923: for (j=0; j<inout_matrix->row_nnz[i]-jstart; j++) {
924: inout_matrix->values[i][j] = inout_matrix->values[i][j+jstart];
925: }
926: }
928: inout_matrix->row_nnz[i] -= jstart;
930: inout_matrix->icols[i] = (PetscInt*) realloc((void*) inout_matrix->icols[i], inout_matrix->row_nnz[i]*sizeof(PetscInt));
932: if (inout_matrix->values) {
933: inout_matrix->values[i] = (PetscScalar*) realloc((void*) inout_matrix->values[i], inout_matrix->row_nnz[i]*sizeof(PetscScalar));
934: }
935: inout_matrix->nnz -= jstart;
936: }
937: }
938: return(0);
939: }