Actual source code: ex114.c
petsc-3.14.6 2021-03-30
2: static char help[] = "Tests MatGetRowMax(), MatGetRowMin(), MatGetRowMaxAbs()\n";
4: #include <petscmat.h>
6: #define M 5
7: #define N 6
9: int main(int argc,char **args)
10: {
11: Mat A;
12: Vec min,max,maxabs,e;
13: PetscInt m,n,j,imin[M],imax[M],imaxabs[M],indices[N],row,testcase=0;
14: PetscScalar values[N];
16: PetscMPIInt size,rank;
17: PetscReal enorm;
19: PetscInitialize(&argc,&args,(char*)0,help);if (ierr) return ierr;
20: MPI_Comm_size(PETSC_COMM_WORLD,&size);
21: MPI_Comm_rank(PETSC_COMM_WORLD,&rank);
22: PetscOptionsGetInt(NULL,NULL,"-testcase",&testcase,NULL);
24: MatCreate(PETSC_COMM_WORLD,&A);
25: if (testcase == 1) { /* proc[0] holds entire A and other processes have no entry */
26: if (!rank) {
27: MatSetSizes(A,M,N,PETSC_DECIDE,PETSC_DECIDE);
28: } else {
29: MatSetSizes(A,0,0,PETSC_DECIDE,PETSC_DECIDE);
30: }
31: testcase = 0;
32: } else {
33: MatSetSizes(A,PETSC_DECIDE,PETSC_DECIDE,M,N);
34: }
35: MatSetFromOptions(A);
36: MatSetUp(A);
38: if (!rank) { /* proc[0] sets matrix A */
39: for (j=0; j<N; j++) indices[j] = j;
40: switch (testcase) {
41: case 1: /* see testcast 0 */
42: break;
43: case 2:
44: row = 0;
45: values[0] = -2.0; values[1] = -2.0; values[2] = -2.0; values[3] = -4.0; values[4] = 1.0; values[5] = 1.0;
46: MatSetValues(A,1,&row,N,indices,values,INSERT_VALUES);
47: row = 2;
48: indices[0] = 0; indices[1] = 3; indices[2] = 5;
49: values[0] = -2.0; values[1] = -2.0; values[2] = -2.0;
50: MatSetValues(A,1,&row,3,indices,values,INSERT_VALUES);
51: row = 3;
52: indices[0] = 0; indices[1] = 1; indices[2] = 4;
53: values[0] = -2.0; values[1] = -2.0; values[2] = -2.0;
54: MatSetValues(A,1,&row,3,indices,values,INSERT_VALUES);
55: row = 4;
56: indices[0] = 0; indices[1] = 1; indices[2] = 2;
57: values[0] = -2.0; values[1] = -2.0; values[2] = -2.0;
58: MatSetValues(A,1,&row,3,indices,values,INSERT_VALUES);
59: break;
60: case 3:
61: row = 0;
62: values[0] = -2.0; values[1] = -2.0; values[2] = -2.0;
63: MatSetValues(A,1,&row,3,indices+1,values,INSERT_VALUES);
64: row = 1;
65: values[0] = -2.0; values[1] = -2.0; values[2] = -2.0;
66: MatSetValues(A,1,&row,3,indices,values,INSERT_VALUES);
67: row = 2;
68: values[0] = -2.0; values[1] = -2.0; values[2] = -2.0;
69: MatSetValues(A,1,&row,3,indices,values,INSERT_VALUES);
70: row = 3;
71: values[0] = -2.0; values[1] = -2.0; values[2] = -2.0; values[3] = -1.0;
72: MatSetValues(A,1,&row,4,indices,values,INSERT_VALUES);
73: row = 4;
74: values[0] = -2.0; values[1] = -2.0; values[2] = -2.0; values[3] = -1.0;
75: MatSetValues(A,1,&row,4,indices,values,INSERT_VALUES);
76: break;
78: default:
79: row = 0;
80: values[0] = -1.0; values[1] = 0.0; values[2] = 1.0; values[3] = 3.0; values[4] = 4.0; values[5] = -5.0;
81: MatSetValues(A,1,&row,N,indices,values,INSERT_VALUES);
82: row = 1;
83: MatSetValues(A,1,&row,3,indices,values,INSERT_VALUES);
84: row = 3;
85: MatSetValues(A,1,&row,1,indices+4,values+4,INSERT_VALUES);
86: row = 4;
87: MatSetValues(A,1,&row,2,indices+4,values+4,INSERT_VALUES);
88: }
89: }
90: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
91: MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
92: MatView(A,PETSC_VIEWER_STDOUT_WORLD);
94: MatGetLocalSize(A, &m,&n);
95: VecCreate(PETSC_COMM_WORLD,&min);
96: VecSetSizes(min,m,PETSC_DECIDE);
97: VecSetFromOptions(min);
98: VecDuplicate(min,&max);
99: VecDuplicate(min,&maxabs);
100: VecDuplicate(min,&e);
102: /* MatGetRowMax() */
103: PetscPrintf(PETSC_COMM_WORLD,"\n MatGetRowMax\n");
104: MatGetRowMax(A,max,NULL);
105: MatGetRowMax(A,max,imax);
106: VecView(max,PETSC_VIEWER_STDOUT_WORLD);
107: VecGetLocalSize(max,&n);
108: PetscIntView(n,imax,PETSC_VIEWER_STDOUT_WORLD);
110: /* MatGetRowMin() */
111: MatScale(A,-1.0);
112: MatView(A,PETSC_VIEWER_STDOUT_WORLD);
113: PetscPrintf(PETSC_COMM_WORLD,"\n MatGetRowMin\n");
114: MatGetRowMin(A,min,NULL);
115: MatGetRowMin(A,min,imin);
117: VecWAXPY(e,1.0,max,min); /* e = max + min */
118: VecNorm(e,NORM_INFINITY,&enorm);
119: if (enorm > PETSC_MACHINE_EPSILON) SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_PLIB,"max+min > PETSC_MACHINE_EPSILON ");
120: for (j = 0; j < n; j++) {
121: if (imin[j] != imax[j]) SETERRQ3(PETSC_COMM_SELF,PETSC_ERR_PLIB,"imin[%D] %D != imax %D",j,imin[j],imax[j]);
122: }
124: /* MatGetRowMaxAbs() */
125: PetscPrintf(PETSC_COMM_WORLD,"\n MatGetRowMaxAbs\n");
126: MatGetRowMaxAbs(A,maxabs,NULL);
127: MatGetRowMaxAbs(A,maxabs,imaxabs);
128: VecView(maxabs,PETSC_VIEWER_STDOUT_WORLD);
129: PetscIntView(n,imaxabs,PETSC_VIEWER_STDOUT_WORLD);
131: /* MatGetRowMinAbs() */
132: PetscPrintf(PETSC_COMM_WORLD,"\n MatGetRowMinAbs\n");
133: MatGetRowMinAbs(A,min,NULL);
134: MatGetRowMinAbs(A,min,imin);
135: VecView(min,PETSC_VIEWER_STDOUT_WORLD);
136: PetscIntView(n,imin,PETSC_VIEWER_STDOUT_WORLD);
138: if (size == 1) {
139: /* Test MatGetRowMax, MatGetRowMin and MatGetRowMaxAbs for SeqDense and MPIBAIJ matrix */
140: Mat Adense;
141: Vec max_d,maxabs_d;
142: VecDuplicate(min,&max_d);
143: VecDuplicate(min,&maxabs_d);
145: MatScale(A,-1.0);
146: MatConvert(A,MATDENSE,MAT_INITIAL_MATRIX,&Adense);
147: PetscPrintf(PETSC_COMM_WORLD,"MatGetRowMax for seqdense matrix\n");
148: MatGetRowMax(Adense,max_d,imax);
150: VecWAXPY(e,-1.0,max,max_d); /* e = -max + max_d */
151: VecNorm(e,NORM_INFINITY,&enorm);
152: if (enorm > PETSC_MACHINE_EPSILON) SETERRQ1(PETSC_COMM_WORLD,PETSC_ERR_PLIB,"norm(-max + max_d) %g > PETSC_MACHINE_EPSILON",(double)enorm);
154: MatScale(Adense,-1.0);
155: PetscPrintf(PETSC_COMM_WORLD,"MatGetRowMin for seqdense matrix\n");
156: MatGetRowMin(Adense,min,imin);
158: VecWAXPY(e,1.0,max,min); /* e = max + min */
159: VecNorm(e,NORM_INFINITY,&enorm);
160: if (enorm > PETSC_MACHINE_EPSILON) SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_PLIB,"max+min > PETSC_MACHINE_EPSILON ");
161: for (j = 0; j < n; j++) {
162: if (imin[j] != imax[j]) {
163: SETERRQ3(PETSC_COMM_SELF,PETSC_ERR_PLIB,"imin[%D] %D != imax %D for seqdense matrix",j,imin[j],imax[j]);
164: }
165: }
167: PetscPrintf(PETSC_COMM_WORLD,"MatGetRowMaxAbs for seqdense matrix\n");
168: MatGetRowMaxAbs(Adense,maxabs_d,imaxabs);
169: VecWAXPY(e,-1.0,maxabs,maxabs_d); /* e = -maxabs + maxabs_d */
170: VecNorm(e,NORM_INFINITY,&enorm);
171: if (enorm > PETSC_MACHINE_EPSILON) SETERRQ1(PETSC_COMM_WORLD,PETSC_ERR_PLIB,"norm(-maxabs + maxabs_d) %g > PETSC_MACHINE_EPSILON",(double)enorm);
173: MatDestroy(&Adense);
174: VecDestroy(&max_d);
175: VecDestroy(&maxabs_d);
176: }
178: { /* BAIJ matrix */
179: Mat B;
180: Vec maxabsB,maxabsB2;
181: PetscInt bs=2,*imaxabsB,*imaxabsB2,rstart,rend,cstart,cend,ncols,col,Brows[2],Bcols[2];
182: const PetscInt *cols;
183: const PetscScalar *vals,*vals2;
184: PetscScalar Bvals[4];
186: PetscMalloc2(M,&imaxabsB,bs*M,&imaxabsB2);
188: /* bs = 1 */
189: MatConvert(A,MATMPIBAIJ,MAT_INITIAL_MATRIX,&B);
190: VecDuplicate(min,&maxabsB);
191: MatGetRowMaxAbs(B,maxabsB,NULL);
192: MatGetRowMaxAbs(B,maxabsB,imaxabsB);
193: PetscPrintf(PETSC_COMM_WORLD,"\n MatGetRowMaxAbs for BAIJ matrix\n");
194: VecWAXPY(e,-1.0,maxabs,maxabsB); /* e = -maxabs + maxabsB */
195: VecNorm(e,NORM_INFINITY,&enorm);
196: if (enorm > PETSC_MACHINE_EPSILON) SETERRQ1(PETSC_COMM_WORLD,PETSC_ERR_PLIB,"norm(-maxabs + maxabs_d) %g > PETSC_MACHINE_EPSILON",(double)enorm);
198: for (j = 0; j < n; j++) {
199: if (imaxabs[j] != imaxabsB[j]) SETERRQ3(PETSC_COMM_SELF,PETSC_ERR_PLIB,"imaxabs[%D] %D != imaxabsB %D",j,imin[j],imax[j]);
200: }
201: MatDestroy(&B);
203: /* Test bs = 2: Create B with bs*bs block structure of A */
204: VecCreate(PETSC_COMM_WORLD,&maxabsB2);
205: VecSetSizes(maxabsB2,bs*m,PETSC_DECIDE);
206: VecSetFromOptions(maxabsB2);
208: MatGetOwnershipRange(A,&rstart,&rend);
209: MatGetOwnershipRangeColumn(A,&cstart,&cend);
210: MatCreate(PETSC_COMM_WORLD,&B);
211: MatSetSizes(B,bs*(rend-rstart),bs*(cend-cstart),PETSC_DECIDE,PETSC_DECIDE);
212: MatSetFromOptions(B);
213: MatSetUp(B);
215: for (row=rstart; row<rend; row++) {
216: MatGetRow(A,row,&ncols,&cols,&vals);
217: for (col=0; col<ncols; col++) {
218: for (j=0; j<bs; j++) {
219: Brows[j] = bs*row + j;
220: Bcols[j] = bs*cols[col]+j;
221: }
222: for (j=0; j<bs*bs; j++) Bvals[j] = vals[col];
223: MatSetValues(B,bs,Brows,bs,Bcols,Bvals,INSERT_VALUES);
224: }
225: MatRestoreRow(A,row,&ncols,&cols,&vals);
226: }
227: MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);
228: MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);
230: MatGetRowMaxAbs(B,maxabsB2,imaxabsB2);
232: /* Check maxabsB2 and imaxabsB2 */
233: VecGetArrayRead(maxabsB,&vals);
234: VecGetArrayRead(maxabsB2,&vals2);
235: for (row=0; row<m; row++) {
236: if (PetscAbsScalar(vals[row] - vals2[bs*row]) > PETSC_MACHINE_EPSILON)
237: SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_PLIB,"row %D maxabsB != maxabsB2",row);
238: }
239: VecRestoreArrayRead(maxabsB,&vals);
240: VecRestoreArrayRead(maxabsB2,&vals2);
242: for (col=0; col<n; col++) {
243: if (imaxabsB[col] != imaxabsB2[bs*col]/bs)
244: SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_PLIB,"col %D imaxabsB != imaxabsB2",col);
245: }
246: VecDestroy(&maxabsB);
247: MatDestroy(&B);
248: VecDestroy(&maxabsB2);
249: PetscFree2(imaxabsB,imaxabsB2);
250: }
252: VecDestroy(&min);
253: VecDestroy(&max);
254: VecDestroy(&maxabs);
255: VecDestroy(&e);
256: MatDestroy(&A);
257: PetscFinalize();
258: return ierr;
259: }
261: /*TEST
263: test:
264: output_file: output/ex114.out
266: test:
267: suffix: 2
268: args: -testcase 1
269: output_file: output/ex114.out
271: test:
272: suffix: 3
273: args: -testcase 2
274: output_file: output/ex114_3.out
276: test:
277: suffix: 4
278: args: -testcase 3
279: output_file: output/ex114_4.out
281: test:
282: suffix: 5
283: nsize: 3
284: args: -testcase 0
285: output_file: output/ex114_5.out
287: test:
288: suffix: 6
289: nsize: 3
290: args: -testcase 1
291: output_file: output/ex114_6.out
293: test:
294: suffix: 7
295: nsize: 3
296: args: -testcase 2
297: output_file: output/ex114_7.out
299: test:
300: suffix: 8
301: nsize: 3
302: args: -testcase 3
303: output_file: output/ex114_8.out
305: TEST*/