Actual source code: ex114.c
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];
15: PetscMPIInt size,rank;
16: PetscReal enorm;
18: PetscInitialize(&argc,&args,(char*)0,help);
19: MPI_Comm_size(PETSC_COMM_WORLD,&size);
20: MPI_Comm_rank(PETSC_COMM_WORLD,&rank);
21: PetscOptionsGetInt(NULL,NULL,"-testcase",&testcase,NULL);
23: MatCreate(PETSC_COMM_WORLD,&A);
24: if (testcase == 1) { /* proc[0] holds entire A and other processes have no entry */
25: if (rank == 0) {
26: MatSetSizes(A,M,N,PETSC_DECIDE,PETSC_DECIDE);
27: } else {
28: MatSetSizes(A,0,0,PETSC_DECIDE,PETSC_DECIDE);
29: }
30: testcase = 0;
31: } else {
32: MatSetSizes(A,PETSC_DECIDE,PETSC_DECIDE,M,N);
33: }
34: MatSetFromOptions(A);
35: MatSetUp(A);
37: if (rank == 0) { /* proc[0] sets matrix A */
38: for (j=0; j<N; j++) indices[j] = j;
39: switch (testcase) {
40: case 1: /* see testcast 0 */
41: break;
42: case 2:
43: row = 0;
44: values[0] = -2.0; values[1] = -2.0; values[2] = -2.0; values[3] = -4.0; values[4] = 1.0; values[5] = 1.0;
45: MatSetValues(A,1,&row,N,indices,values,INSERT_VALUES);
46: row = 2;
47: indices[0] = 0; indices[1] = 3; indices[2] = 5;
48: values[0] = -2.0; values[1] = -2.0; values[2] = -2.0;
49: MatSetValues(A,1,&row,3,indices,values,INSERT_VALUES);
50: row = 3;
51: indices[0] = 0; indices[1] = 1; indices[2] = 4;
52: values[0] = -2.0; values[1] = -2.0; values[2] = -2.0;
53: MatSetValues(A,1,&row,3,indices,values,INSERT_VALUES);
54: row = 4;
55: indices[0] = 0; indices[1] = 1; indices[2] = 2;
56: values[0] = -2.0; values[1] = -2.0; values[2] = -2.0;
57: MatSetValues(A,1,&row,3,indices,values,INSERT_VALUES);
58: break;
59: case 3:
60: row = 0;
61: values[0] = -2.0; values[1] = -2.0; values[2] = -2.0;
62: MatSetValues(A,1,&row,3,indices+1,values,INSERT_VALUES);
63: row = 1;
64: values[0] = -2.0; values[1] = -2.0; values[2] = -2.0;
65: MatSetValues(A,1,&row,3,indices,values,INSERT_VALUES);
66: row = 2;
67: values[0] = -2.0; values[1] = -2.0; values[2] = -2.0;
68: MatSetValues(A,1,&row,3,indices,values,INSERT_VALUES);
69: row = 3;
70: values[0] = -2.0; values[1] = -2.0; values[2] = -2.0; values[3] = -1.0;
71: MatSetValues(A,1,&row,4,indices,values,INSERT_VALUES);
72: row = 4;
73: values[0] = -2.0; values[1] = -2.0; values[2] = -2.0; values[3] = -1.0;
74: MatSetValues(A,1,&row,4,indices,values,INSERT_VALUES);
75: break;
77: default:
78: row = 0;
79: values[0] = -1.0; values[1] = 0.0; values[2] = 1.0; values[3] = 3.0; values[4] = 4.0; values[5] = -5.0;
80: MatSetValues(A,1,&row,N,indices,values,INSERT_VALUES);
81: row = 1;
82: MatSetValues(A,1,&row,3,indices,values,INSERT_VALUES);
83: row = 3;
84: MatSetValues(A,1,&row,1,indices+4,values+4,INSERT_VALUES);
85: row = 4;
86: MatSetValues(A,1,&row,2,indices+4,values+4,INSERT_VALUES);
87: }
88: }
89: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
90: MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
91: MatView(A,PETSC_VIEWER_STDOUT_WORLD);
93: MatGetLocalSize(A, &m,&n);
94: VecCreate(PETSC_COMM_WORLD,&min);
95: VecSetSizes(min,m,PETSC_DECIDE);
96: VecSetFromOptions(min);
97: VecDuplicate(min,&max);
98: VecDuplicate(min,&maxabs);
99: VecDuplicate(min,&e);
101: /* MatGetRowMax() */
102: PetscPrintf(PETSC_COMM_WORLD,"\n MatGetRowMax\n");
103: MatGetRowMax(A,max,NULL);
104: MatGetRowMax(A,max,imax);
105: VecView(max,PETSC_VIEWER_STDOUT_WORLD);
106: VecGetLocalSize(max,&n);
107: PetscIntView(n,imax,PETSC_VIEWER_STDOUT_WORLD);
109: /* MatGetRowMin() */
110: MatScale(A,-1.0);
111: MatView(A,PETSC_VIEWER_STDOUT_WORLD);
112: PetscPrintf(PETSC_COMM_WORLD,"\n MatGetRowMin\n");
113: MatGetRowMin(A,min,NULL);
114: MatGetRowMin(A,min,imin);
116: VecWAXPY(e,1.0,max,min); /* e = max + min */
117: VecNorm(e,NORM_INFINITY,&enorm);
119: for (j = 0; j < n; j++) {
121: }
123: /* MatGetRowMaxAbs() */
124: PetscPrintf(PETSC_COMM_WORLD,"\n MatGetRowMaxAbs\n");
125: MatGetRowMaxAbs(A,maxabs,NULL);
126: MatGetRowMaxAbs(A,maxabs,imaxabs);
127: VecView(maxabs,PETSC_VIEWER_STDOUT_WORLD);
128: PetscIntView(n,imaxabs,PETSC_VIEWER_STDOUT_WORLD);
130: /* MatGetRowMinAbs() */
131: PetscPrintf(PETSC_COMM_WORLD,"\n MatGetRowMinAbs\n");
132: MatGetRowMinAbs(A,min,NULL);
133: MatGetRowMinAbs(A,min,imin);
134: VecView(min,PETSC_VIEWER_STDOUT_WORLD);
135: PetscIntView(n,imin,PETSC_VIEWER_STDOUT_WORLD);
137: if (size == 1) {
138: /* Test MatGetRowMax, MatGetRowMin and MatGetRowMaxAbs for SeqDense and MPIBAIJ matrix */
139: Mat Adense;
140: Vec max_d,maxabs_d;
141: VecDuplicate(min,&max_d);
142: VecDuplicate(min,&maxabs_d);
144: MatScale(A,-1.0);
145: MatConvert(A,MATDENSE,MAT_INITIAL_MATRIX,&Adense);
146: PetscPrintf(PETSC_COMM_WORLD,"MatGetRowMax for seqdense matrix\n");
147: MatGetRowMax(Adense,max_d,imax);
149: VecWAXPY(e,-1.0,max,max_d); /* e = -max + max_d */
150: VecNorm(e,NORM_INFINITY,&enorm);
153: MatScale(Adense,-1.0);
154: PetscPrintf(PETSC_COMM_WORLD,"MatGetRowMin for seqdense matrix\n");
155: MatGetRowMin(Adense,min,imin);
157: VecWAXPY(e,1.0,max,min); /* e = max + min */
158: VecNorm(e,NORM_INFINITY,&enorm);
160: for (j = 0; j < n; j++) {
161: if (imin[j] != imax[j]) {
162: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"imin[%" PetscInt_FMT "] %" PetscInt_FMT " != imax %" PetscInt_FMT " for seqdense matrix",j,imin[j],imax[j]);
163: }
164: }
166: PetscPrintf(PETSC_COMM_WORLD,"MatGetRowMaxAbs for seqdense matrix\n");
167: MatGetRowMaxAbs(Adense,maxabs_d,imaxabs);
168: VecWAXPY(e,-1.0,maxabs,maxabs_d); /* e = -maxabs + maxabs_d */
169: VecNorm(e,NORM_INFINITY,&enorm);
172: MatDestroy(&Adense);
173: VecDestroy(&max_d);
174: VecDestroy(&maxabs_d);
175: }
177: { /* BAIJ matrix */
178: Mat B;
179: Vec maxabsB,maxabsB2;
180: PetscInt bs=2,*imaxabsB,*imaxabsB2,rstart,rend,cstart,cend,ncols,col,Brows[2],Bcols[2];
181: const PetscInt *cols;
182: const PetscScalar *vals,*vals2;
183: PetscScalar Bvals[4];
185: PetscMalloc2(M,&imaxabsB,bs*M,&imaxabsB2);
187: /* bs = 1 */
188: MatConvert(A,MATMPIBAIJ,MAT_INITIAL_MATRIX,&B);
189: VecDuplicate(min,&maxabsB);
190: MatGetRowMaxAbs(B,maxabsB,NULL);
191: MatGetRowMaxAbs(B,maxabsB,imaxabsB);
192: PetscPrintf(PETSC_COMM_WORLD,"\n MatGetRowMaxAbs for BAIJ matrix\n");
193: VecWAXPY(e,-1.0,maxabs,maxabsB); /* e = -maxabs + maxabsB */
194: VecNorm(e,NORM_INFINITY,&enorm);
197: for (j = 0; j < n; j++) {
199: }
200: MatDestroy(&B);
202: /* Test bs = 2: Create B with bs*bs block structure of A */
203: VecCreate(PETSC_COMM_WORLD,&maxabsB2);
204: VecSetSizes(maxabsB2,bs*m,PETSC_DECIDE);
205: VecSetFromOptions(maxabsB2);
207: MatGetOwnershipRange(A,&rstart,&rend);
208: MatGetOwnershipRangeColumn(A,&cstart,&cend);
209: MatCreate(PETSC_COMM_WORLD,&B);
210: MatSetSizes(B,bs*(rend-rstart),bs*(cend-cstart),PETSC_DECIDE,PETSC_DECIDE);
211: MatSetFromOptions(B);
212: MatSetUp(B);
214: for (row=rstart; row<rend; row++) {
215: MatGetRow(A,row,&ncols,&cols,&vals);
216: for (col=0; col<ncols; col++) {
217: for (j=0; j<bs; j++) {
218: Brows[j] = bs*row + j;
219: Bcols[j] = bs*cols[col]+j;
220: }
221: for (j=0; j<bs*bs; j++) Bvals[j] = vals[col];
222: MatSetValues(B,bs,Brows,bs,Bcols,Bvals,INSERT_VALUES);
223: }
224: MatRestoreRow(A,row,&ncols,&cols,&vals);
225: }
226: MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);
227: MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);
229: MatGetRowMaxAbs(B,maxabsB2,imaxabsB2);
231: /* Check maxabsB2 and imaxabsB2 */
232: VecGetArrayRead(maxabsB,&vals);
233: VecGetArrayRead(maxabsB2,&vals2);
234: for (row=0; row<m; row++) {
235: if (PetscAbsScalar(vals[row] - vals2[bs*row]) > PETSC_MACHINE_EPSILON)
236: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"row %" PetscInt_FMT " maxabsB != maxabsB2",row);
237: }
238: VecRestoreArrayRead(maxabsB,&vals);
239: VecRestoreArrayRead(maxabsB2,&vals2);
241: for (col=0; col<n; col++) {
242: if (imaxabsB[col] != imaxabsB2[bs*col]/bs)
243: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"col %" PetscInt_FMT " imaxabsB != imaxabsB2",col);
244: }
245: VecDestroy(&maxabsB);
246: MatDestroy(&B);
247: VecDestroy(&maxabsB2);
248: PetscFree2(imaxabsB,imaxabsB2);
249: }
251: VecDestroy(&min);
252: VecDestroy(&max);
253: VecDestroy(&maxabs);
254: VecDestroy(&e);
255: MatDestroy(&A);
256: PetscFinalize();
257: return 0;
258: }
260: /*TEST
262: test:
263: output_file: output/ex114.out
265: test:
266: suffix: 2
267: args: -testcase 1
268: output_file: output/ex114.out
270: test:
271: suffix: 3
272: args: -testcase 2
273: output_file: output/ex114_3.out
275: test:
276: suffix: 4
277: args: -testcase 3
278: output_file: output/ex114_4.out
280: test:
281: suffix: 5
282: nsize: 3
283: args: -testcase 0
284: output_file: output/ex114_5.out
286: test:
287: suffix: 6
288: nsize: 3
289: args: -testcase 1
290: output_file: output/ex114_6.out
292: test:
293: suffix: 7
294: nsize: 3
295: args: -testcase 2
296: output_file: output/ex114_7.out
298: test:
299: suffix: 8
300: nsize: 3
301: args: -testcase 3
302: output_file: output/ex114_8.out
304: TEST*/