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

  4: PETSC_INTERN PetscErrorCode MatConvert_Shell(Mat oldmat, MatType newtype,MatReuse reuse,Mat *newmat)
  5: {
  6:   Mat            mat;
  7:   Vec            in,out;
  9:   PetscInt       i,M,m,*rows,start,end;
 10:   MPI_Comm       comm;
 11:   PetscScalar    *array,zero = 0.0,one = 1.0;

 14:   PetscObjectGetComm((PetscObject)oldmat,&comm);

 16:   MatGetOwnershipRange(oldmat,&start,&end);
 17:   VecCreateMPI(comm,end-start,PETSC_DECIDE,&in);
 18:   VecDuplicate(in,&out);
 19:   VecGetSize(in,&M);
 20:   VecGetLocalSize(in,&m);
 21:   PetscMalloc1(m+1,&rows);
 22:   for (i=0; i<m; i++) rows[i] = start + i;

 24:   MatCreate(comm,&mat);
 25:   MatSetSizes(mat,m,M,M,M);
 26:   MatSetType(mat,newtype);
 27:   MatSetBlockSizesFromMats(mat,oldmat,oldmat);
 28:   MatSetUp(mat);

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

 36:     MatMult(oldmat,in,out);

 38:     VecGetArray(out,&array);
 39:     MatSetValues(mat,m,rows,1,&i,array,INSERT_VALUES);
 40:     VecRestoreArray(out,&array);

 42:   }
 43:   PetscFree(rows);
 44:   VecDestroy(&in);
 45:   VecDestroy(&out);
 46:   MatAssemblyBegin(mat,MAT_FINAL_ASSEMBLY);
 47:   MatAssemblyEnd(mat,MAT_FINAL_ASSEMBLY);

 49:   if (reuse == MAT_INPLACE_MATRIX) {
 50:     MatHeaderReplace(oldmat,&mat);
 51:   } else {
 52:     *newmat = mat;
 53:   }
 54:   return(0);
 55: }