1 #include <petsc/private/matimpl.h> 2 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; 8 9 PetscFunctionBegin; 10 PetscValidHeaderSpecific(lgmap, IS_LTOGM_CLASSID, 1); 11 PetscValidHeaderSpecific(A, MAT_CLASSID, 2); 12 PetscValidType(A, 2); 13 PetscValidLogicalCollectiveBool(A, cols, 3); 14 PetscValidLogicalCollectiveBool(A, trans, 4); 15 PetscAssertPointer(P, 6); 16 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; 32 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; 56 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 } 80