xref: /petsc/src/mat/impls/aij/mpi/mpihashmat.h (revision 9680b7cedbd3b2fa4b943552cd701a6633a2ad02)

/*
   used by MPIAIJ, BAIJ and SBAIJ to reduce code duplication

     define TYPE to AIJ BAIJ or SBAIJ
            TYPE_SBAIJ for SBAIJ matrix

*/

static PetscErrorCode MatSetValues_MPI_Hash(Mat A, PetscInt m, const PetscInt *rows, PetscInt n, const PetscInt *cols, const PetscScalar *values, InsertMode addv)
{
  PetscConcat(Mat_MPI, TYPE) *a = (PetscConcat(Mat_MPI, TYPE) *)A->data;
  const PetscInt rStart         = A->rmap->rstart;
  const PetscInt rEnd           = A->rmap->rend;
  const PetscInt cStart         = A->cmap->rstart;
  const PetscInt cEnd           = A->cmap->rend;
#if defined(TYPE_SBAIJ)
  const PetscInt bs = A->rmap->bs;
#endif
  const PetscBool ignorezeroentries = ((Mat_SeqAIJ *)a->A->data)->ignorezeroentries;

  PetscFunctionBegin;
  for (PetscInt r = 0; r < m; ++r) {
    PetscScalar value;
    if (rows[r] < 0) continue;
    if (rows[r] < rStart || rows[r] >= rEnd) {
      PetscCheck(!A->nooffprocentries, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Setting off process row %" PetscInt_FMT " even though MatSetOption(,MAT_NO_OFF_PROC_ENTRIES,PETSC_TRUE) was set", rows[r]);
      if (!a->donotstash) {
        A->assembled = PETSC_FALSE;
        if (a->roworiented) {
          PetscCall(MatStashValuesRow_Private(&A->stash, rows[r], n, cols, PetscSafePointerPlusOffset(values, r * n), (PetscBool)(ignorezeroentries && (addv == ADD_VALUES))));
        } else {
          PetscCall(MatStashValuesCol_Private(&A->stash, rows[r], n, cols, PetscSafePointerPlusOffset(values, r), m, (PetscBool)(ignorezeroentries && (addv == ADD_VALUES))));
        }
      }
    } else {
      for (PetscInt c = 0; c < n; ++c) {
#if defined(TYPE_SBAIJ)
        if (cols[c] / bs < rows[r] / bs) continue;
#else
        if (cols[c] < 0) continue;
#endif
        value = values ? (a->roworiented ? values[r * n + c] : values[r + m * c]) : 0;
        if (ignorezeroentries && value == 0.0 && (addv == ADD_VALUES) && rows[r] != cols[c]) continue;
        if (cols[c] >= cStart && cols[c] < cEnd) {
          PetscCall(MatSetValue(a->A, rows[r] - rStart, cols[c] - cStart, value, addv));
        } else if (!ignorezeroentries || value != 0.0) {
          /* MPIAIJ never inserts in B if it is 0 */
          PetscCall(MatSetValue(a->B, rows[r] - rStart, cols[c], value, addv));
        }
      }
    }
  }
  PetscFunctionReturn(PETSC_SUCCESS);
}

static PetscErrorCode MatAssemblyBegin_MPI_Hash(Mat A, PETSC_UNUSED MatAssemblyType type)
{
  PetscConcat(Mat_MPI, TYPE) *a = (PetscConcat(Mat_MPI, TYPE) *)A->data;
  PetscInt nstash, reallocs;

  PetscFunctionBegin;
  if (a->donotstash || A->nooffprocentries) PetscFunctionReturn(PETSC_SUCCESS);
  PetscCall(MatStashScatterBegin_Private(A, &A->stash, A->rmap->range));
  PetscCall(MatStashGetInfo_Private(&A->stash, &nstash, &reallocs));
  PetscCall(PetscInfo(A, "Stash has %" PetscInt_FMT " entries, uses %" PetscInt_FMT " mallocs.\n", nstash, reallocs));
  PetscFunctionReturn(PETSC_SUCCESS);
}

static PetscErrorCode MatFinishScatterAndSetValues_MPI_Hash(Mat A)
{
  PetscConcat(Mat_MPI, TYPE) *a = (PetscConcat(Mat_MPI, TYPE) *)A->data;
  PetscMPIInt  n;
  PetscScalar *val;
  PetscInt    *row, *col;
  PetscInt     j, ncols, flg, rstart;

  PetscFunctionBegin;
  if (!a->donotstash && !A->nooffprocentries) {
    while (1) {
      PetscCall(MatStashScatterGetMesg_Private(&A->stash, &n, &row, &col, &val, &flg));
      if (!flg) break;

      for (PetscInt i = 0; i < n;) {
        /* Now identify the consecutive vals belonging to the same row */
        for (j = i, rstart = row[j]; j < n; j++) {
          if (row[j] != rstart) break;
        }
        if (j < n) ncols = j - i;
        else ncols = n - i;
        /* Now assemble all these values with a single function call */
        PetscCall(MatSetValues_MPI_Hash(A, 1, row + i, ncols, col + i, val + i, A->insertmode));
        i = j;
      }
    }
    PetscCall(MatStashScatterEnd_Private(&A->stash));
  }
  PetscFunctionReturn(PETSC_SUCCESS);
}

static PetscErrorCode MatAssemblyEnd_MPI_Hash(Mat A, MatAssemblyType type)
{
  PetscConcat(Mat_MPI, TYPE) *a = (PetscConcat(Mat_MPI, TYPE) *)A->data;

  PetscFunctionBegin;
  PetscCall(MatFinishScatterAndSetValues_MPI_Hash(A));
  if (type != MAT_FINAL_ASSEMBLY) PetscFunctionReturn(PETSC_SUCCESS);

  A->insertmode = NOT_SET_VALUES; /* this was set by the previous calls to MatSetValues() */

  A->ops[0]      = a->cops;
  A->hash_active = PETSC_FALSE;

  PetscCall(MatAssemblyBegin(a->A, MAT_FINAL_ASSEMBLY));
  PetscCall(MatAssemblyEnd(a->A, MAT_FINAL_ASSEMBLY));
  PetscCall(MatAssemblyBegin(a->B, MAT_FINAL_ASSEMBLY));
  PetscCall(MatAssemblyEnd(a->B, MAT_FINAL_ASSEMBLY));
  PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
  PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
  PetscFunctionReturn(PETSC_SUCCESS);
}

static PetscErrorCode MatCopyHashToXAIJ_MPI_Hash(Mat A, Mat B)
{
  PetscConcat(Mat_MPI, TYPE) *a = (PetscConcat(Mat_MPI, TYPE) *)A->data, *b = (PetscConcat(Mat_MPI, TYPE) *)B->data;

  PetscFunctionBegin;
  /* Let's figure there's no harm done in doing the scatters for A now even if A != B */
  PetscCall(MatAssemblyBegin_MPI_Hash(A, /*unused*/ MAT_FINAL_ASSEMBLY));
  PetscCall(MatFinishScatterAndSetValues_MPI_Hash(A));

  PetscCall(MatCopyHashToXAIJ(a->A, b->A));
  PetscCall(MatCopyHashToXAIJ(a->B, b->B));
  PetscCall(MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY));
  PetscCall(MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY));
  PetscFunctionReturn(PETSC_SUCCESS);
}

