1 #include <petsc/private/matimpl.h>
2
MatCreateFromISLocalToGlobalMapping(ISLocalToGlobalMapping lgmap,Mat A,PetscBool cols,PetscBool trans,MatType ptype,Mat * P)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