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