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_MPIAIJ_MPISBAIJ(Mat A, MatType newtype, MatReuse reuse, Mat *newmat)
  8: {
  9:   Mat                M;
 10:   Mat_MPIAIJ        *mpimat = (Mat_MPIAIJ *)A->data;
 11:   Mat_SeqAIJ        *Aa = (Mat_SeqAIJ *)mpimat->A->data, *Ba = (Mat_SeqAIJ *)mpimat->B->data;
 12:   PetscInt          *d_nnz, *o_nnz;
 13:   PetscInt           i, j, nz;
 14:   PetscInt           m, n, lm, ln;
 15:   PetscInt           rstart, rend, bs = PetscAbs(A->rmap->bs);
 16:   const PetscScalar *vwork;
 17:   const PetscInt    *cwork;

 20:   if (reuse != MAT_REUSE_MATRIX) {
 21:     MatGetSize(A, &m, &n);
 22:     MatGetLocalSize(A, &lm, &ln);
 23:     PetscMalloc2(lm / bs, &d_nnz, lm / bs, &o_nnz);

 25:     MatMarkDiagonal_SeqAIJ(mpimat->A);
 26:     for (i = 0; i < lm / bs; i++) {
 27:       if (Aa->i[i * bs + 1] == Aa->diag[i * bs]) { /* misses diagonal entry */
 28:         d_nnz[i] = (Aa->i[i * bs + 1] - Aa->i[i * bs]) / bs;
 29:       } else {
 30:         d_nnz[i] = (Aa->i[i * bs + 1] - Aa->diag[i * bs]) / bs;
 31:       }
 32:       o_nnz[i] = (Ba->i[i * bs + 1] - Ba->i[i * bs]) / bs;
 33:     }

 35:     MatCreate(PetscObjectComm((PetscObject)A), &M);
 36:     MatSetSizes(M, lm, ln, m, n);
 37:     MatSetType(M, MATMPISBAIJ);
 38:     MatSeqSBAIJSetPreallocation(M, bs, 0, d_nnz);
 39:     MatMPISBAIJSetPreallocation(M, bs, 0, d_nnz, 0, o_nnz);
 40:     PetscFree2(d_nnz, o_nnz);
 41:   } else M = *newmat;

 43:   if (bs == 1) {
 44:     MatGetOwnershipRange(A, &rstart, &rend);
 45:     for (i = rstart; i < rend; i++) {
 46:       MatGetRow(A, i, &nz, &cwork, &vwork);
 47:       if (nz) {
 48:         j = 0;
 49:         while (cwork[j] < i) {
 50:           j++;
 51:           nz--;
 52:         }
 53:         MatSetValues(M, 1, &i, nz, cwork + j, vwork + j, INSERT_VALUES);
 54:       }
 55:       MatRestoreRow(A, i, &nz, &cwork, &vwork);
 56:     }
 57:     MatAssemblyBegin(M, MAT_FINAL_ASSEMBLY);
 58:     MatAssemblyEnd(M, MAT_FINAL_ASSEMBLY);
 59:   } else {
 60:     MatSetOption(M, MAT_IGNORE_LOWER_TRIANGULAR, PETSC_TRUE);
 61:     /* reuse may not be equal to MAT_REUSE_MATRIX, but the basic converter will reallocate or replace newmat if this value is not used */
 62:     /* if reuse is equal to MAT_INITIAL_MATRIX, it has been appropriately preallocated before                                          */
 63:     /*                      MAT_INPLACE_MATRIX, it will be replaced with MatHeaderReplace below                                        */
 64:     MatConvert_Basic(A, newtype, MAT_REUSE_MATRIX, &M);
 65:   }

 67:   if (reuse == MAT_INPLACE_MATRIX) {
 68:     MatHeaderReplace(A, &M);
 69:   } else *newmat = M;
 70:   return 0;
 71: }
 72: /* contributed by Dahai Guo <dhguo@ncsa.uiuc.edu> April 2011 */
 73: PETSC_INTERN PetscErrorCode MatConvert_MPIBAIJ_MPISBAIJ(Mat A, MatType newtype, MatReuse reuse, Mat *newmat)
 74: {
 75:   Mat                M;
 76:   Mat_MPIBAIJ       *mpimat = (Mat_MPIBAIJ *)A->data;
 77:   Mat_SeqBAIJ       *Aa = (Mat_SeqBAIJ *)mpimat->A->data, *Ba = (Mat_SeqBAIJ *)mpimat->B->data;
 78:   PetscInt          *d_nnz, *o_nnz;
 79:   PetscInt           i, nz;
 80:   PetscInt           m, n, lm, ln;
 81:   PetscInt           rstart, rend;
 82:   const PetscScalar *vwork;
 83:   const PetscInt    *cwork;
 84:   PetscInt           bs = A->rmap->bs;

 86:   if (reuse != MAT_REUSE_MATRIX) {
 87:     MatGetSize(A, &m, &n);
 88:     MatGetLocalSize(A, &lm, &ln);
 89:     PetscMalloc2(lm / bs, &d_nnz, lm / bs, &o_nnz);

 91:     MatMarkDiagonal_SeqBAIJ(mpimat->A);
 92:     for (i = 0; i < lm / bs; i++) {
 93:       d_nnz[i] = Aa->i[i + 1] - Aa->diag[i];
 94:       o_nnz[i] = Ba->i[i + 1] - Ba->i[i];
 95:     }

 97:     MatCreate(PetscObjectComm((PetscObject)A), &M);
 98:     MatSetSizes(M, lm, ln, m, n);
 99:     MatSetType(M, MATMPISBAIJ);
100:     MatSeqSBAIJSetPreallocation(M, bs, 0, d_nnz);
101:     MatMPISBAIJSetPreallocation(M, bs, 0, d_nnz, 0, o_nnz);

103:     PetscFree2(d_nnz, o_nnz);
104:   } else M = *newmat;

106:   MatGetOwnershipRange(A, &rstart, &rend);
107:   MatSetOption(M, MAT_IGNORE_LOWER_TRIANGULAR, PETSC_TRUE);
108:   for (i = rstart; i < rend; i++) {
109:     MatGetRow(A, i, &nz, &cwork, &vwork);
110:     MatSetValues(M, 1, &i, nz, cwork, vwork, INSERT_VALUES);
111:     MatRestoreRow(A, i, &nz, &cwork, &vwork);
112:   }
113:   MatAssemblyBegin(M, MAT_FINAL_ASSEMBLY);
114:   MatAssemblyEnd(M, MAT_FINAL_ASSEMBLY);

116:   if (reuse == MAT_INPLACE_MATRIX) {
117:     MatHeaderReplace(A, &M);
118:   } else *newmat = M;
119:   return 0;
120: }