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:   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);

     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));

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

     36:   cia[0] = oshift;
     37:   for (i = 0; i < n; i++) cia[i + 1] = cia[i] + collengths[i];
     38:   PetscCall(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:   PetscCall(PetscFree(collengths));
     53:   *ia    = cia;
     54:   *ja    = cja;
     55:   *spidx = cspidx;
     56:   PetscFunctionReturn(PETSC_SUCCESS);
     57: }

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