xref: /petsc/src/mat/impls/sell/seq/fdsell.c (revision 503c0ea9b45bcfbcebbb1ea5341243bbc69f0bea)
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(0);
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++) {
37     cia[i+1] = cia[i] + collengths[i];
38   }
39   PetscCall(PetscArrayzero(collengths,n));
40 
41   for (i=0; i<totalslices; i++) { /* loop over slices */
42     for (j=a->sliidx[i],row=0; j<a->sliidx[i+1]; j++,row=((row+1)&0x07)) {
43       isnonzero = (PetscBool)((j-a->sliidx[i])/8 < a->rlen[8*i+row]);
44       if (isnonzero) {
45         col = a->colidx[j];
46         cspidx[cia[col]+collengths[col]-oshift] = j; /* index of a->colidx */
47         cja[cia[col]+collengths[col]-oshift] = 8*i+row +oshift; /* row index */
48         collengths[col]++;
49       }
50     }
51   }
52 
53   PetscCall(PetscFree(collengths));
54   *ia    = cia; *ja = cja;
55   *spidx = cspidx;
56   PetscFunctionReturn(0);
57 }
58 
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 
63   if (!ia) PetscFunctionReturn(0);
64   PetscCall(PetscFree(*ia));
65   PetscCall(PetscFree(*ja));
66   PetscCall(PetscFree(*spidx));
67   PetscFunctionReturn(0);
68 }
69