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: }