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