xref: /petsc/src/mat/impls/sell/seq/fdsell.c (revision f4f49eeac7efa77fffa46b7ff95a3ed169f659ed)
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 */
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   *nn = n;
20   if (!ia) PetscFunctionReturn(PETSC_SUCCESS);
21 
22   PetscCall(PetscCalloc1(n + 1, &collengths));
23   PetscCall(PetscMalloc1(n + 1, &cia));
24   PetscCall(PetscMalloc1(a->nz + 1, &cja));
25   PetscCall(PetscMalloc1(a->nz + 1, &cspidx));
26 
27   totalslices = A->rmap->n / 8 + ((A->rmap->n & 0x07) ? 1 : 0); /* floor(n/8) */
28   for (i = 0; i < totalslices; i++) {                           /* loop over slices */
29     for (j = a->sliidx[i], row = 0; j < a->sliidx[i + 1]; j++, row = ((row + 1) & 0x07)) {
30       isnonzero = (PetscBool)((j - a->sliidx[i]) / 8 < a->rlen[8 * i + row]);
31       if (isnonzero) collengths[a->colidx[j]]++;
32     }
33   }
34 
35   cia[0] = oshift;
36   for (i = 0; i < n; i++) cia[i + 1] = cia[i] + collengths[i];
37   PetscCall(PetscArrayzero(collengths, n));
38 
39   for (i = 0; i < totalslices; i++) { /* loop over slices */
40     for (j = a->sliidx[i], row = 0; j < a->sliidx[i + 1]; j++, row = ((row + 1) & 0x07)) {
41       isnonzero = (PetscBool)((j - a->sliidx[i]) / 8 < a->rlen[8 * i + row]);
42       if (isnonzero) {
43         col                                         = a->colidx[j];
44         cspidx[cia[col] + collengths[col] - oshift] = j;                    /* index of a->colidx */
45         cja[cia[col] + collengths[col] - oshift]    = 8 * i + row + oshift; /* row index */
46         collengths[col]++;
47       }
48     }
49   }
50 
51   PetscCall(PetscFree(collengths));
52   *ia    = cia;
53   *ja    = cja;
54   *spidx = cspidx;
55   PetscFunctionReturn(PETSC_SUCCESS);
56 }
57 
58 PetscErrorCode MatRestoreColumnIJ_SeqSELL_Color(Mat A, PetscInt oshift, PetscBool symmetric, PetscBool inodecompressed, PetscInt *n, const PetscInt *ia[], const PetscInt *ja[], PetscInt *spidx[], PetscBool *done)
59 {
60   PetscFunctionBegin;
61   if (!ia) PetscFunctionReturn(PETSC_SUCCESS);
62   PetscCall(PetscFree(*ia));
63   PetscCall(PetscFree(*ja));
64   PetscCall(PetscFree(*spidx));
65   PetscFunctionReturn(PETSC_SUCCESS);
66 }
67