1 2 #include <petsc-private/matimpl.h> 3 4 #undef __FUNCT__ 5 #define __FUNCT__ "MatConvert_Basic" 6 /* 7 MatConvert_Basic - Converts from any input format to another format. For 8 parallel formats, the new matrix distribution is determined by PETSc. 9 10 Does not do preallocation so in general will be slow 11 */ 12 PetscErrorCode MatConvert_Basic(Mat mat, MatType newtype,MatReuse reuse,Mat *newmat) 13 { 14 Mat M; 15 const PetscScalar *vwork; 16 PetscErrorCode ierr; 17 PetscInt i,j,nz,m,n,rstart,rend,lm,ln,prbs,pcbs,cstart,cend,*dnz,*onz; 18 const PetscInt *cwork; 19 PetscBool isseqsbaij,ismpisbaij,isseqbaij,ismpibaij,isseqdense,ismpidense; 20 21 PetscFunctionBegin; 22 ierr = MatGetSize(mat,&m,&n);CHKERRQ(ierr); 23 ierr = MatGetLocalSize(mat,&lm,&ln);CHKERRQ(ierr); 24 25 if (ln == n) ln = PETSC_DECIDE; /* try to preserve column ownership */ 26 27 ierr = MatCreate(PetscObjectComm((PetscObject)mat),&M);CHKERRQ(ierr); 28 ierr = MatSetSizes(M,lm,ln,m,n);CHKERRQ(ierr); 29 ierr = MatSetBlockSizesFromMats(M,mat,mat);CHKERRQ(ierr); 30 ierr = MatSetType(M,newtype);CHKERRQ(ierr); 31 ierr = MatGetOwnershipRange(mat,&rstart,&rend);CHKERRQ(ierr); 32 33 ierr = PetscObjectTypeCompare((PetscObject)M,MATSEQSBAIJ,&isseqsbaij);CHKERRQ(ierr); 34 ierr = PetscObjectTypeCompare((PetscObject)M,MATMPISBAIJ,&ismpisbaij);CHKERRQ(ierr); 35 if (isseqsbaij || ismpisbaij) {ierr = MatSetOption(M,MAT_IGNORE_LOWER_TRIANGULAR,PETSC_TRUE);CHKERRQ(ierr);} 36 ierr = PetscObjectTypeCompare((PetscObject)M,MATSEQBAIJ,&isseqbaij);CHKERRQ(ierr); 37 ierr = PetscObjectTypeCompare((PetscObject)M,MATMPIBAIJ,&ismpibaij);CHKERRQ(ierr); 38 ierr = PetscObjectTypeCompare((PetscObject)M,MATSEQDENSE,&isseqdense);CHKERRQ(ierr); 39 ierr = PetscObjectTypeCompare((PetscObject)M,MATMPIDENSE,&ismpidense);CHKERRQ(ierr); 40 41 if (isseqdense) { 42 ierr = MatSeqDenseSetPreallocation(M,NULL);CHKERRQ(ierr); 43 } else if (ismpidense) { 44 ierr = MatMPIDenseSetPreallocation(M,NULL);CHKERRQ(ierr); 45 } else { 46 /* Preallocation block sizes. (S)BAIJ matrices will have one index per block. */ 47 prbs = (isseqbaij || ismpibaij || isseqsbaij || ismpisbaij) ? PetscAbs(M->rmap->bs) : 1; 48 pcbs = (isseqbaij || ismpibaij || isseqsbaij || ismpisbaij) ? PetscAbs(M->cmap->bs) : 1; 49 50 ierr = PetscMalloc2(lm/prbs,&dnz,lm/prbs,&onz);CHKERRQ(ierr); 51 ierr = MatGetOwnershipRangeColumn(mat,&cstart,&cend);CHKERRQ(ierr); 52 for (i=0; i<lm; i+=prbs) { 53 ierr = MatGetRow(mat,rstart+i,&nz,&cwork,NULL);CHKERRQ(ierr); 54 dnz[i] = 0; 55 onz[i] = 0; 56 for (j=0; j<nz; j+=pcbs) { 57 if ((isseqsbaij || ismpisbaij) && cwork[j] < rstart+i) continue; 58 if (cstart <= cwork[j] && cwork[j] < cend) dnz[i]++; 59 else onz[i]++; 60 } 61 ierr = MatRestoreRow(mat,rstart+i,&nz,&cwork,NULL);CHKERRQ(ierr); 62 } 63 ierr = MatXAIJSetPreallocation(M,PETSC_DECIDE,dnz,onz,dnz,onz);CHKERRQ(ierr); 64 ierr = PetscFree2(dnz,onz);CHKERRQ(ierr); 65 } 66 67 for (i=rstart; i<rend; i++) { 68 ierr = MatGetRow(mat,i,&nz,&cwork,&vwork);CHKERRQ(ierr); 69 ierr = MatSetValues(M,1,&i,nz,cwork,vwork,INSERT_VALUES);CHKERRQ(ierr); 70 ierr = MatRestoreRow(mat,i,&nz,&cwork,&vwork);CHKERRQ(ierr); 71 } 72 ierr = MatAssemblyBegin(M,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 73 ierr = MatAssemblyEnd(M,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 74 75 if (reuse == MAT_REUSE_MATRIX) { 76 ierr = MatHeaderReplace(mat,M);CHKERRQ(ierr); 77 } else { 78 *newmat = M; 79 } 80 PetscFunctionReturn(0); 81 } 82