static PetscErrorCode MatDestroy_MPI_Hash(Mat A)
{
  PetscConcat(Mat_MPI, TYPE) *a = (PetscConcat(Mat_MPI, TYPE) *)A->data;

  PetscFunctionBegin;
  PetscCall(MatStashDestroy_Private(&A->stash));
  PetscCall(MatDestroy(&a->A));
  PetscCall(MatDestroy(&a->B));
  PetscCall((*a->cops.destroy)(A));
  PetscFunctionReturn(PETSC_SUCCESS);
}

static PetscErrorCode MatZeroEntries_MPI_Hash(PETSC_UNUSED Mat A)
{
  PetscFunctionBegin;
  PetscFunctionReturn(PETSC_SUCCESS);
}

static PetscErrorCode MatSetRandom_MPI_Hash(Mat A, PETSC_UNUSED PetscRandom r)
{
  SETERRQ(PetscObjectComm((PetscObject)A), PETSC_ERR_ARG_WRONGSTATE, "Must set preallocation first");
}

static PetscErrorCode MatSetUp_MPI_Hash(Mat A)
{
  PetscConcat(Mat_MPI, TYPE) *a = (PetscConcat(Mat_MPI, TYPE) *)A->data;
  PetscMPIInt size;
#if !defined(TYPE_AIJ)
  PetscInt bs;
#endif

  PetscFunctionBegin;
  PetscCall(PetscInfo(A, "Using hash-based MatSetValues() for MATMPISBAIJ because no preallocation provided\n"));
  PetscCall(PetscLayoutSetUp(A->rmap));
  PetscCall(PetscLayoutSetUp(A->cmap));
  if (A->rmap->bs < 1) A->rmap->bs = 1;
  if (A->cmap->bs < 1) A->cmap->bs = 1;
  PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)A), &size));

