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