1 2 #include "src/mat/matimpl.h" /*I "petscmat.h" I*/ 3 #include "vecimpl.h" 4 5 #undef __FUNCT__ 6 #define __FUNCT__ "MatConvert_Shell" 7 int MatConvert_Shell(Mat oldmat,const MatType newtype,Mat *newmat) { 8 Mat mat; 9 Vec in,out; 10 int ierr,i,M,m,*rows,start,end; 11 MPI_Comm comm; 12 PetscScalar *array,zero = 0.0,one = 1.0; 13 14 PetscFunctionBegin; 15 comm = oldmat->comm; 16 17 ierr = MatGetOwnershipRange(oldmat,&start,&end);CHKERRQ(ierr); 18 ierr = VecCreateMPI(comm,end-start,PETSC_DECIDE,&in);CHKERRQ(ierr); 19 ierr = VecDuplicate(in,&out);CHKERRQ(ierr); 20 ierr = VecGetSize(in,&M);CHKERRQ(ierr); 21 ierr = VecGetLocalSize(in,&m);CHKERRQ(ierr); 22 ierr = PetscMalloc((m+1)*sizeof(int),&rows);CHKERRQ(ierr); 23 for (i=0; i<m; i++) {rows[i] = start + i;} 24 25 ierr = MatCreate(comm,m,M,M,M,&mat);CHKERRQ(ierr); 26 ierr = MatSetType(mat,newtype);CHKERRQ(ierr); 27 28 for (i=0; i<M; i++) { 29 ierr = VecSet(&zero,in);CHKERRQ(ierr); 30 ierr = VecSetValues(in,1,&i,&one,INSERT_VALUES);CHKERRQ(ierr); 31 ierr = VecAssemblyBegin(in);CHKERRQ(ierr); 32 ierr = VecAssemblyEnd(in);CHKERRQ(ierr); 33 34 ierr = MatMult(oldmat,in,out);CHKERRQ(ierr); 35 36 ierr = VecGetArray(out,&array);CHKERRQ(ierr); 37 ierr = MatSetValues(mat,m,rows,1,&i,array,INSERT_VALUES);CHKERRQ(ierr); 38 ierr = VecRestoreArray(out,&array);CHKERRQ(ierr); 39 40 } 41 ierr = PetscFree(rows);CHKERRQ(ierr); 42 ierr = VecDestroy(in);CHKERRQ(ierr); 43 ierr = VecDestroy(out);CHKERRQ(ierr); 44 ierr = MatAssemblyBegin(mat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 45 ierr = MatAssemblyEnd(mat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 46 47 /* Fake support for "inplace" convert. */ 48 if (*newmat == oldmat) { 49 ierr = MatDestroy(oldmat);CHKERRQ(ierr); 50 } 51 *newmat = mat; 52 PetscFunctionReturn(0); 53 } 54 55 56 57