#if !defined(TYPE_AIJ)
  PetscCall(MatGetBlockSize(A, &bs));
  /* these values are set in MatMPISBAIJSetPreallocation() */
  a->bs2 = bs * bs;
  a->mbs = A->rmap->n / bs;
  a->nbs = A->cmap->n / bs;
  a->Mbs = A->rmap->N / bs;
  a->Nbs = A->cmap->N / bs;

  for (PetscInt i = 0; i <= a->size; i++) a->rangebs[i] = A->rmap->range[i] / bs;
  a->rstartbs = A->rmap->rstart / bs;
  a->rendbs   = A->rmap->rend / bs;
  a->cstartbs = A->cmap->rstart / bs;
  a->cendbs   = A->cmap->rend / bs;
  PetscCall(MatStashCreate_Private(PetscObjectComm((PetscObject)A), A->rmap->bs, &A->bstash));
#endif

  PetscCall(MatCreate(PETSC_COMM_SELF, &a->A));
  PetscCall(MatSetSizes(a->A, A->rmap->n, A->cmap->n, A->rmap->n, A->cmap->n));
  PetscCall(MatSetBlockSizesFromMats(a->A, A, A));
#if defined(SUB_TYPE_CUSPARSE)
  PetscCall(MatSetType(a->A, MATSEQAIJCUSPARSE));
#else
  PetscCall(MatSetType(a->A, PetscConcat(MATSEQ, TYPE)));
#endif
  PetscCall(MatSetUp(a->A));

  PetscCall(MatCreate(PETSC_COMM_SELF, &a->B));
  PetscCall(MatSetSizes(a->B, A->rmap->n, size > 1 ? A->cmap->N : 0, A->rmap->n, size > 1 ? A->cmap->N : 0));
  PetscCall(MatSetBlockSizesFromMats(a->B, A, A));
#if defined(TYPE_SBAIJ)
  PetscCall(MatSetType(a->B, MATSEQBAIJ));
#else
  #if defined(SUB_TYPE_CUSPARSE)
  PetscCall(MatSetType(a->B, MATSEQAIJCUSPARSE));
  #else
  PetscCall(MatSetType(a->B, PetscConcat(MATSEQ, TYPE)));
  #endif
#endif
  PetscCall(MatSetUp(a->B));

  /* keep a record of the operations so they can be reset when the hash handling is complete */
  a->cops                  = A->ops[0];
  A->ops->assemblybegin    = MatAssemblyBegin_MPI_Hash;
  A->ops->assemblyend      = MatAssemblyEnd_MPI_Hash;
  A->ops->setvalues        = MatSetValues_MPI_Hash;
  A->ops->destroy          = MatDestroy_MPI_Hash;
  A->ops->zeroentries      = MatZeroEntries_MPI_Hash;
  A->ops->setrandom        = MatSetRandom_MPI_Hash;
  A->ops->copyhashtoxaij   = MatCopyHashToXAIJ_MPI_Hash;
  A->ops->setvaluesblocked = NULL;

  A->preallocated = PETSC_TRUE;
  A->hash_active  = PETSC_TRUE;
  PetscFunctionReturn(PETSC_SUCCESS);
}