1 #include <petsc/private/matimpl.h> /*I "petscmat.h" I*/ 2 3 PetscErrorCode MatConvert_Shell(Mat oldmat, MatType newtype,MatReuse reuse,Mat *newmat) 4 { 5 Mat mat; 6 Vec in,out; 7 PetscErrorCode ierr; 8 PetscInt i,m,n,M,N,*rows,start; 9 MPI_Comm comm; 10 PetscScalar *array; 11 12 PetscFunctionBegin; 13 ierr = PetscObjectGetComm((PetscObject)oldmat,&comm);CHKERRQ(ierr); 14 15 ierr = MatGetOwnershipRange(oldmat,&start,NULL);CHKERRQ(ierr); 16 ierr = MatCreateVecs(oldmat,&in,&out);CHKERRQ(ierr); 17 ierr = MatGetLocalSize(oldmat,&m,&n);CHKERRQ(ierr); 18 ierr = MatGetSize(oldmat,&M,&N);CHKERRQ(ierr); 19 ierr = PetscMalloc1(m,&rows);CHKERRQ(ierr); 20 for (i=0; i<m; i++) rows[i] = start + i; 21 22 ierr = MatCreate(comm,&mat);CHKERRQ(ierr); 23 ierr = MatSetSizes(mat,m,n,M,N);CHKERRQ(ierr); 24 ierr = MatSetType(mat,newtype);CHKERRQ(ierr); 25 ierr = MatSetBlockSizesFromMats(mat,oldmat,oldmat);CHKERRQ(ierr); 26 ierr = MatSetUp(mat);CHKERRQ(ierr); 27 ierr = VecSetOption(in,VEC_IGNORE_OFF_PROC_ENTRIES,PETSC_TRUE);CHKERRQ(ierr); 28 29 for (i=0; i<N; i++) { 30 ierr = VecZeroEntries(in);CHKERRQ(ierr); 31 ierr = VecSetValue(in,i,1.,INSERT_VALUES);CHKERRQ(ierr); 32 ierr = VecAssemblyBegin(in);CHKERRQ(ierr); 33 ierr = VecAssemblyEnd(in);CHKERRQ(ierr); 34 35 ierr = MatMult(oldmat,in,out);CHKERRQ(ierr); 36 37 ierr = VecGetArray(out,&array);CHKERRQ(ierr); 38 ierr = MatSetValues(mat,m,rows,1,&i,array,INSERT_VALUES);CHKERRQ(ierr); 39 ierr = VecRestoreArray(out,&array);CHKERRQ(ierr); 40 41 } 42 ierr = PetscFree(rows);CHKERRQ(ierr); 43 ierr = VecDestroy(&in);CHKERRQ(ierr); 44 ierr = VecDestroy(&out);CHKERRQ(ierr); 45 ierr = MatAssemblyBegin(mat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 46 ierr = MatAssemblyEnd(mat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 47 48 if (reuse == MAT_INPLACE_MATRIX) { 49 ierr = MatHeaderReplace(oldmat,&mat);CHKERRQ(ierr); 50 } else { 51 *newmat = mat; 52 } 53 PetscFunctionReturn(0); 54 } 55 56 static PetscErrorCode MatGetDiagonal_CF(Mat A,Vec X) 57 { 58 Mat B; 59 PetscErrorCode ierr; 60 61 PetscFunctionBegin; 62 ierr = MatShellGetContext(A,&B);CHKERRQ(ierr); 63 if (!B) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_PLIB,"Missing user matrix"); 64 ierr = MatGetDiagonal(B,X);CHKERRQ(ierr); 65 PetscFunctionReturn(0); 66 } 67 68 static PetscErrorCode MatMult_CF(Mat A,Vec X,Vec Y) 69 { 70 Mat B; 71 PetscErrorCode ierr; 72 73 PetscFunctionBegin; 74 ierr = MatShellGetContext(A,&B);CHKERRQ(ierr); 75 if (!B) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_PLIB,"Missing user matrix"); 76 ierr = MatMult(B,X,Y);CHKERRQ(ierr); 77 PetscFunctionReturn(0); 78 } 79 80 static PetscErrorCode MatMultTranspose_CF(Mat A,Vec X,Vec Y) 81 { 82 Mat B; 83 PetscErrorCode ierr; 84 85 PetscFunctionBegin; 86 ierr = MatShellGetContext(A,&B);CHKERRQ(ierr); 87 if (!B) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_PLIB,"Missing user matrix"); 88 ierr = MatMultTranspose(B,X,Y);CHKERRQ(ierr); 89 PetscFunctionReturn(0); 90 } 91 92 static PetscErrorCode MatDestroy_CF(Mat A) 93 { 94 Mat B; 95 PetscErrorCode ierr; 96 97 PetscFunctionBegin; 98 ierr = MatShellGetContext(A,&B);CHKERRQ(ierr); 99 if (!B) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_PLIB,"Missing user matrix"); 100 ierr = MatDestroy(&B);CHKERRQ(ierr); 101 PetscFunctionReturn(0); 102 } 103 104 PetscErrorCode MatConvertFrom_Shell(Mat A, MatType newtype,MatReuse reuse,Mat *B) 105 { 106 Mat M; 107 PetscBool flg; 108 PetscErrorCode ierr; 109 110 PetscFunctionBegin; 111 ierr = PetscStrcmp(newtype,MATSHELL,&flg);CHKERRQ(ierr); 112 if (!flg) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Only conversion to MATSHELL"); 113 if (reuse == MAT_INITIAL_MATRIX) { 114 ierr = PetscObjectReference((PetscObject)A);CHKERRQ(ierr); 115 ierr = MatCreateShell(PetscObjectComm((PetscObject)A),A->rmap->n,A->cmap->n,A->rmap->N,A->cmap->N,A,&M);CHKERRQ(ierr); 116 ierr = MatShellSetOperation(M,MATOP_MULT, (void (*)(void))MatMult_CF);CHKERRQ(ierr); 117 ierr = MatShellSetOperation(M,MATOP_MULT_TRANSPOSE,(void (*)(void))MatMultTranspose_CF);CHKERRQ(ierr); 118 ierr = MatShellSetOperation(M,MATOP_GET_DIAGONAL, (void (*)(void))MatGetDiagonal_CF);CHKERRQ(ierr); 119 ierr = MatShellSetOperation(M,MATOP_DESTROY, (void (*)(void))MatDestroy_CF);CHKERRQ(ierr); 120 *B = M; 121 } else SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_SUP,"Not implemented"); 122 PetscFunctionReturn(0); 123 } 124