2: #include <petsc-private/matimpl.h>        /*I "petscmat.h" I*/
  3: #include <petsc-private/vecimpl.h>  
  4: 
  7: PetscErrorCode MatConvert_Shell(Mat oldmat, const MatType newtype,MatReuse reuse,Mat *newmat)
  8: {
  9:   Mat            mat;
 10:   Vec            in,out;
 12:   PetscInt       i,M,m,*rows,start,end;
 13:   MPI_Comm       comm;
 14:   PetscScalar    *array,zero = 0.0,one = 1.0;

 17:   comm = ((PetscObject)oldmat)->comm;

 19:   MatGetOwnershipRange(oldmat,&start,&end);
 20:   VecCreateMPI(comm,end-start,PETSC_DECIDE,&in);
 21:   VecDuplicate(in,&out);
 22:   VecGetSize(in,&M);
 23:   VecGetLocalSize(in,&m);
 24:   PetscMalloc((m+1)*sizeof(PetscInt),&rows);
 25:   for (i=0; i<m; i++) {rows[i] = start + i;}

 27:   MatCreate(comm,&mat);
 28:   MatSetSizes(mat,m,M,M,M);
 29:   MatSetType(mat,newtype);
 30:   MatSetBlockSize(mat,oldmat->rmap->bs);
 31:   MatSetUp(mat);

 33:   for (i=0; i<M; i++) {
 34:     VecSet(in,zero);
 35:     VecSetValues(in,1,&i,&one,INSERT_VALUES);
 36:     VecAssemblyBegin(in);
 37:     VecAssemblyEnd(in);

 39:     MatMult(oldmat,in,out);
 40: 
 41:     VecGetArray(out,&array);
 42:     MatSetValues(mat,m,rows,1,&i,array,INSERT_VALUES);
 43:     VecRestoreArray(out,&array);

 45:   }
 46:   PetscFree(rows);
 47:   VecDestroy(&in);
 48:   VecDestroy(&out);
 49:   MatAssemblyBegin(mat,MAT_FINAL_ASSEMBLY);
 50:   MatAssemblyEnd(mat,MAT_FINAL_ASSEMBLY);

 52:   if (reuse == MAT_REUSE_MATRIX) {
 53:     MatHeaderReplace(oldmat,mat);
 54:   } else {
 55:     *newmat = mat;
 56:   }
 57:   return(0);
 58: }