Actual source code: shellcnv.c

      1: #define PETSCMAT_DLL

     3:  #include src/mat/matimpl.h
     4:  #include vecimpl.h
      5:
      8: PetscErrorCode MatConvert_Shell(Mat oldmat, MatType newtype,MatReuse reuse,Mat *newmat)
      9: {
     10:   Mat            mat;
     11:   Vec            in,out;
     13:   PetscInt       i,M,m,*rows,start,end;
     14:   MPI_Comm       comm;
     15:   PetscScalar    *array,zero = 0.0,one = 1.0;

     18:   comm = oldmat->comm;

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

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

     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: }