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