Actual source code: ex88.c
petsc-3.10.5 2019-03-28
2: static char help[] = "Tests MatShift(), MatScale(), and MatDiagonalScale() for SHELL and NEST matrices\n\n";
4: #include <petscmat.h>
6: typedef struct _n_User *User;
7: struct _n_User {
8: Mat B;
9: };
11: static PetscErrorCode MatView_User(Mat A,PetscViewer viewer)
12: {
13: User user;
17: MatShellGetContext(A,&user);
18: MatView(user->B,viewer);
19: return(0);
20: }
22: static PetscErrorCode MatMult_User(Mat A,Vec X,Vec Y)
23: {
24: User user;
28: MatShellGetContext(A,&user);
29: MatMult(user->B,X,Y);
30: return(0);
31: }
33: static PetscErrorCode MatMultTranspose_User(Mat A,Vec X,Vec Y)
34: {
35: User user;
39: MatShellGetContext(A,&user);
40: MatMultTranspose(user->B,X,Y);
41: return(0);
42: }
44: static PetscErrorCode MatGetDiagonal_User(Mat A,Vec X)
45: {
46: User user;
50: MatShellGetContext(A,&user);
51: MatGetDiagonal(user->B,X);
52: return(0);
53: }
55: static PetscErrorCode TestMatrix(Mat A,Vec X,Vec Y,Vec Z)
56: {
58: Vec W1,W2,diff;
59: Mat E;
60: const char *mattypename;
61: PetscViewer viewer = PETSC_VIEWER_STDOUT_WORLD;
62: PetscScalar diag[2] = { 2.9678190300000000e+08, 1.4173141580000000e+09};
63: PetscScalar multadd[2] = {-6.8966198500000000e+08,-2.0310609940000000e+09};
64: PetscScalar multtadd[2] = {-9.1052873900000000e+08,-1.8101942400000000e+09};
65: PetscReal nrm,norm1,norminf;
68: PetscObjectGetType((PetscObject)A,&mattypename);
69: PetscViewerASCIIPrintf(viewer,"\nMatrix of type: %s\n",mattypename);
70: VecDuplicate(X,&W1);
71: VecDuplicate(X,&W2);
72: MatScale(A,31);
73: MatShift(A,37);
74: MatDiagonalScale(A,X,Y);
75: MatScale(A,41);
76: MatDiagonalScale(A,Y,Z);
77: MatComputeExplicitOperator(A,&E);
79: MatView(E,viewer);
80: PetscViewerASCIIPrintf(viewer,"Testing MatMult + MatMultTranspose\n");
81: MatMult(A,Z,W1);
82: MatMultTranspose(A,W1,W2);
83: VecView(W2,viewer);
85: PetscViewerASCIIPrintf(viewer,"Testing MatMultAdd\n");
86: VecCreateSeqWithArray(PETSC_COMM_SELF,1,2,multadd,&diff);
87: VecSet(W1,-1.0);
88: MatMultAdd(A,W1,W1,W2);
89: VecView(W2,viewer);
90: VecAXPY(W2,-1.0,diff);
91: VecNorm(W2,NORM_2,&nrm);
92: #if defined(PETSC_USE_REAL_DOUBLE) || defined(PETSC_USE_REAL___FLOAT128)
93: if (nrm > PETSC_SMALL) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"MatMultAdd(A,x,x,y) produces incorrect result");
94: #endif
96: VecSet(W2,-1.0);
97: MatMultAdd(A,W1,W2,W2);
98: VecView(W2,viewer);
99: VecAXPY(W2,-1.0,diff);
100: VecNorm(W2,NORM_2,&nrm);
101: #if defined(PETSC_USE_REAL_DOUBLE) || defined(PETSC_USE_REAL___FLOAT128)
102: if (nrm > PETSC_SMALL) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"MatMultAdd(A,x,y,y) produces incorrect result");
103: #endif
104: VecDestroy(&diff);
106: PetscViewerASCIIPrintf(viewer,"Testing MatMultTranposeAdd\n");
107: VecCreateSeqWithArray(PETSC_COMM_SELF,1,2,multtadd,&diff);
109: VecSet(W1,-1.0);
110: MatMultTransposeAdd(A,W1,W1,W2);
111: VecView(W2,viewer);
112: VecAXPY(W2,-1.0,diff);
113: VecNorm(W2,NORM_2,&nrm);
114: #if defined(PETSC_USE_REAL_DOUBLE) || defined(PETSC_USE_REAL___FLOAT128)
115: if (nrm > PETSC_SMALL) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"MatMultTranposeAdd(A,x,x,y) produces incorrect result");
116: #endif
118: VecSet(W2,-1.0);
119: MatMultTransposeAdd(A,W1,W2,W2);
120: VecView(W2,viewer);
121: VecAXPY(W2,-1.0,diff);
122: VecNorm(W2,NORM_2,&nrm);
123: #if defined(PETSC_USE_REAL_DOUBLE) || defined(PETSC_USE_REAL___FLOAT128)
124: if (nrm > PETSC_SMALL) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"MatMultTranposeAdd(A,x,y,y) produces incorrect result");
125: #endif
126: VecDestroy(&diff);
128: PetscViewerASCIIPrintf(viewer,"Testing MatGetDiagonal\n");
129: MatGetDiagonal(A,W2);
130: VecView(W2,viewer);
131: VecCreateSeqWithArray(PETSC_COMM_SELF,1,2,diag,&diff);
132: VecAXPY(diff,-1.0,W2);
133: VecNorm(diff,NORM_2,&nrm);
134: if (nrm > PETSC_SMALL) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"MatGetDiagonal() produces incorrect result");
135: VecDestroy(&diff);
136: MatNorm(A,NORM_1,&norm1);
137: MatNorm(A,NORM_INFINITY,&norminf);
138: PetscViewerASCIIPrintf(viewer,"Norm_1: %g, Norm_infty %g\n",(double)norm1,(double)norminf);
139: /* MATSHELL does not support MatDiagonalSet after MatScale */
140: if (strncmp(mattypename, "shell", 5)) {
141: MatDiagonalSet(A,X,INSERT_VALUES);
142: MatGetDiagonal(A,W1);
143: VecView(W1,viewer);
144: } else {
145: PetscViewerASCIIPrintf(viewer,"MatDiagonalSet not tested on MATSHELL\n");
146: }
148: MatDestroy(&E);
149: VecDestroy(&W1);
150: VecDestroy(&W2);
151: return(0);
152: }
154: int main(int argc,char **args)
155: {
156: const PetscScalar xvals[] = {11,13},yvals[] = {17,19},zvals[] = {23,29};
157: const PetscInt inds[] = {0,1};
158: PetscScalar avals[] = {2,3,5,7};
159: Mat A,S,D[4],N;
160: Vec X,Y,Z;
161: User user;
162: PetscInt i;
163: PetscErrorCode ierr;
165: PetscInitialize(&argc,&args,(char*)0,help);if (ierr) return ierr;
166: MatCreateSeqAIJ(PETSC_COMM_WORLD,2,2,2,NULL,&A);
167: MatSetUp(A);
168: MatSetValues(A,2,inds,2,inds,avals,INSERT_VALUES);
169: VecCreateSeq(PETSC_COMM_WORLD,2,&X);
170: VecDuplicate(X,&Y);
171: VecDuplicate(X,&Z);
172: VecSetValues(X,2,inds,xvals,INSERT_VALUES);
173: VecSetValues(Y,2,inds,yvals,INSERT_VALUES);
174: VecSetValues(Z,2,inds,zvals,INSERT_VALUES);
175: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
176: MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
177: VecAssemblyBegin(X);
178: VecAssemblyBegin(Y);
179: VecAssemblyBegin(Z);
180: VecAssemblyEnd(X);
181: VecAssemblyEnd(Y);
182: VecAssemblyEnd(Z);
184: PetscNew(&user);
185: user->B = A;
187: MatCreateShell(PETSC_COMM_WORLD,2,2,2,2,user,&S);
188: MatSetUp(S);
189: MatShellSetOperation(S,MATOP_VIEW,(void (*)(void))MatView_User);
190: MatShellSetOperation(S,MATOP_MULT,(void (*)(void))MatMult_User);
191: MatShellSetOperation(S,MATOP_MULT_TRANSPOSE,(void (*)(void))MatMultTranspose_User);
192: MatShellSetOperation(S,MATOP_GET_DIAGONAL,(void (*)(void))MatGetDiagonal_User);
194: for (i=0; i<4; i++) {
195: MatCreateSeqDense(PETSC_COMM_WORLD,1,1,&avals[i],&D[i]);
196: }
197: MatCreateNest(PETSC_COMM_WORLD,2,NULL,2,NULL,D,&N);
198: MatSetUp(N);
200: TestMatrix(S,X,Y,Z);
201: TestMatrix(A,X,Y,Z);
202: TestMatrix(N,X,Y,Z);
204: for (i=0; i<4; i++) {MatDestroy(&D[i]);}
205: MatDestroy(&A);
206: MatDestroy(&S);
207: MatDestroy(&N);
208: VecDestroy(&X);
209: VecDestroy(&Y);
210: VecDestroy(&Z);
211: PetscFree(user);
212: PetscFinalize();
213: return ierr;
214: }
217: /*TEST
219: test:
221: TEST*/