1: #include <../src/mat/impls/sell/seq/sell.h>
  2: #include <../src/mat/impls/aij/seq/aij.h>
  3: #include <petsc/private/isimpl.h>

  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;

 20:   *nn = n;
 21:   if (!ia) return(0);

 23:   PetscCalloc1(n+1,&collengths);
 24:   PetscMalloc1(n+1,&cia);
 25:   PetscMalloc1(a->nz+1,&cja);
 26:   PetscMalloc1(a->nz+1,&cspidx);

 28:   totalslices = A->rmap->n/8+((A->rmap->n & 0x07)?1:0); /* floor(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:   }

 36:   cia[0] = oshift;
 37:   for (i=0; i<n; i++) {
 38:     cia[i+1] = cia[i] + collengths[i];
 39:   }
 40:   PetscArrayzero(collengths,n);

 42:   for (i=0; i<totalslices; i++) { /* loop over slices */
 43:     for (j=a->sliidx[i],row=0; j<a->sliidx[i+1]; j++,row=((row+1)&0x07)) {
 44:       isnonzero = (PetscBool)((j-a->sliidx[i])/8 < a->rlen[8*i+row]);
 45:       if (isnonzero) {
 46:         col = a->colidx[j];
 47:         cspidx[cia[col]+collengths[col]-oshift] = j; /* index of a->colidx */
 48:         cja[cia[col]+collengths[col]-oshift] = 8*i+row +oshift; /* row index */
 49:         collengths[col]++;
 50:       }
 51:     }
 52:   }

 54:   PetscFree(collengths);
 55:   *ia    = cia; *ja = cja;
 56:   *spidx = cspidx;
 57:   return(0);
 58: }

 60: PetscErrorCode MatRestoreColumnIJ_SeqSELL_Color(Mat A,PetscInt oshift,PetscBool symmetric,PetscBool inodecompressed,PetscInt *n,const PetscInt *ia[],const PetscInt *ja[],PetscInt *spidx[],PetscBool *done)
 61: {


 66:   if (!ia) return(0);
 67:   PetscFree(*ia);
 68:   PetscFree(*ja);
 69:   PetscFree(*spidx);
 70:   return(0);
 71: }