Actual source code: ex75.c
petsc-3.13.6 2020-09-29
2: static char help[] = "Tests various routines in MatMPISBAIJ format.\n";
4: #include <petscmat.h>
6: int main(int argc,char **args)
7: {
8: Vec x,y,u,s1,s2;
9: Mat A,sA,sB;
10: PetscRandom rctx;
11: PetscReal r1,r2,rnorm,tol = PETSC_SQRT_MACHINE_EPSILON;
12: PetscScalar one=1.0, neg_one=-1.0, value[3], four=4.0,alpha=0.1;
13: PetscInt n,col[3],n1,block,row,i,j,i2,j2,Ii,J,rstart,rend,bs=1,mbs=16,d_nz=3,o_nz=3,prob=2;
15: PetscMPIInt size,rank;
16: PetscBool flg;
17: MatType type;
19: PetscInitialize(&argc,&args,(char*)0,help);if (ierr) return ierr;
20: PetscOptionsGetInt(NULL,NULL,"-mbs",&mbs,NULL);
21: PetscOptionsGetInt(NULL,NULL,"-bs",&bs,NULL);
23: MPI_Comm_rank(PETSC_COMM_WORLD,&rank);
24: MPI_Comm_size(PETSC_COMM_WORLD,&size);
26: n = mbs*bs;
28: /* Assemble MPISBAIJ matrix sA */
29: MatCreate(PETSC_COMM_WORLD,&sA);
30: MatSetSizes(sA,PETSC_DECIDE,PETSC_DECIDE,n,n);
31: MatSetType(sA,MATSBAIJ);
32: MatSetFromOptions(sA);
33: MatGetType(sA,&type);
34: MatMPISBAIJSetPreallocation(sA,bs,d_nz,NULL,o_nz,NULL);
35: MatSeqSBAIJSetPreallocation(sA,bs,d_nz,NULL);
36: MatSetOption(sA,MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_FALSE);
38: if (bs == 1) {
39: if (prob == 1) { /* tridiagonal matrix */
40: value[0] = -1.0; value[1] = 2.0; value[2] = -1.0;
41: for (i=1; i<n-1; i++) {
42: col[0] = i-1; col[1] = i; col[2] = i+1;
43: MatSetValues(sA,1,&i,3,col,value,INSERT_VALUES);
44: }
45: i = n - 1; col[0]=0; col[1] = n - 2; col[2] = n - 1;
46: value[0]= 0.1; value[1]=-1; value[2]=2;
47: MatSetValues(sA,1,&i,3,col,value,INSERT_VALUES);
49: i = 0; col[0] = 0; col[1] = 1; col[2]=n-1;
50: value[0] = 2.0; value[1] = -1.0; value[2]=0.1;
51: MatSetValues(sA,1,&i,3,col,value,INSERT_VALUES);
52: } else if (prob ==2) { /* matrix for the five point stencil */
53: n1 = (int) PetscSqrtReal((PetscReal)n);
54: if (n1*n1 != n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"n must be a perfect square of n1");
56: for (i=0; i<n1; i++) {
57: for (j=0; j<n1; j++) {
58: Ii = j + n1*i;
59: if (i>0) {J = Ii - n1; MatSetValues(sA,1,&Ii,1,&J,&neg_one,INSERT_VALUES);}
60: if (i<n1-1) {J = Ii + n1; MatSetValues(sA,1,&Ii,1,&J,&neg_one,INSERT_VALUES);}
61: if (j>0) {J = Ii - 1; MatSetValues(sA,1,&Ii,1,&J,&neg_one,INSERT_VALUES);}
62: if (j<n1-1) {J = Ii + 1; MatSetValues(sA,1,&Ii,1,&J,&neg_one,INSERT_VALUES);}
63: MatSetValues(sA,1,&Ii,1,&Ii,&four,INSERT_VALUES);
64: }
65: }
66: }
67: /* end of if (bs == 1) */
68: } else { /* bs > 1 */
69: for (block=0; block<n/bs; block++) {
70: /* diagonal blocks */
71: value[0] = -1.0; value[1] = 4.0; value[2] = -1.0;
72: for (i=1+block*bs; i<bs-1+block*bs; i++) {
73: col[0] = i-1; col[1] = i; col[2] = i+1;
74: MatSetValues(sA,1,&i,3,col,value,INSERT_VALUES);
75: }
76: i = bs - 1+block*bs; col[0] = bs - 2+block*bs; col[1] = bs - 1+block*bs;
77: value[0]=-1.0; value[1]=4.0;
78: MatSetValues(sA,1,&i,2,col,value,INSERT_VALUES);
80: i = 0+block*bs; col[0] = 0+block*bs; col[1] = 1+block*bs;
81: value[0]=4.0; value[1] = -1.0;
82: MatSetValues(sA,1,&i,2,col,value,INSERT_VALUES);
83: }
84: /* off-diagonal blocks */
85: value[0]=-1.0;
86: for (i=0; i<(n/bs-1)*bs; i++) {
87: col[0]=i+bs;
88: MatSetValues(sA,1,&i,1,col,value,INSERT_VALUES);
89: col[0]=i; row=i+bs;
90: MatSetValues(sA,1,&row,1,col,value,INSERT_VALUES);
91: }
92: }
93: MatAssemblyBegin(sA,MAT_FINAL_ASSEMBLY);
94: MatAssemblyEnd(sA,MAT_FINAL_ASSEMBLY);
96: /* Test MatView() */
97: MatCreateBAIJ(PETSC_COMM_WORLD,bs,PETSC_DECIDE,PETSC_DECIDE,n,n,d_nz,NULL,o_nz,NULL,&A);
98: MatSetOption(A,MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_FALSE);
100: if (bs == 1) {
101: if (prob == 1) { /* tridiagonal matrix */
102: value[0] = -1.0; value[1] = 2.0; value[2] = -1.0;
103: for (i=1; i<n-1; i++) {
104: col[0] = i-1; col[1] = i; col[2] = i+1;
105: MatSetValues(A,1,&i,3,col,value,INSERT_VALUES);
106: }
107: i = n - 1; col[0]=0; col[1] = n - 2; col[2] = n - 1;
108: value[0]= 0.1; value[1]=-1; value[2]=2;
109: MatSetValues(A,1,&i,3,col,value,INSERT_VALUES);
111: i = 0; col[0] = 0; col[1] = 1; col[2]=n-1;
112: value[0] = 2.0; value[1] = -1.0; value[2]=0.1;
113: MatSetValues(A,1,&i,3,col,value,INSERT_VALUES);
114: } else if (prob ==2) { /* matrix for the five point stencil */
115: n1 = (int) PetscSqrtReal((PetscReal)n);
116: for (i=0; i<n1; i++) {
117: for (j=0; j<n1; j++) {
118: Ii = j + n1*i;
119: if (i>0) {J = Ii - n1; MatSetValues(A,1,&Ii,1,&J,&neg_one,INSERT_VALUES);}
120: if (i<n1-1) {J = Ii + n1; MatSetValues(A,1,&Ii,1,&J,&neg_one,INSERT_VALUES);}
121: if (j>0) {J = Ii - 1; MatSetValues(A,1,&Ii,1,&J,&neg_one,INSERT_VALUES);}
122: if (j<n1-1) {J = Ii + 1; MatSetValues(A,1,&Ii,1,&J,&neg_one,INSERT_VALUES);}
123: MatSetValues(A,1,&Ii,1,&Ii,&four,INSERT_VALUES);
124: }
125: }
126: }
127: /* end of if (bs == 1) */
128: } else { /* bs > 1 */
129: for (block=0; block<n/bs; block++) {
130: /* diagonal blocks */
131: value[0] = -1.0; value[1] = 4.0; value[2] = -1.0;
132: for (i=1+block*bs; i<bs-1+block*bs; i++) {
133: col[0] = i-1; col[1] = i; col[2] = i+1;
134: MatSetValues(A,1,&i,3,col,value,INSERT_VALUES);
135: }
136: i = bs - 1+block*bs; col[0] = bs - 2+block*bs; col[1] = bs - 1+block*bs;
137: value[0]=-1.0; value[1]=4.0;
138: MatSetValues(A,1,&i,2,col,value,INSERT_VALUES);
140: i = 0+block*bs; col[0] = 0+block*bs; col[1] = 1+block*bs;
141: value[0]=4.0; value[1] = -1.0;
142: MatSetValues(A,1,&i,2,col,value,INSERT_VALUES);
143: }
144: /* off-diagonal blocks */
145: value[0]=-1.0;
146: for (i=0; i<(n/bs-1)*bs; i++) {
147: col[0]=i+bs;
148: MatSetValues(A,1,&i,1,col,value,INSERT_VALUES);
149: col[0]=i; row=i+bs;
150: MatSetValues(A,1,&row,1,col,value,INSERT_VALUES);
151: }
152: }
153: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
154: MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
156: /* Test MatGetSize(), MatGetLocalSize() */
157: MatGetSize(sA, &i,&j);
158: MatGetSize(A, &i2,&j2);
159: i -= i2; j -= j2;
160: if (i || j) {
161: PetscSynchronizedPrintf(PETSC_COMM_WORLD,"[%d], Error: MatGetSize()\n",rank);
162: PetscSynchronizedFlush(PETSC_COMM_WORLD,PETSC_STDOUT);
163: }
165: MatGetLocalSize(sA, &i,&j);
166: MatGetLocalSize(A, &i2,&j2);
167: i2 -= i; j2 -= j;
168: if (i2 || j2) {
169: PetscSynchronizedPrintf(PETSC_COMM_WORLD,"[%d], Error: MatGetLocalSize()\n",rank);
170: PetscSynchronizedFlush(PETSC_COMM_WORLD,PETSC_STDOUT);
171: }
173: /* vectors */
174: /*--------------------*/
175: /* i is obtained from MatGetLocalSize() */
176: VecCreate(PETSC_COMM_WORLD,&x);
177: VecSetSizes(x,i,PETSC_DECIDE);
178: VecSetFromOptions(x);
179: VecDuplicate(x,&y);
180: VecDuplicate(x,&u);
181: VecDuplicate(x,&s1);
182: VecDuplicate(x,&s2);
184: PetscRandomCreate(PETSC_COMM_WORLD,&rctx);
185: PetscRandomSetFromOptions(rctx);
186: VecSetRandom(x,rctx);
187: VecSet(u,one);
189: /* Test MatNorm() */
190: MatNorm(A,NORM_FROBENIUS,&r1);
191: MatNorm(sA,NORM_FROBENIUS,&r2);
192: rnorm = PetscAbsReal(r1-r2)/r2;
193: if (rnorm > tol && !rank) {
194: PetscPrintf(PETSC_COMM_SELF,"Error: MatNorm_FROBENIUS(), Anorm=%16.14e, sAnorm=%16.14e bs=%D\n",r1,r2,bs);
195: }
196: MatNorm(A,NORM_INFINITY,&r1);
197: MatNorm(sA,NORM_INFINITY,&r2);
198: rnorm = PetscAbsReal(r1-r2)/r2;
199: if (rnorm > tol && !rank) {
200: PetscPrintf(PETSC_COMM_WORLD,"Error: MatNorm_INFINITY(), Anorm=%16.14e, sAnorm=%16.14e bs=%D\n",r1,r2,bs);
201: }
202: MatNorm(A,NORM_1,&r1);
203: MatNorm(sA,NORM_1,&r2);
204: rnorm = PetscAbsReal(r1-r2)/r2;
205: if (rnorm > tol && !rank) {
206: PetscPrintf(PETSC_COMM_WORLD,"Error: MatNorm_1(), Anorm=%16.14e, sAnorm=%16.14e bs=%D\n",r1,r2,bs);
207: }
209: /* Test MatGetOwnershipRange() */
210: MatGetOwnershipRange(sA,&rstart,&rend);
211: MatGetOwnershipRange(A,&i2,&j2);
212: i2 -= rstart; j2 -= rend;
213: if (i2 || j2) {
214: PetscSynchronizedPrintf(PETSC_COMM_WORLD,"[%d], Error: MaGetOwnershipRange()\n",rank);
215: PetscSynchronizedFlush(PETSC_COMM_WORLD,PETSC_STDOUT);
216: }
218: /* Test MatDiagonalScale() */
219: MatDiagonalScale(A,x,x);
220: MatDiagonalScale(sA,x,x);
221: MatMultEqual(A,sA,10,&flg);
222: if (!flg) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_NOTSAMETYPE,"Error in MatDiagonalScale");
224: /* Test MatGetDiagonal(), MatScale() */
225: MatGetDiagonal(A,s1);
226: MatGetDiagonal(sA,s2);
227: VecNorm(s1,NORM_1,&r1);
228: VecNorm(s2,NORM_1,&r2);
229: r1 -= r2;
230: if (r1<-tol || r1>tol) {
231: PetscSynchronizedPrintf(PETSC_COMM_WORLD,"[%d], Error: MatDiagonalScale() or MatGetDiagonal(), r1=%g \n",rank,(double)r1);
232: PetscSynchronizedFlush(PETSC_COMM_WORLD,PETSC_STDOUT);
233: }
235: MatScale(A,alpha);
236: MatScale(sA,alpha);
238: /* Test MatGetRowMaxAbs() */
239: MatGetRowMaxAbs(A,s1,NULL);
240: MatGetRowMaxAbs(sA,s2,NULL);
242: VecNorm(s1,NORM_1,&r1);
243: VecNorm(s2,NORM_1,&r2);
244: r1 -= r2;
245: if (r1<-tol || r1>tol) {
246: PetscPrintf(PETSC_COMM_SELF,"Error: MatGetRowMaxAbs() \n");
247: }
249: /* Test MatMult(), MatMultAdd() */
250: MatMultEqual(A,sA,10,&flg);
251: if (!flg) {
252: PetscSynchronizedPrintf(PETSC_COMM_WORLD,"[%d], Error: MatMult() or MatScale()\n",rank);
253: PetscSynchronizedFlush(PETSC_COMM_WORLD,PETSC_STDOUT);
254: }
256: MatMultAddEqual(A,sA,10,&flg);
257: if (!flg) {
258: PetscSynchronizedPrintf(PETSC_COMM_WORLD,"[%d], Error: MatMultAdd()\n",rank);
259: PetscSynchronizedFlush(PETSC_COMM_WORLD,PETSC_STDOUT);
260: }
262: /* Test MatMultTranspose(), MatMultTransposeAdd() */
263: for (i=0; i<10; i++) {
264: VecSetRandom(x,rctx);
265: MatMultTranspose(A,x,s1);
266: MatMultTranspose(sA,x,s2);
267: VecNorm(s1,NORM_1,&r1);
268: VecNorm(s2,NORM_1,&r2);
269: r1 -= r2;
270: if (r1<-tol || r1>tol) {
271: PetscSynchronizedPrintf(PETSC_COMM_WORLD,"[%d], Error: MatMult() or MatScale(), err=%g\n",rank,(double)r1);
272: PetscSynchronizedFlush(PETSC_COMM_WORLD,PETSC_STDOUT);
273: }
274: }
275: for (i=0; i<10; i++) {
276: VecSetRandom(x,rctx);
277: VecSetRandom(y,rctx);
278: MatMultTransposeAdd(A,x,y,s1);
279: MatMultTransposeAdd(sA,x,y,s2);
280: VecNorm(s1,NORM_1,&r1);
281: VecNorm(s2,NORM_1,&r2);
282: r1 -= r2;
283: if (r1<-tol || r1>tol) {
284: PetscSynchronizedPrintf(PETSC_COMM_WORLD,"[%d], Error: MatMultAdd(), err=%g \n",rank,(double)r1);
285: PetscSynchronizedFlush(PETSC_COMM_WORLD,PETSC_STDOUT);
286: }
287: }
289: /* Test MatDuplicate() */
290: MatDuplicate(sA,MAT_COPY_VALUES,&sB);
291: MatEqual(sA,sB,&flg);
292: if (!flg) {
293: PetscPrintf(PETSC_COMM_WORLD," Error in MatDuplicate(), sA != sB \n");
294: }
295: MatMultEqual(sA,sB,5,&flg);
296: if (!flg) {
297: PetscSynchronizedPrintf(PETSC_COMM_WORLD,"[%d], Error: MatDuplicate() or MatMult()\n",rank);
298: PetscSynchronizedFlush(PETSC_COMM_WORLD,PETSC_STDOUT);
299: }
300: MatMultAddEqual(sA,sB,5,&flg);
301: if (!flg) {
302: PetscSynchronizedPrintf(PETSC_COMM_WORLD,"[%d], Error: MatDuplicate() or MatMultAdd(()\n",rank);
303: PetscSynchronizedFlush(PETSC_COMM_WORLD,PETSC_STDOUT);
304: }
305: MatDestroy(&sB);
306: VecDestroy(&u);
307: VecDestroy(&x);
308: VecDestroy(&y);
309: VecDestroy(&s1);
310: VecDestroy(&s2);
311: MatDestroy(&sA);
312: MatDestroy(&A);
313: PetscRandomDestroy(&rctx);
315: PetscFinalize();
316: return ierr;
317: }
319: /*TEST
321: test:
322: nsize: {{1 3}}
323: args: -bs {{1 2 3 5 7 8}} -mat_ignore_lower_triangular
325: TEST*/