xref: /petsc/src/mat/tests/ex114.c (revision 030f984af8d8bb4c203755d35bded3c05b3d83ce) 
1 
2 static char help[] = "Tests MatGetRowMax(), MatGetRowMin(), MatGetRowMaxAbs()\n";
3 
4 #include <petscmat.h>
5 
6 #define M 5
7 #define N 6
8 
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   PetscErrorCode ierr;
16   PetscMPIInt    size,rank;
17   PetscReal      enorm;
18 
19   ierr = PetscInitialize(&argc,&args,(char*)0,help);if (ierr) return ierr;
20   ierr = MPI_Comm_size(PETSC_COMM_WORLD,&size);CHKERRMPI(ierr);
21   ierr = MPI_Comm_rank(PETSC_COMM_WORLD,&rank);CHKERRMPI(ierr);
22   ierr = PetscOptionsGetInt(NULL,NULL,"-testcase",&testcase,NULL);CHKERRQ(ierr);
23 
24   ierr = MatCreate(PETSC_COMM_WORLD,&A);CHKERRQ(ierr);
25   if (testcase == 1) { /* proc[0] holds entire A and other processes have no entry */
26     if (!rank) {
27       ierr = MatSetSizes(A,M,N,PETSC_DECIDE,PETSC_DECIDE);CHKERRQ(ierr);
28     } else {
29       ierr = MatSetSizes(A,0,0,PETSC_DECIDE,PETSC_DECIDE);CHKERRQ(ierr);
30     }
31     testcase = 0;
32   } else {
33     ierr = MatSetSizes(A,PETSC_DECIDE,PETSC_DECIDE,M,N);CHKERRQ(ierr);
34   }
35   ierr = MatSetFromOptions(A);CHKERRQ(ierr);
36   ierr = MatSetUp(A);CHKERRQ(ierr);
37 
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       ierr = MatSetValues(A,1,&row,N,indices,values,INSERT_VALUES);CHKERRQ(ierr);
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       ierr = MatSetValues(A,1,&row,3,indices,values,INSERT_VALUES);CHKERRQ(ierr);
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       ierr = MatSetValues(A,1,&row,3,indices,values,INSERT_VALUES);CHKERRQ(ierr);
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       ierr = MatSetValues(A,1,&row,3,indices,values,INSERT_VALUES);CHKERRQ(ierr);
59       break;
60     case 3:
61       row = 0;
62       values[0]  = -2.0; values[1] = -2.0; values[2] = -2.0;
63       ierr = MatSetValues(A,1,&row,3,indices+1,values,INSERT_VALUES);CHKERRQ(ierr);
64       row = 1;
65       values[0]  = -2.0; values[1] = -2.0; values[2] = -2.0;
66       ierr = MatSetValues(A,1,&row,3,indices,values,INSERT_VALUES);CHKERRQ(ierr);
67       row = 2;
68       values[0]  = -2.0; values[1] = -2.0; values[2]  = -2.0;
69       ierr = MatSetValues(A,1,&row,3,indices,values,INSERT_VALUES);CHKERRQ(ierr);
70       row = 3;
71       values[0]  = -2.0; values[1] = -2.0; values[2]  = -2.0; values[3] = -1.0;
72       ierr = MatSetValues(A,1,&row,4,indices,values,INSERT_VALUES);CHKERRQ(ierr);
73       row = 4;
74       values[0]  = -2.0; values[1] = -2.0; values[2]  = -2.0; values[3] = -1.0;
75       ierr = MatSetValues(A,1,&row,4,indices,values,INSERT_VALUES);CHKERRQ(ierr);
76       break;
77 
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       ierr = MatSetValues(A,1,&row,N,indices,values,INSERT_VALUES);CHKERRQ(ierr);
82       row  = 1;
83       ierr = MatSetValues(A,1,&row,3,indices,values,INSERT_VALUES);CHKERRQ(ierr);
84       row  = 3;
85       ierr = MatSetValues(A,1,&row,1,indices+4,values+4,INSERT_VALUES);CHKERRQ(ierr);
86       row  = 4;
87       ierr = MatSetValues(A,1,&row,2,indices+4,values+4,INSERT_VALUES);CHKERRQ(ierr);
88     }
89   }
90   ierr = MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
91   ierr = MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
92   ierr = MatView(A,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);
93 
94   ierr = MatGetLocalSize(A, &m,&n);CHKERRQ(ierr);
95   ierr = VecCreate(PETSC_COMM_WORLD,&min);CHKERRQ(ierr);
96   ierr = VecSetSizes(min,m,PETSC_DECIDE);CHKERRQ(ierr);
97   ierr = VecSetFromOptions(min);CHKERRQ(ierr);
98   ierr = VecDuplicate(min,&max);CHKERRQ(ierr);
99   ierr = VecDuplicate(min,&maxabs);CHKERRQ(ierr);
100   ierr = VecDuplicate(min,&e);CHKERRQ(ierr);
101 
102   /* MatGetRowMax() */
103   ierr = PetscPrintf(PETSC_COMM_WORLD,"\n MatGetRowMax\n");CHKERRQ(ierr);
104   ierr = MatGetRowMax(A,max,NULL);CHKERRQ(ierr);
105   ierr = MatGetRowMax(A,max,imax);CHKERRQ(ierr);
106   ierr = VecView(max,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);
107   ierr = VecGetLocalSize(max,&n);CHKERRQ(ierr);
108   ierr = PetscIntView(n,imax,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);
109 
110   /* MatGetRowMin() */
111   ierr = MatScale(A,-1.0);CHKERRQ(ierr);
112   ierr = MatView(A,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);
113   ierr = PetscPrintf(PETSC_COMM_WORLD,"\n MatGetRowMin\n");CHKERRQ(ierr);
114   ierr = MatGetRowMin(A,min,NULL);CHKERRQ(ierr);
115   ierr = MatGetRowMin(A,min,imin);CHKERRQ(ierr);
116 
117   ierr = VecWAXPY(e,1.0,max,min);CHKERRQ(ierr); /* e = max + min */
118   ierr = VecNorm(e,NORM_INFINITY,&enorm);CHKERRQ(ierr);
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   }
123 
124   /* MatGetRowMaxAbs() */
125   ierr = PetscPrintf(PETSC_COMM_WORLD,"\n MatGetRowMaxAbs\n");CHKERRQ(ierr);
126   ierr = MatGetRowMaxAbs(A,maxabs,NULL);CHKERRQ(ierr);
127   ierr = MatGetRowMaxAbs(A,maxabs,imaxabs);CHKERRQ(ierr);
128   ierr = VecView(maxabs,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);
129   ierr = PetscIntView(n,imaxabs,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);
130 
131   /* MatGetRowMinAbs() */
132   ierr = PetscPrintf(PETSC_COMM_WORLD,"\n MatGetRowMinAbs\n");CHKERRQ(ierr);
133   ierr = MatGetRowMinAbs(A,min,NULL);CHKERRQ(ierr);
134   ierr = MatGetRowMinAbs(A,min,imin);CHKERRQ(ierr);
135   ierr = VecView(min,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);
136   ierr = PetscIntView(n,imin,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);
137 
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     ierr = VecDuplicate(min,&max_d);CHKERRQ(ierr);
143     ierr = VecDuplicate(min,&maxabs_d);CHKERRQ(ierr);
144 
145     ierr = MatScale(A,-1.0);CHKERRQ(ierr);
146     ierr = MatConvert(A,MATDENSE,MAT_INITIAL_MATRIX,&Adense);CHKERRQ(ierr);
147     ierr = PetscPrintf(PETSC_COMM_WORLD,"MatGetRowMax for seqdense matrix\n");CHKERRQ(ierr);
148     ierr = MatGetRowMax(Adense,max_d,imax);CHKERRQ(ierr);
149 
150     ierr = VecWAXPY(e,-1.0,max,max_d);CHKERRQ(ierr); /* e = -max + max_d */
151     ierr = VecNorm(e,NORM_INFINITY,&enorm);CHKERRQ(ierr);
152     if (enorm > PETSC_MACHINE_EPSILON) SETERRQ1(PETSC_COMM_WORLD,PETSC_ERR_PLIB,"norm(-max + max_d) %g > PETSC_MACHINE_EPSILON",(double)enorm);
153 
154     ierr = MatScale(Adense,-1.0);CHKERRQ(ierr);
155     ierr = PetscPrintf(PETSC_COMM_WORLD,"MatGetRowMin for seqdense matrix\n");CHKERRQ(ierr);
156     ierr = MatGetRowMin(Adense,min,imin);CHKERRQ(ierr);
157 
158     ierr = VecWAXPY(e,1.0,max,min);CHKERRQ(ierr); /* e = max + min */
159     ierr = VecNorm(e,NORM_INFINITY,&enorm);CHKERRQ(ierr);
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     }
166 
167     ierr = PetscPrintf(PETSC_COMM_WORLD,"MatGetRowMaxAbs for seqdense matrix\n");CHKERRQ(ierr);
168     ierr = MatGetRowMaxAbs(Adense,maxabs_d,imaxabs);CHKERRQ(ierr);
169     ierr = VecWAXPY(e,-1.0,maxabs,maxabs_d);CHKERRQ(ierr); /* e = -maxabs + maxabs_d */
170     ierr = VecNorm(e,NORM_INFINITY,&enorm);CHKERRQ(ierr);
171     if (enorm > PETSC_MACHINE_EPSILON) SETERRQ1(PETSC_COMM_WORLD,PETSC_ERR_PLIB,"norm(-maxabs + maxabs_d) %g > PETSC_MACHINE_EPSILON",(double)enorm);
172 
173     ierr = MatDestroy(&Adense);CHKERRQ(ierr);
174     ierr = VecDestroy(&max_d);CHKERRQ(ierr);
175     ierr = VecDestroy(&maxabs_d);CHKERRQ(ierr);
176   }
177 
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];
185 
186     ierr = PetscMalloc2(M,&imaxabsB,bs*M,&imaxabsB2);CHKERRQ(ierr);
187 
188     /* bs = 1 */
189     ierr = MatConvert(A,MATMPIBAIJ,MAT_INITIAL_MATRIX,&B);CHKERRQ(ierr);
190     ierr = VecDuplicate(min,&maxabsB);CHKERRQ(ierr);
191     ierr = MatGetRowMaxAbs(B,maxabsB,NULL);CHKERRQ(ierr);
192     ierr = MatGetRowMaxAbs(B,maxabsB,imaxabsB);CHKERRQ(ierr);
193     ierr = PetscPrintf(PETSC_COMM_WORLD,"\n MatGetRowMaxAbs for BAIJ matrix\n");CHKERRQ(ierr);
194     ierr = VecWAXPY(e,-1.0,maxabs,maxabsB);CHKERRQ(ierr); /* e = -maxabs + maxabsB */
195     ierr = VecNorm(e,NORM_INFINITY,&enorm);CHKERRQ(ierr);
196     if (enorm > PETSC_MACHINE_EPSILON) SETERRQ1(PETSC_COMM_WORLD,PETSC_ERR_PLIB,"norm(-maxabs + maxabs_d) %g > PETSC_MACHINE_EPSILON",(double)enorm);
197 
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     ierr = MatDestroy(&B);CHKERRQ(ierr);
202 
203     /* Test bs = 2: Create B with bs*bs block structure of A */
204     ierr = VecCreate(PETSC_COMM_WORLD,&maxabsB2);CHKERRQ(ierr);
205     ierr = VecSetSizes(maxabsB2,bs*m,PETSC_DECIDE);CHKERRQ(ierr);
206     ierr = VecSetFromOptions(maxabsB2);CHKERRQ(ierr);
207 
208     ierr = MatGetOwnershipRange(A,&rstart,&rend);CHKERRQ(ierr);
209     ierr = MatGetOwnershipRangeColumn(A,&cstart,&cend);CHKERRQ(ierr);
210     ierr = MatCreate(PETSC_COMM_WORLD,&B);CHKERRQ(ierr);
211     ierr = MatSetSizes(B,bs*(rend-rstart),bs*(cend-cstart),PETSC_DECIDE,PETSC_DECIDE);CHKERRQ(ierr);
212     ierr = MatSetFromOptions(B);CHKERRQ(ierr);
213     ierr = MatSetUp(B);CHKERRQ(ierr);
214 
215     for (row=rstart; row<rend; row++) {
216       ierr = MatGetRow(A,row,&ncols,&cols,&vals);CHKERRQ(ierr);
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         ierr = MatSetValues(B,bs,Brows,bs,Bcols,Bvals,INSERT_VALUES);CHKERRQ(ierr);
224       }
225       ierr = MatRestoreRow(A,row,&ncols,&cols,&vals);CHKERRQ(ierr);
226     }
227     ierr = MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
228     ierr = MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
229 
230     ierr = MatGetRowMaxAbs(B,maxabsB2,imaxabsB2);CHKERRQ(ierr);
231 
232     /* Check maxabsB2 and imaxabsB2 */
233     ierr = VecGetArrayRead(maxabsB,&vals);CHKERRQ(ierr);
234     ierr = VecGetArrayRead(maxabsB2,&vals2);CHKERRQ(ierr);
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     ierr = VecRestoreArrayRead(maxabsB,&vals);CHKERRQ(ierr);
240     ierr = VecRestoreArrayRead(maxabsB2,&vals2);CHKERRQ(ierr);
241 
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     ierr = VecDestroy(&maxabsB);CHKERRQ(ierr);
247     ierr = MatDestroy(&B);CHKERRQ(ierr);
248     ierr = VecDestroy(&maxabsB2);CHKERRQ(ierr);
249     ierr = PetscFree2(imaxabsB,imaxabsB2);CHKERRQ(ierr);
250   }
251 
252   ierr = VecDestroy(&min);CHKERRQ(ierr);
253   ierr = VecDestroy(&max);CHKERRQ(ierr);
254   ierr = VecDestroy(&maxabs);CHKERRQ(ierr);
255   ierr = VecDestroy(&e);CHKERRQ(ierr);
256   ierr = MatDestroy(&A);CHKERRQ(ierr);
257   ierr = PetscFinalize();
258   return ierr;
259 }
260 
261 /*TEST
262 
263    test:
264       output_file: output/ex114.out
265 
266    test:
267       suffix: 2
268       args: -testcase 1
269       output_file: output/ex114.out
270 
271    test:
272       suffix: 3
273       args: -testcase 2
274       output_file: output/ex114_3.out
275 
276    test:
277       suffix: 4
278       args: -testcase 3
279       output_file: output/ex114_4.out
280 
281    test:
282       suffix: 5
283       nsize: 3
284       args: -testcase 0
285       output_file: output/ex114_5.out
286 
287    test:
288       suffix: 6
289       nsize: 3
290       args: -testcase 1
291       output_file: output/ex114_6.out
292 
293    test:
294       suffix: 7
295       nsize: 3
296       args: -testcase 2
297       output_file: output/ex114_7.out
298 
299    test:
300       suffix: 8
301       nsize: 3
302       args: -testcase 3
303       output_file: output/ex114_8.out
304 
305 TEST*/
306