1 #include <../src/mat/impls/sbaij/mpi/mpisbaij.h> /*I "petscmat.h" I*/
2 #include <../src/mat/impls/aij/mpi/mpiaij.h>
3 #include <petsc/private/matimpl.h>
4 #include <petscmat.h>
5
6 PETSC_INTERN PetscErrorCode MatConvert_SeqAIJ_SeqBAIJ_Preallocate(Mat, PetscInt **);
7 PETSC_INTERN PetscErrorCode MatConvert_SeqAIJ_SeqSBAIJ_Preallocate(Mat, PetscInt **);
8
9 /* The code is virtually identical to MatConvert_MPIAIJ_MPIBAIJ() */
MatConvert_MPIAIJ_MPISBAIJ(Mat A,MatType newtype,MatReuse reuse,Mat * newmat)10 PETSC_INTERN PetscErrorCode MatConvert_MPIAIJ_MPISBAIJ(Mat A, MatType newtype, MatReuse reuse, Mat *newmat)
11 {
12 Mat M;
13 Mat_MPIAIJ *mpimat = (Mat_MPIAIJ *)A->data;
14 PetscInt *d_nnz, *o_nnz;
15 PetscInt m, n, lm, ln, bs = A->rmap->bs;
16
17 PetscFunctionBegin;
18 if (reuse != MAT_REUSE_MATRIX) {
19 const PetscBool3 symmetric = A->symmetric, hermitian = A->hermitian, spd = A->spd;
20
21 PetscCall(MatDisAssemble_MPIAIJ(A, PETSC_FALSE));
22 PetscCall(MatGetSize(A, &m, &n));
23 PetscCall(MatGetLocalSize(A, &lm, &ln));
24 PetscCall(MatConvert_SeqAIJ_SeqSBAIJ_Preallocate(mpimat->A, &d_nnz));
25 PetscCall(MatConvert_SeqAIJ_SeqBAIJ_Preallocate(mpimat->B, &o_nnz));
26 PetscCall(MatCreate(PetscObjectComm((PetscObject)A), &M));
27 PetscCall(MatSetSizes(M, lm, ln, m, n));
28 PetscCall(MatSetType(M, MATMPISBAIJ));
29 PetscCall(MatSeqSBAIJSetPreallocation(M, bs, 0, d_nnz));
30 PetscCall(MatMPISBAIJSetPreallocation(M, bs, 0, d_nnz, 0, o_nnz));
31 PetscCall(PetscFree(d_nnz));
32 PetscCall(PetscFree(o_nnz));
33 PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
34 PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
35 if (symmetric != PETSC_BOOL3_UNKNOWN) PetscCall(MatSetOption(A, MAT_SYMMETRIC, PetscBool3ToBool(symmetric)));
36 if (hermitian != PETSC_BOOL3_UNKNOWN) PetscCall(MatSetOption(A, MAT_HERMITIAN, PetscBool3ToBool(hermitian)));
37 if (spd != PETSC_BOOL3_UNKNOWN) PetscCall(MatSetOption(A, MAT_SPD, PetscBool3ToBool(spd)));
38 } else M = *newmat;
39
40 /* reuse may not be equal to MAT_REUSE_MATRIX, but the basic converter will reallocate or replace newmat if this value is not used */
41 /* if reuse is equal to MAT_INITIAL_MATRIX, it has been appropriately preallocated before */
42 /* MAT_INPLACE_MATRIX, it will be replaced with MatHeaderReplace below */
43 PetscCall(MatConvert_Basic(A, newtype, MAT_REUSE_MATRIX, &M));
44
45 if (reuse == MAT_INPLACE_MATRIX) PetscCall(MatHeaderReplace(A, &M));
46 else *newmat = M;
47 PetscFunctionReturn(PETSC_SUCCESS);
48 }
49
50 /* contributed by Dahai Guo <dhguo@ncsa.uiuc.edu> April 2011 */
MatConvert_MPIBAIJ_MPISBAIJ(Mat A,MatType newtype,MatReuse reuse,Mat * newmat)51 PETSC_INTERN PetscErrorCode MatConvert_MPIBAIJ_MPISBAIJ(Mat A, MatType newtype, MatReuse reuse, Mat *newmat)
52 {
53 Mat M;
54 Mat_MPIBAIJ *mpimat = (Mat_MPIBAIJ *)A->data;
55 Mat_SeqBAIJ *Aa = (Mat_SeqBAIJ *)mpimat->A->data, *Ba = (Mat_SeqBAIJ *)mpimat->B->data;
56 PetscInt *d_nnz, *o_nnz;
57 PetscInt i, nz;
58 PetscInt m, n, lm, ln;
59 PetscInt rstart, rend;
60 const PetscScalar *vwork;
61 const PetscInt *cwork;
62 PetscInt bs = A->rmap->bs;
63 const PetscInt *adiag;
64
65 PetscFunctionBegin;
66 if (reuse != MAT_REUSE_MATRIX) {
67 PetscCall(MatGetSize(A, &m, &n));
68 PetscCall(MatGetLocalSize(A, &lm, &ln));
69 PetscCall(PetscMalloc2(lm / bs, &d_nnz, lm / bs, &o_nnz));
70
71 PetscCall(MatGetDiagonalMarkers_SeqBAIJ(mpimat->A, &adiag, NULL));
72 for (i = 0; i < lm / bs; i++) {
73 d_nnz[i] = Aa->i[i + 1] - adiag[i];
74 o_nnz[i] = Ba->i[i + 1] - Ba->i[i];
75 }
76
77 PetscCall(MatCreate(PetscObjectComm((PetscObject)A), &M));
78 PetscCall(MatSetSizes(M, lm, ln, m, n));
79 PetscCall(MatSetType(M, MATMPISBAIJ));
80 PetscCall(MatSeqSBAIJSetPreallocation(M, bs, 0, d_nnz));
81 PetscCall(MatMPISBAIJSetPreallocation(M, bs, 0, d_nnz, 0, o_nnz));
82
83 PetscCall(PetscFree2(d_nnz, o_nnz));
84 } else M = *newmat;
85
86 PetscCall(MatGetOwnershipRange(A, &rstart, &rend));
87 PetscCall(MatSetOption(M, MAT_IGNORE_LOWER_TRIANGULAR, PETSC_TRUE));
88 for (i = rstart; i < rend; i++) {
89 PetscCall(MatGetRow(A, i, &nz, &cwork, &vwork));
90 PetscCall(MatSetValues(M, 1, &i, nz, cwork, vwork, INSERT_VALUES));
91 PetscCall(MatRestoreRow(A, i, &nz, &cwork, &vwork));
92 }
93 PetscCall(MatAssemblyBegin(M, MAT_FINAL_ASSEMBLY));
94 PetscCall(MatAssemblyEnd(M, MAT_FINAL_ASSEMBLY));
95
96 if (reuse == MAT_INPLACE_MATRIX) PetscCall(MatHeaderReplace(A, &M));
97 else *newmat = M;
98 PetscFunctionReturn(PETSC_SUCCESS);
99 }
100