Actual source code: mpiaijsbaij.c
2: #include <../src/mat/impls/sbaij/mpi/mpisbaij.h>
3: #include <../src/mat/impls/aij/mpi/mpiaij.h>
4: #include <petsc/private/matimpl.h>
5: #include <petscmat.h>
7: PETSC_INTERN PetscErrorCode MatConvert_SeqAIJ_SeqBAIJ_Preallocate(Mat, PetscInt **);
8: PETSC_INTERN PetscErrorCode MatConvert_SeqAIJ_SeqSBAIJ_Preallocate(Mat, PetscInt **);
10: /* The code is virtually identical to MatConvert_MPIAIJ_MPIBAIJ() */
11: PETSC_INTERN PetscErrorCode MatConvert_MPIAIJ_MPISBAIJ(Mat A, MatType newtype, MatReuse reuse, Mat *newmat)
12: {
13: Mat M;
14: Mat_MPIAIJ *mpimat = (Mat_MPIAIJ *)A->data;
15: PetscInt *d_nnz, *o_nnz;
16: PetscInt m, n, lm, ln, bs = PetscAbs(A->rmap->bs);
18: if (reuse != MAT_REUSE_MATRIX) {
19: MatDisAssemble_MPIAIJ(A);
20: MatGetSize(A, &m, &n);
21: MatGetLocalSize(A, &lm, &ln);
22: MatConvert_SeqAIJ_SeqSBAIJ_Preallocate(mpimat->A, &d_nnz);
23: MatConvert_SeqAIJ_SeqBAIJ_Preallocate(mpimat->B, &o_nnz);
24: MatCreate(PetscObjectComm((PetscObject)A), &M);
25: MatSetSizes(M, lm, ln, m, n);
26: MatSetType(M, MATMPISBAIJ);
27: MatSeqSBAIJSetPreallocation(M, bs, 0, d_nnz);
28: MatMPISBAIJSetPreallocation(M, bs, 0, d_nnz, 0, o_nnz);
29: PetscFree(d_nnz);
30: PetscFree(o_nnz);
31: MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY);
32: 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: MatConvert_Basic(A, newtype, MAT_REUSE_MATRIX, &M);
40: if (reuse == MAT_INPLACE_MATRIX) {
41: MatHeaderReplace(A, &M);
42: } else *newmat = M;
43: return 0;
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: if (reuse != MAT_REUSE_MATRIX) {
61: MatGetSize(A, &m, &n);
62: MatGetLocalSize(A, &lm, &ln);
63: PetscMalloc2(lm / bs, &d_nnz, lm / bs, &o_nnz);
65: MatMarkDiagonal_SeqBAIJ(mpimat->A);
66: for (i = 0; i < lm / bs; i++) {
67: d_nnz[i] = Aa->i[i + 1] - Aa->diag[i];
68: o_nnz[i] = Ba->i[i + 1] - Ba->i[i];
69: }
71: MatCreate(PetscObjectComm((PetscObject)A), &M);
72: MatSetSizes(M, lm, ln, m, n);
73: MatSetType(M, MATMPISBAIJ);
74: MatSeqSBAIJSetPreallocation(M, bs, 0, d_nnz);
75: MatMPISBAIJSetPreallocation(M, bs, 0, d_nnz, 0, o_nnz);
77: PetscFree2(d_nnz, o_nnz);
78: } else M = *newmat;
80: MatGetOwnershipRange(A, &rstart, &rend);
81: MatSetOption(M, MAT_IGNORE_LOWER_TRIANGULAR, PETSC_TRUE);
82: for (i = rstart; i < rend; i++) {
83: MatGetRow(A, i, &nz, &cwork, &vwork);
84: MatSetValues(M, 1, &i, nz, cwork, vwork, INSERT_VALUES);
85: MatRestoreRow(A, i, &nz, &cwork, &vwork);
86: }
87: MatAssemblyBegin(M, MAT_FINAL_ASSEMBLY);
88: MatAssemblyEnd(M, MAT_FINAL_ASSEMBLY);
90: if (reuse == MAT_INPLACE_MATRIX) {
91: MatHeaderReplace(A, &M);
92: } else *newmat = M;
93: return 0;
94: }