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