Actual source code: shellcnv.c
petsc-3.13.6 2020-09-29
1: #include <petsc/private/matimpl.h>
3: PetscErrorCode MatConvert_Shell(Mat oldmat, MatType newtype,MatReuse reuse,Mat *newmat)
4: {
5: Mat mat;
6: Vec in,out;
7: MPI_Comm comm;
8: PetscScalar *array;
9: PetscInt *dnnz,*onnz,*dnnzu,*onnzu;
10: PetscInt cst,Nbs,mbs,nbs,rbs,cbs;
11: PetscInt im,i,m,n,M,N,*rows,start;
15: PetscObjectGetComm((PetscObject)oldmat,&comm);
17: MatGetOwnershipRange(oldmat,&start,NULL);
18: MatGetOwnershipRangeColumn(oldmat,&cst,NULL);
19: MatCreateVecs(oldmat,&in,&out);
20: MatGetLocalSize(oldmat,&m,&n);
21: MatGetSize(oldmat,&M,&N);
22: PetscMalloc1(m,&rows);
24: MatCreate(comm,&mat);
25: MatSetSizes(mat,m,n,M,N);
26: MatSetType(mat,newtype);
27: MatSetBlockSizesFromMats(mat,oldmat,oldmat);
28: MatGetBlockSizes(mat,&rbs,&cbs);
29: mbs = m/rbs;
30: nbs = n/cbs;
31: Nbs = N/cbs;
32: cst = cst/cbs;
33: PetscMalloc4(mbs,&dnnz,mbs,&onnz,mbs,&dnnzu,mbs,&onnzu);
34: for (i=0; i<mbs; i++) {
35: dnnz[i] = nbs;
36: onnz[i] = Nbs - nbs;
37: dnnzu[i] = PetscMax(nbs - i,0);
38: onnzu[i] = PetscMax(Nbs - (cst + nbs),0);
39: }
40: MatXAIJSetPreallocation(mat,PETSC_DECIDE,dnnz,onnz,dnnzu,onnzu);
41: PetscFree4(dnnz,onnz,dnnzu,onnzu);
42: VecSetOption(in,VEC_IGNORE_OFF_PROC_ENTRIES,PETSC_TRUE);
43: MatSetUp(mat);
44: for (i=0; i<N; i++) {
45: PetscInt j;
47: VecZeroEntries(in);
48: VecSetValue(in,i,1.,INSERT_VALUES);
49: VecAssemblyBegin(in);
50: VecAssemblyEnd(in);
51: MatMult(oldmat,in,out);
52: VecGetArray(out,&array);
53: for (j=0, im = 0; j<m; j++) {
54: if (PetscAbsScalar(array[j]) == 0.0) continue;
55: rows[im] = j+start;
56: array[im] = array[j];
57: im++;
58: }
59: MatSetValues(mat,im,rows,1,&i,array,INSERT_VALUES);
60: VecRestoreArray(out,&array);
61: }
62: PetscFree(rows);
63: VecDestroy(&in);
64: VecDestroy(&out);
65: MatAssemblyBegin(mat,MAT_FINAL_ASSEMBLY);
66: MatAssemblyEnd(mat,MAT_FINAL_ASSEMBLY);
67: if (reuse == MAT_INPLACE_MATRIX) {
68: MatHeaderReplace(oldmat,&mat);
69: } else {
70: *newmat = mat;
71: }
72: return(0);
73: }
75: static PetscErrorCode MatGetDiagonal_CF(Mat A,Vec X)
76: {
77: Mat B;
81: MatShellGetContext(A,&B);
82: if (!B) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_PLIB,"Missing user matrix");
83: MatGetDiagonal(B,X);
84: return(0);
85: }
87: static PetscErrorCode MatMult_CF(Mat A,Vec X,Vec Y)
88: {
89: Mat B;
93: MatShellGetContext(A,&B);
94: if (!B) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_PLIB,"Missing user matrix");
95: MatMult(B,X,Y);
96: return(0);
97: }
99: static PetscErrorCode MatMultTranspose_CF(Mat A,Vec X,Vec Y)
100: {
101: Mat B;
105: MatShellGetContext(A,&B);
106: if (!B) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_PLIB,"Missing user matrix");
107: MatMultTranspose(B,X,Y);
108: return(0);
109: }
111: static PetscErrorCode MatDestroy_CF(Mat A)
112: {
113: Mat B;
117: MatShellGetContext(A,&B);
118: if (!B) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_PLIB,"Missing user matrix");
119: MatDestroy(&B);
120: return(0);
121: }
123: PetscErrorCode MatConvertFrom_Shell(Mat A, MatType newtype,MatReuse reuse,Mat *B)
124: {
125: Mat M;
126: PetscBool flg;
130: PetscStrcmp(newtype,MATSHELL,&flg);
131: if (!flg) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Only conversion to MATSHELL");
132: if (reuse == MAT_INITIAL_MATRIX) {
133: PetscObjectReference((PetscObject)A);
134: MatCreateShell(PetscObjectComm((PetscObject)A),A->rmap->n,A->cmap->n,A->rmap->N,A->cmap->N,A,&M);
135: MatSetBlockSizesFromMats(M,A,A);
136: MatShellSetOperation(M,MATOP_MULT, (void (*)(void))MatMult_CF);
137: MatShellSetOperation(M,MATOP_MULT_TRANSPOSE,(void (*)(void))MatMultTranspose_CF);
138: MatShellSetOperation(M,MATOP_GET_DIAGONAL, (void (*)(void))MatGetDiagonal_CF);
139: MatShellSetOperation(M,MATOP_DESTROY, (void (*)(void))MatDestroy_CF);
140: *B = M;
141: } else SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_SUP,"Not implemented");
142: return(0);
143: }