1 #include <../src/mat/impls/sell/seq/sell.h>
2 #include <../src/mat/impls/aij/seq/aij.h>
3 #include <petsc/private/isimpl.h>
4
5 /*
6 MatGetColumnIJ_SeqSELL_Color() and MatRestoreColumnIJ_SeqSELL_Color() are customized from
7 MatGetColumnIJ_SeqSELL() and MatRestoreColumnIJ_SeqSELL() by adding an output
8 spidx[], index of a->a, to be used in MatTransposeColoringCreate_SeqSELL() and MatFDColoringCreate_SeqSELL()
9 */
MatGetColumnIJ_SeqSELL_Color(Mat A,PetscInt oshift,PetscBool symmetric,PetscBool inodecompressed,PetscInt * nn,const PetscInt * ia[],const PetscInt * ja[],PetscInt * spidx[],PetscBool * done)10 PetscErrorCode MatGetColumnIJ_SeqSELL_Color(Mat A, PetscInt oshift, PetscBool symmetric, PetscBool inodecompressed, PetscInt *nn, const PetscInt *ia[], const PetscInt *ja[], PetscInt *spidx[], PetscBool *done)
11 {
12 Mat_SeqSELL *a = (Mat_SeqSELL *)A->data;
13 PetscInt i, j, *collengths, *cia, *cja, n = A->cmap->n, totalslices;
14 PetscInt row, col;
15 PetscInt *cspidx;
16 PetscBool isnonzero;
17
18 PetscFunctionBegin;
19 PetscCheck(n >= 0, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Expected cmap->n %" PetscInt_FMT " >= 0", n);
20 *nn = n;
21 if (!ia) PetscFunctionReturn(PETSC_SUCCESS);
22
23 PetscCall(PetscCalloc1(n, &collengths));
24 PetscCall(PetscMalloc1(n + 1, &cia));
25 PetscCall(PetscMalloc1(a->nz + 1, &cja));
26 PetscCall(PetscMalloc1(a->nz + 1, &cspidx));
27
28 totalslices = PetscCeilInt(A->rmap->n, 8);
29 for (i = 0; i < totalslices; i++) { /* loop over slices */
30 for (j = a->sliidx[i], row = 0; j < a->sliidx[i + 1]; j++, row = ((row + 1) & 0x07)) {
31 isnonzero = (PetscBool)((j - a->sliidx[i]) / 8 < a->rlen[8 * i + row]);
32 if (isnonzero) collengths[a->colidx[j]]++;
33 }
34 }
35
36 cia[0] = oshift;
37 for (i = 0; i < n; i++) cia[i + 1] = cia[i] + collengths[i];
38 PetscCall(PetscArrayzero(collengths, n));
39
40 for (i = 0; i < totalslices; i++) { /* loop over slices */
41 for (j = a->sliidx[i], row = 0; j < a->sliidx[i + 1]; j++, row = ((row + 1) & 0x07)) {
42 isnonzero = (PetscBool)((j - a->sliidx[i]) / 8 < a->rlen[8 * i + row]);
43 if (isnonzero) {
44 col = a->colidx[j];
45 cspidx[cia[col] + collengths[col] - oshift] = j; /* index of a->colidx */
46 cja[cia[col] + collengths[col] - oshift] = 8 * i + row + oshift; /* row index */
47 collengths[col]++;
48 }
49 }
50 }
51
52 PetscCall(PetscFree(collengths));
53 *ia = cia;
54 *ja = cja;
55 *spidx = cspidx;
56 PetscFunctionReturn(PETSC_SUCCESS);
57 }
58
MatRestoreColumnIJ_SeqSELL_Color(Mat A,PetscInt oshift,PetscBool symmetric,PetscBool inodecompressed,PetscInt * n,const PetscInt * ia[],const PetscInt * ja[],PetscInt * spidx[],PetscBool * done)59 PetscErrorCode MatRestoreColumnIJ_SeqSELL_Color(Mat A, PetscInt oshift, PetscBool symmetric, PetscBool inodecompressed, PetscInt *n, const PetscInt *ia[], const PetscInt *ja[], PetscInt *spidx[], PetscBool *done)
60 {
61 PetscFunctionBegin;
62 if (!ia) PetscFunctionReturn(PETSC_SUCCESS);
63 PetscCall(PetscFree(*ia));
64 PetscCall(PetscFree(*ja));
65 PetscCall(PetscFree(*spidx));
66 PetscFunctionReturn(PETSC_SUCCESS);
67 }
68