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