1 2 #include <petsc/private/matimpl.h> /*I "petscmat.h" I*/ 3 4 PetscErrorCode MatConvert_Shell(Mat oldmat, MatType newtype,MatReuse reuse,Mat *newmat) 5 { 6 Mat mat; 7 Vec in,out; 8 PetscErrorCode ierr; 9 PetscInt i,m,n,M,N,*rows,start; 10 MPI_Comm comm; 11 PetscScalar *array; 12 13 PetscFunctionBegin; 14 ierr = PetscObjectGetComm((PetscObject)oldmat,&comm);CHKERRQ(ierr); 15 16 ierr = MatGetOwnershipRange(oldmat,&start,NULL);CHKERRQ(ierr); 17 ierr = MatCreateVecs(oldmat,&in,&out);CHKERRQ(ierr); 18 ierr = MatGetLocalSize(oldmat,&m,&n);CHKERRQ(ierr); 19 ierr = MatGetSize(oldmat,&M,&N);CHKERRQ(ierr); 20 ierr = PetscMalloc1(m,&rows);CHKERRQ(ierr); 21 for (i=0; i<m; i++) rows[i] = start + i; 22 23 ierr = MatCreate(comm,&mat);CHKERRQ(ierr); 24 ierr = MatSetSizes(mat,m,n,M,N);CHKERRQ(ierr); 25 ierr = MatSetType(mat,newtype);CHKERRQ(ierr); 26 ierr = MatSetBlockSizesFromMats(mat,oldmat,oldmat);CHKERRQ(ierr); 27 ierr = MatSetUp(mat);CHKERRQ(ierr); 28 ierr = VecSetOption(in,VEC_IGNORE_OFF_PROC_ENTRIES,PETSC_TRUE);CHKERRQ(ierr); 29 30 for (i=0; i<N; i++) { 31 ierr = VecZeroEntries(in);CHKERRQ(ierr); 32 ierr = VecSetValue(in,i,1.,INSERT_VALUES);CHKERRQ(ierr); 33 ierr = VecAssemblyBegin(in);CHKERRQ(ierr); 34 ierr = VecAssemblyEnd(in);CHKERRQ(ierr); 35 36 ierr = MatMult(oldmat,in,out);CHKERRQ(ierr); 37 38 ierr = VecGetArray(out,&array);CHKERRQ(ierr); 39 ierr = MatSetValues(mat,m,rows,1,&i,array,INSERT_VALUES);CHKERRQ(ierr); 40 ierr = VecRestoreArray(out,&array);CHKERRQ(ierr); 41 42 } 43 ierr = PetscFree(rows);CHKERRQ(ierr); 44 ierr = VecDestroy(&in);CHKERRQ(ierr); 45 ierr = VecDestroy(&out);CHKERRQ(ierr); 46 ierr = MatAssemblyBegin(mat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 47 ierr = MatAssemblyEnd(mat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 48 49 if (reuse == MAT_INPLACE_MATRIX) { 50 ierr = MatHeaderReplace(oldmat,&mat);CHKERRQ(ierr); 51 } else { 52 *newmat = mat; 53 } 54 PetscFunctionReturn(0); 55 } 56