2: #include <petsc-private/matimpl.h>

  6: /*
  7:   MatConvert_Basic - Converts from any input format to another format. For
  8:   parallel formats, the new matrix distribution is determined by PETSc.

 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;

 22:   MatGetSize(mat,&m,&n);
 23:   MatGetLocalSize(mat,&lm,&ln);

 25:   if (ln == n) ln = PETSC_DECIDE; /* try to preserve column ownership */

 27:   MatCreate(PetscObjectComm((PetscObject)mat),&M);
 28:   MatSetSizes(M,lm,ln,m,n);
 29:   MatSetBlockSizesFromMats(M,mat,mat);
 30:   MatSetType(M,newtype);
 31:   MatGetOwnershipRange(mat,&rstart,&rend);

 33:   PetscObjectTypeCompare((PetscObject)M,MATSEQSBAIJ,&isseqsbaij);
 34:   PetscObjectTypeCompare((PetscObject)M,MATMPISBAIJ,&ismpisbaij);
 35:   if (isseqsbaij || ismpisbaij) {MatSetOption(M,MAT_IGNORE_LOWER_TRIANGULAR,PETSC_TRUE);}
 36:   PetscObjectTypeCompare((PetscObject)M,MATSEQBAIJ,&isseqbaij);
 37:   PetscObjectTypeCompare((PetscObject)M,MATMPIBAIJ,&ismpibaij);
 38:   PetscObjectTypeCompare((PetscObject)M,MATSEQDENSE,&isseqdense);
 39:   PetscObjectTypeCompare((PetscObject)M,MATMPIDENSE,&ismpidense);

 41:   if (isseqdense) {
 42:     MatSeqDenseSetPreallocation(M,NULL);
 43:   } else if (ismpidense) {
 44:     MatMPIDenseSetPreallocation(M,NULL);
 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;

 50:     PetscMalloc2(lm/prbs,&dnz,lm/prbs,&onz);
 51:     MatGetOwnershipRangeColumn(mat,&cstart,&cend);
 52:     for (i=0; i<lm; i+=prbs) {
 53:       MatGetRow(mat,rstart+i,&nz,&cwork,NULL);
 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:       MatRestoreRow(mat,rstart+i,&nz,&cwork,NULL);
 62:     }
 63:     MatXAIJSetPreallocation(M,PETSC_DECIDE,dnz,onz,dnz,onz);
 64:     PetscFree2(dnz,onz);
 65:   }

 67:   for (i=rstart; i<rend; i++) {
 68:     MatGetRow(mat,i,&nz,&cwork,&vwork);
 69:     MatSetValues(M,1,&i,nz,cwork,vwork,INSERT_VALUES);
 70:     MatRestoreRow(mat,i,&nz,&cwork,&vwork);
 71:   }
 72:   MatAssemblyBegin(M,MAT_FINAL_ASSEMBLY);
 73:   MatAssemblyEnd(M,MAT_FINAL_ASSEMBLY);

 75:   if (reuse == MAT_REUSE_MATRIX) {
 76:     MatHeaderReplace(mat,M);
 77:   } else {
 78:     *newmat = M;
 79:   }
 80:   return(0);
 81: }