Actual source code: isltog.c

  1: #include <petsc/private/matimpl.h>

  3: PetscErrorCode MatCreateFromISLocalToGlobalMapping(ISLocalToGlobalMapping lgmap, Mat A, PetscBool cols, PetscBool trans, MatType ptype, Mat *P)
  4: {
  5:   PetscBool       matfree = PETSC_FALSE;
  6:   const PetscInt *idxs;
  7:   PetscInt        msize, *pidxs, c = 0;

  9:   PetscFunctionBegin;
 15:   PetscAssertPointer(P, 6);

 17:   if (!ptype) PetscCall(MatGetType(A, &ptype));
 18:   PetscCall(PetscStrcmpAny(ptype, &matfree, MATSHELL, MATSCATTER, ""));
 19:   PetscCall(ISLocalToGlobalMappingGetIndices(lgmap, &idxs));
 20:   PetscCall(ISLocalToGlobalMappingGetSize(lgmap, &msize));
 21:   PetscCall(PetscMalloc1(msize, &pidxs));
 22:   for (PetscInt i = 0; i < msize; i++)
 23:     if (idxs[i] >= 0) pidxs[c++] = idxs[i];
 24:   PetscCall(ISLocalToGlobalMappingRestoreIndices(lgmap, &idxs));
 25:   msize = c;
 26:   if (matfree) {
 27:     Vec        v, lv;
 28:     VecType    vtype;
 29:     IS         is;
 30:     VecScatter sct;
 31:     PetscBool  matshell;

 33:     if (cols) PetscCall(MatCreateVecs(A, &v, NULL));
 34:     else PetscCall(MatCreateVecs(A, NULL, &v));
 35:     PetscCall(VecGetType(v, &vtype));
 36:     PetscCall(VecCreate(PETSC_COMM_SELF, &lv));
 37:     PetscCall(VecSetSizes(lv, msize, PETSC_DECIDE));
 38:     PetscCall(VecSetType(lv, vtype));
 39:     PetscCall(ISCreateGeneral(PetscObjectComm((PetscObject)A), msize, pidxs, PETSC_USE_POINTER, &is));
 40:     if (trans) PetscCall(VecScatterCreate(lv, NULL, v, is, &sct));
 41:     else PetscCall(VecScatterCreate(v, is, lv, NULL, &sct));
 42:     PetscCall(MatCreateScatter(PetscObjectComm((PetscObject)A), sct, P));
 43:     PetscCall(PetscStrcmp(ptype, MATSHELL, &matshell));
 44:     if (matshell) {
 45:       Mat tP;
 46:       PetscCall(MatConvert(*P, ptype, MAT_INITIAL_MATRIX, &tP));
 47:       PetscCall(MatDestroy(P));
 48:       *P = tP;
 49:     }
 50:     PetscCall(ISDestroy(&is));
 51:     PetscCall(VecScatterDestroy(&sct));
 52:     PetscCall(VecDestroy(&lv));
 53:     PetscCall(VecDestroy(&v));
 54:   } else {
 55:     PetscInt lar, lac, rst;

 57:     PetscCall(MatGetLocalSize(A, &lar, &lac));
 58:     PetscCall(MatCreate(PetscObjectComm((PetscObject)A), P));
 59:     PetscCall(MatSetType(*P, ptype));
 60:     PetscCall(MatSetSizes(*P, msize, cols ? lac : lar, PETSC_DECIDE, PETSC_DECIDE));
 61:     PetscCall(MatSeqAIJSetPreallocation(*P, 1, NULL));
 62:     PetscCall(MatMPIAIJSetPreallocation(*P, 1, NULL, 1, NULL));
 63: #if defined(PETSC_HAVE_HYPRE)
 64:     PetscCall(MatHYPRESetPreallocation(*P, 1, NULL, 1, NULL));
 65: #endif
 66:     PetscCall(MatGetOwnershipRange(*P, &rst, NULL));
 67:     for (PetscInt i = 0; i < msize; i++) { PetscCall(MatSetValue(*P, i + rst, pidxs[i], 1.0, INSERT_VALUES)); }
 68:     PetscCall(MatAssemblyBegin(*P, MAT_FINAL_ASSEMBLY));
 69:     PetscCall(MatAssemblyEnd(*P, MAT_FINAL_ASSEMBLY));
 70:     if (trans) {
 71:       Mat tP;
 72:       PetscCall(MatTranspose(*P, MAT_INITIAL_MATRIX, &tP));
 73:       PetscCall(MatDestroy(P));
 74:       *P = tP;
 75:     }
 76:   }
 77:   PetscCall(PetscFree(pidxs));
 78:   PetscFunctionReturn(PETSC_SUCCESS);
 79: }