1 2 #include <../src/mat/impls/sbaij/mpi/mpisbaij.h> /*I "petscmat.h" I*/ 3 #include <../src/mat/impls/aij/mpi/mpiaij.h> 4 #include <petsc/private/matimpl.h> 5 #include <petscmat.h> 6 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; 18 19 PetscFunctionBegin; 20 PetscCheck(A->symmetric == PETSC_BOOL3_TRUE || A->hermitian == PETSC_BOOL3_TRUE,PetscObjectComm((PetscObject)A),PETSC_ERR_USER,"Matrix must be symmetric or Hermitian. Call MatSetOption(mat,MAT_SYMMETRIC,PETSC_TRUE) or MatSetOption(mat,MAT_HERMITIAN,PETSC_TRUE)"); 21 if (reuse != MAT_REUSE_MATRIX) { 22 PetscCall(MatGetSize(A,&m,&n)); 23 PetscCall(MatGetLocalSize(A,&lm,&ln)); 24 PetscCall(PetscMalloc2(lm/bs,&d_nnz,lm/bs,&o_nnz)); 25 26 PetscCall(MatMarkDiagonal_SeqAIJ(mpimat->A)); 27 for (i=0; i<lm/bs; i++) { 28 if (Aa->i[i*bs+1] == Aa->diag[i*bs]) { /* misses diagonal entry */ 29 d_nnz[i] = (Aa->i[i*bs+1] - Aa->i[i*bs])/bs; 30 } else { 31 d_nnz[i] = (Aa->i[i*bs+1] - Aa->diag[i*bs])/bs; 32 } 33 o_nnz[i] = (Ba->i[i*bs+1] - Ba->i[i*bs])/bs; 34 } 35 36 PetscCall(MatCreate(PetscObjectComm((PetscObject)A),&M)); 37 PetscCall(MatSetSizes(M,lm,ln,m,n)); 38 PetscCall(MatSetType(M,MATMPISBAIJ)); 39 PetscCall(MatSeqSBAIJSetPreallocation(M,bs,0,d_nnz)); 40 PetscCall(MatMPISBAIJSetPreallocation(M,bs,0,d_nnz,0,o_nnz)); 41 PetscCall(PetscFree2(d_nnz,o_nnz)); 42 } else M = *newmat; 43 44 if (bs == 1) { 45 PetscCall(MatGetOwnershipRange(A,&rstart,&rend)); 46 for (i=rstart; i<rend; i++) { 47 PetscCall(MatGetRow(A,i,&nz,&cwork,&vwork)); 48 if (nz) { 49 j = 0; 50 while (cwork[j] < i) { 51 j++; nz--; 52 } 53 PetscCall(MatSetValues(M,1,&i,nz,cwork+j,vwork+j,INSERT_VALUES)); 54 } 55 PetscCall(MatRestoreRow(A,i,&nz,&cwork,&vwork)); 56 } 57 PetscCall(MatAssemblyBegin(M,MAT_FINAL_ASSEMBLY)); 58 PetscCall(MatAssemblyEnd(M,MAT_FINAL_ASSEMBLY)); 59 } else { 60 PetscCall(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 PetscCall(MatConvert_Basic(A,newtype,MAT_REUSE_MATRIX,&M)); 65 } 66 67 if (reuse == MAT_INPLACE_MATRIX) { 68 PetscCall(MatHeaderReplace(A,&M)); 69 } else *newmat = M; 70 PetscFunctionReturn(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; 85 86 PetscFunctionBegin; 87 if (reuse != MAT_REUSE_MATRIX) { 88 PetscCall(MatGetSize(A,&m,&n)); 89 PetscCall(MatGetLocalSize(A,&lm,&ln)); 90 PetscCall(PetscMalloc2(lm/bs,&d_nnz,lm/bs,&o_nnz)); 91 92 PetscCall(MatMarkDiagonal_SeqBAIJ(mpimat->A)); 93 for (i=0; i<lm/bs; i++) { 94 d_nnz[i] = Aa->i[i+1] - Aa->diag[i]; 95 o_nnz[i] = Ba->i[i+1] - Ba->i[i]; 96 } 97 98 PetscCall(MatCreate(PetscObjectComm((PetscObject)A),&M)); 99 PetscCall(MatSetSizes(M,lm,ln,m,n)); 100 PetscCall(MatSetType(M,MATMPISBAIJ)); 101 PetscCall(MatSeqSBAIJSetPreallocation(M,bs,0,d_nnz)); 102 PetscCall(MatMPISBAIJSetPreallocation(M,bs,0,d_nnz,0,o_nnz)); 103 104 PetscCall(PetscFree2(d_nnz,o_nnz)); 105 } else M = *newmat; 106 107 PetscCall(MatGetOwnershipRange(A,&rstart,&rend)); 108 PetscCall(MatSetOption(M,MAT_IGNORE_LOWER_TRIANGULAR,PETSC_TRUE)); 109 for (i=rstart; i<rend; i++) { 110 PetscCall(MatGetRow(A,i,&nz,&cwork,&vwork)); 111 PetscCall(MatSetValues(M,1,&i,nz,cwork,vwork,INSERT_VALUES)); 112 PetscCall(MatRestoreRow(A,i,&nz,&cwork,&vwork)); 113 } 114 PetscCall(MatAssemblyBegin(M,MAT_FINAL_ASSEMBLY)); 115 PetscCall(MatAssemblyEnd(M,MAT_FINAL_ASSEMBLY)); 116 117 if (reuse == MAT_INPLACE_MATRIX) { 118 PetscCall(MatHeaderReplace(A,&M)); 119 } else *newmat = M; 120 PetscFunctionReturn(0); 121 } 122