Actual source code: fdsell.c

  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;

 18:   *nn = n;
 19:   if (!ia) return 0;

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

 26:   totalslices = A->rmap->n/8+((A->rmap->n & 0x07)?1:0); /* floor(n/8) */
 27:   for (i=0; i<totalslices; i++) { /* loop over slices */
 28:     for (j=a->sliidx[i],row=0; j<a->sliidx[i+1]; j++,row=((row+1)&0x07)) {
 29:       isnonzero = (PetscBool)((j-a->sliidx[i])/8 < a->rlen[8*i+row]);
 30:       if (isnonzero) collengths[a->colidx[j]]++;
 31:     }
 32:   }

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

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

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

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

 61:   if (!ia) return 0;
 62:   PetscFree(*ia);
 63:   PetscFree(*ja);
 64:   PetscFree(*spidx);
 65:   return 0;
 66: }