1ae8d29abSPierre Jolivet 2ae8d29abSPierre Jolivet #include <../src/mat/impls/baij/mpi/mpibaij.h> /*I "petscmat.h" I*/ 3ae8d29abSPierre Jolivet #include <../src/mat/impls/aij/mpi/mpiaij.h> 4ae8d29abSPierre Jolivet #include <petsc/private/matimpl.h> 5ae8d29abSPierre Jolivet #include <petscmat.h> 6ae8d29abSPierre Jolivet 7ae8d29abSPierre Jolivet PETSC_INTERN PetscErrorCode MatConvert_MPIAIJ_MPIBAIJ(Mat A, MatType newtype,MatReuse reuse,Mat *newmat) 8ae8d29abSPierre Jolivet { 9ae8d29abSPierre Jolivet Mat M; 10ae8d29abSPierre Jolivet Mat_MPIAIJ *mpimat = (Mat_MPIAIJ*)A->data; 11ae8d29abSPierre Jolivet Mat_SeqAIJ *Aa = (Mat_SeqAIJ*)mpimat->A->data,*Ba = (Mat_SeqAIJ*)mpimat->B->data; 12ae8d29abSPierre Jolivet PetscInt *d_nnz,*o_nnz; 13ae8d29abSPierre Jolivet PetscInt i,m,n,lm,ln,bs=PetscAbs(A->rmap->bs); 14ae8d29abSPierre Jolivet 15ae8d29abSPierre Jolivet PetscFunctionBegin; 16ae8d29abSPierre Jolivet if (reuse != MAT_REUSE_MATRIX) { 17*9566063dSJacob Faibussowitsch PetscCall(MatGetSize(A,&m,&n)); 18*9566063dSJacob Faibussowitsch PetscCall(MatGetLocalSize(A,&lm,&ln)); 19*9566063dSJacob Faibussowitsch PetscCall(PetscMalloc2(lm/bs,&d_nnz,lm/bs,&o_nnz)); 20ae8d29abSPierre Jolivet 21ae8d29abSPierre Jolivet for (i=0; i<lm/bs; i++) { 22ae8d29abSPierre Jolivet d_nnz[i] = (Aa->i[i*bs+1] - Aa->i[i*bs])/bs; 23ae8d29abSPierre Jolivet o_nnz[i] = (Ba->i[i*bs+1] - Ba->i[i*bs])/bs; 24ae8d29abSPierre Jolivet } 25ae8d29abSPierre Jolivet 26*9566063dSJacob Faibussowitsch PetscCall(MatCreate(PetscObjectComm((PetscObject)A),&M)); 27*9566063dSJacob Faibussowitsch PetscCall(MatSetSizes(M,lm,ln,m,n)); 28*9566063dSJacob Faibussowitsch PetscCall(MatSetType(M,MATMPIBAIJ)); 29*9566063dSJacob Faibussowitsch PetscCall(MatSeqBAIJSetPreallocation(M,bs,0,d_nnz)); 30*9566063dSJacob Faibussowitsch PetscCall(MatMPIBAIJSetPreallocation(M,bs,0,d_nnz,0,o_nnz)); 31*9566063dSJacob Faibussowitsch PetscCall(PetscFree2(d_nnz,o_nnz)); 32ae8d29abSPierre Jolivet } else M = *newmat; 33ae8d29abSPierre Jolivet 34ae8d29abSPierre Jolivet /* reuse may not be equal to MAT_REUSE_MATRIX, but the basic converter will reallocate or replace newmat if this value is not used */ 35ae8d29abSPierre Jolivet /* if reuse is equal to MAT_INITIAL_MATRIX, it has been appropriately preallocated before */ 36ae8d29abSPierre Jolivet /* MAT_INPLACE_MATRIX, it will be replaced with MatHeaderReplace below */ 37*9566063dSJacob Faibussowitsch PetscCall(MatConvert_Basic(A,newtype,MAT_REUSE_MATRIX,&M)); 38ae8d29abSPierre Jolivet 39ae8d29abSPierre Jolivet if (reuse == MAT_INPLACE_MATRIX) { 40*9566063dSJacob Faibussowitsch PetscCall(MatHeaderReplace(A,&M)); 41ae8d29abSPierre Jolivet } else *newmat = M; 42ae8d29abSPierre Jolivet PetscFunctionReturn(0); 43ae8d29abSPierre Jolivet } 44