Actual source code: ex89.c
petsc-3.12.5 2020-03-29
1: static char help[] ="Tests MatPtAP() for MPIMAIJ and MPIAIJ \n ";
3: #include <petscdmda.h>
5: int main(int argc,char **argv)
6: {
8: DM coarsedm,finedm;
9: PetscMPIInt size,rank;
10: PetscInt M,N,Z,i,nrows;
11: PetscScalar one = 1.0;
12: PetscReal fill=2.0;
13: Mat A,P,C;
14: PetscScalar *array,none = -1.0,alpha;
15: Vec x,v1,v2,v3,v4;
16: PetscReal norm,norm_tmp,norm_tmp1,tol=100.*PETSC_MACHINE_EPSILON;
17: PetscRandom rdm;
18: PetscBool Test_3D=PETSC_FALSE,flg;
19: const PetscInt *ia,*ja;
20: PetscInt dof;
21: MPI_Comm comm;
23: PetscInitialize(&argc,&argv,NULL,help);if (ierr) return ierr;
24: comm = PETSC_COMM_WORLD;
25: MPI_Comm_rank(comm,&rank);
26: MPI_Comm_size(comm,&size);
27: M = 10; N = 10; Z = 10;
28: dof = 10;
30: PetscOptionsGetBool(NULL,NULL,"-test_3D",&Test_3D,NULL);
31: PetscOptionsGetInt(NULL,NULL,"-M",&M,NULL);
32: PetscOptionsGetInt(NULL,NULL,"-N",&N,NULL);
33: PetscOptionsGetInt(NULL,NULL,"-Z",&Z,NULL);
34: /* Set up distributed array for fine grid */
35: if (!Test_3D) {
36: DMDACreate2d(comm,DM_BOUNDARY_NONE,DM_BOUNDARY_NONE,DMDA_STENCIL_STAR,M,N,PETSC_DECIDE,PETSC_DECIDE,dof,1,NULL,NULL,&coarsedm);
37: } else {
38: DMDACreate3d(comm,DM_BOUNDARY_NONE,DM_BOUNDARY_NONE,DM_BOUNDARY_NONE,DMDA_STENCIL_STAR,M,N,Z,PETSC_DECIDE,PETSC_DECIDE,PETSC_DECIDE,dof,1,NULL,NULL,NULL,&coarsedm);
39: }
40: DMSetFromOptions(coarsedm);
41: DMSetUp(coarsedm);
43: /* This makes sure the coarse DMDA has the same partition as the fine DMDA */
44: DMRefine(coarsedm,PetscObjectComm((PetscObject)coarsedm),&finedm);
46: /*------------------------------------------------------------*/
47: DMSetMatType(finedm,MATAIJ);
48: DMCreateMatrix(finedm,&A);
50: /* set val=one to A */
51: if (size == 1) {
52: MatGetRowIJ(A,0,PETSC_FALSE,PETSC_FALSE,&nrows,&ia,&ja,&flg);
53: if (flg) {
54: MatSeqAIJGetArray(A,&array);
55: for (i=0; i<ia[nrows]; i++) array[i] = one;
56: MatSeqAIJRestoreArray(A,&array);
57: }
58: MatRestoreRowIJ(A,0,PETSC_FALSE,PETSC_FALSE,&nrows,&ia,&ja,&flg);
59: } else {
60: Mat AA,AB;
61: MatMPIAIJGetSeqAIJ(A,&AA,&AB,NULL);
62: MatGetRowIJ(AA,0,PETSC_FALSE,PETSC_FALSE,&nrows,&ia,&ja,&flg);
63: if (flg) {
64: MatSeqAIJGetArray(AA,&array);
65: for (i=0; i<ia[nrows]; i++) array[i] = one;
66: MatSeqAIJRestoreArray(AA,&array);
67: }
68: MatRestoreRowIJ(AA,0,PETSC_FALSE,PETSC_FALSE,&nrows,&ia,&ja,&flg);
69: MatGetRowIJ(AB,0,PETSC_FALSE,PETSC_FALSE,&nrows,&ia,&ja,&flg);
70: if (flg) {
71: MatSeqAIJGetArray(AB,&array);
72: for (i=0; i<ia[nrows]; i++) array[i] = one;
73: MatSeqAIJRestoreArray(AB,&array);
74: }
75: MatRestoreRowIJ(AB,0,PETSC_FALSE,PETSC_FALSE,&nrows,&ia,&ja,&flg);
76: }
77: /* Create interpolation between the fine and coarse grids */
78: DMCreateInterpolation(coarsedm,finedm,&P,NULL);
79: /* Create vectors v1 and v2 that are compatible with A */
80: MatCreateVecs(A,&v1,NULL);
81: VecDuplicate(v1,&v2);
82: PetscRandomCreate(comm,&rdm);
83: PetscRandomSetFromOptions(rdm);
85: /* Test P^T * A * P - MatPtAP() */
86: /*------------------------------*/
87: MatPtAP(A,P,MAT_INITIAL_MATRIX,fill,&C);
88: /* Test MAT_REUSE_MATRIX - reuse symbolic C */
89: alpha=1.0;
90: for (i=0; i<1; i++) {
91: alpha -= 0.1;
92: MatScale(A,alpha);
93: MatPtAP(A,P,MAT_REUSE_MATRIX,fill,&C);
94: }
96: /* Free intermediate data structures created for reuse of C=Pt*A*P */
97: MatFreeIntermediateDataStructures(C);
99: /* Create vector x that is compatible with P */
100: MatCreateVecs(P,&x,NULL);
101: VecDuplicate(x,&v3);
102: VecDuplicate(x,&v4);
104: norm = 0.0;
105: for (i=0; i<10; i++) {
106: VecSetRandom(x,rdm);
107: MatMult(P,x,v1);
108: MatMult(A,v1,v2); /* v2 = A*P*x */
110: MatMultTranspose(P,v2,v3); /* v3 = Pt*A*P*x */
111: MatMult(C,x,v4); /* v3 = C*x */
112: VecAXPY(v4,none,v3);
113: VecNorm(v4,NORM_1,&norm_tmp);
114: VecNorm(v3,NORM_1,&norm_tmp1);
116: norm_tmp /= norm_tmp1;
117: if (norm_tmp > norm) norm = norm_tmp;
118: }
119: if (norm >= tol && !rank) {
120: PetscPrintf(PETSC_COMM_SELF,"Error: MatPtAP(), |v3 - v4|/|v3|: %g\n",(double)norm);
121: }
123: MatDestroy(&C);
124: MatDestroy(&A);
125: MatDestroy(&P);
126: VecDestroy(&v3);
127: VecDestroy(&v4);
128: VecDestroy(&x);
129: VecDestroy(&v1);
130: VecDestroy(&v2);
131: DMDestroy(&finedm);
132: DMDestroy(&coarsedm);
133: PetscRandomDestroy(&rdm);
134: PetscFinalize();
135: return ierr;
136: }
138: /*TEST
140: test:
141: args: -M 10 -N 10 -Z 10
142: output_file: output/ex89_1.out
144: test:
145: suffix: allatonce
146: nsize: 4
147: args: -M 10 -N 10 -Z 10 -matmaijptap_via allatonce
148: output_file: output/ex89_1.out
150: test:
151: suffix: allatonce_merged
152: nsize: 4
153: args: -M 10 -M 5 -M 10 -matmaijptap_via allatonce_merged
154: output_file: output/ex96_1.out
156: test:
157: suffix: allatonce_3D
158: nsize: 4
159: args: -M 10 -M 5 -M 10 -test_3D 1 -matmaijptap_via allatonce
160: output_file: output/ex96_1.out
162: test:
163: suffix: allatonce_merged_3D
164: nsize: 4
165: args: -M 10 -M 5 -M 10 -test_3D 1 -matmaijptap_via allatonce_merged
166: output_file: output/ex96_1.out
168: TEST*/