Actual source code: ex75.c
petsc-3.14.6 2021-03-30
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=1;
15: PetscMPIInt size,rank;
16: PetscBool flg;
18: PetscInitialize(&argc,&args,(char*)0,help);if (ierr) return ierr;
19: PetscOptionsGetInt(NULL,NULL,"-mbs",&mbs,NULL);
20: PetscOptionsGetInt(NULL,NULL,"-bs",&bs,NULL);
21: PetscOptionsGetInt(NULL,NULL,"-prob",&prob,NULL);
23: MPI_Comm_rank(PETSC_COMM_WORLD,&rank);
24: MPI_Comm_size(PETSC_COMM_WORLD,&size);
26: /* Create a BAIJ matrix A */
27: n = mbs*bs;
28: MatCreate(PETSC_COMM_WORLD,&A);
29: MatSetSizes(A,PETSC_DECIDE,PETSC_DECIDE,n,n);
30: MatSetType(A,MATBAIJ);
31: MatSetFromOptions(A);
32: MatMPIBAIJSetPreallocation(A,bs,d_nz,NULL,o_nz,NULL);
33: MatSeqBAIJSetPreallocation(A,bs,d_nz,NULL);
34: MatSetOption(A,MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_FALSE);
36: if (bs == 1) {
37: if (prob == 1) { /* tridiagonal matrix */
38: value[0] = -1.0; value[1] = 0.0; value[2] = -1.0;
39: for (i=1; i<n-1; i++) {
40: col[0] = i-1; col[1] = i; col[2] = i+1;
41: MatSetValues(A,1,&i,3,col,value,INSERT_VALUES);
42: }
43: i = n - 1; col[0]=0; col[1] = n - 2; col[2] = n - 1;
44: value[0]= 0.1; value[1]=-1.0; value[2]=0.0;
45: MatSetValues(A,1,&i,3,col,value,INSERT_VALUES);
47: i = 0; col[0] = 0; col[1] = 1; col[2]=n-1;
48: value[0] = 0.0; value[1] = -1.0; value[2]=0.1;
49: MatSetValues(A,1,&i,3,col,value,INSERT_VALUES);
50: } else if (prob ==2) { /* matrix for the five point stencil */
51: n1 = (int) PetscSqrtReal((PetscReal)n);
52: for (i=0; i<n1; i++) {
53: for (j=0; j<n1; j++) {
54: Ii = j + n1*i;
55: if (i>0) {J = Ii - n1; MatSetValues(A,1,&Ii,1,&J,&neg_one,INSERT_VALUES);}
56: if (i<n1-1) {J = Ii + n1; MatSetValues(A,1,&Ii,1,&J,&neg_one,INSERT_VALUES);}
57: if (j>0) {J = Ii - 1; MatSetValues(A,1,&Ii,1,&J,&neg_one,INSERT_VALUES);}
58: if (j<n1-1) {J = Ii + 1; MatSetValues(A,1,&Ii,1,&J,&neg_one,INSERT_VALUES);}
59: MatSetValues(A,1,&Ii,1,&Ii,&four,INSERT_VALUES);
60: }
61: }
62: }
63: /* end of if (bs == 1) */
64: } else { /* bs > 1 */
65: for (block=0; block<n/bs; block++) {
66: /* diagonal blocks */
67: value[0] = -1.0; value[1] = 4.0; value[2] = -1.0;
68: for (i=1+block*bs; i<bs-1+block*bs; i++) {
69: col[0] = i-1; col[1] = i; col[2] = i+1;
70: MatSetValues(A,1,&i,3,col,value,INSERT_VALUES);
71: }
72: i = bs - 1+block*bs; col[0] = bs - 2+block*bs; col[1] = bs - 1+block*bs;
73: value[0]=-1.0; value[1]=4.0;
74: MatSetValues(A,1,&i,2,col,value,INSERT_VALUES);
76: i = 0+block*bs; col[0] = 0+block*bs; col[1] = 1+block*bs;
77: value[0]=4.0; value[1] = -1.0;
78: MatSetValues(A,1,&i,2,col,value,INSERT_VALUES);
79: }
80: /* off-diagonal blocks */
81: value[0]=-1.0;
82: for (i=0; i<(n/bs-1)*bs; i++) {
83: col[0]=i+bs;
84: MatSetValues(A,1,&i,1,col,value,INSERT_VALUES);
85: col[0]=i; row=i+bs;
86: MatSetValues(A,1,&row,1,col,value,INSERT_VALUES);
87: }
88: }
89: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
90: MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
91: MatSetOption(A,MAT_SYMMETRIC,PETSC_TRUE);
93: /* Get SBAIJ matrix sA from A */
94: MatConvert(A,MATSBAIJ,MAT_INITIAL_MATRIX,&sA);
96: /* Test MatGetSize(), MatGetLocalSize() */
97: MatGetSize(sA, &i,&j);
98: MatGetSize(A, &i2,&j2);
99: i -= i2; j -= j2;
100: if (i || j) {
101: PetscSynchronizedPrintf(PETSC_COMM_WORLD,"[%d], Error: MatGetSize()\n",rank);
102: PetscSynchronizedFlush(PETSC_COMM_WORLD,PETSC_STDOUT);
103: }
105: MatGetLocalSize(sA, &i,&j);
106: MatGetLocalSize(A, &i2,&j2);
107: i2 -= i; j2 -= j;
108: if (i2 || j2) {
109: PetscSynchronizedPrintf(PETSC_COMM_WORLD,"[%d], Error: MatGetLocalSize()\n",rank);
110: PetscSynchronizedFlush(PETSC_COMM_WORLD,PETSC_STDOUT);
111: }
113: /* vectors */
114: /*--------------------*/
115: /* i is obtained from MatGetLocalSize() */
116: VecCreate(PETSC_COMM_WORLD,&x);
117: VecSetSizes(x,i,PETSC_DECIDE);
118: VecSetFromOptions(x);
119: VecDuplicate(x,&y);
120: VecDuplicate(x,&u);
121: VecDuplicate(x,&s1);
122: VecDuplicate(x,&s2);
124: PetscRandomCreate(PETSC_COMM_WORLD,&rctx);
125: PetscRandomSetFromOptions(rctx);
126: VecSetRandom(x,rctx);
127: VecSet(u,one);
129: /* Test MatNorm() */
130: MatNorm(A,NORM_FROBENIUS,&r1);
131: MatNorm(sA,NORM_FROBENIUS,&r2);
132: rnorm = PetscAbsReal(r1-r2)/r2;
133: if (rnorm > tol && !rank) {
134: PetscPrintf(PETSC_COMM_SELF,"Error: MatNorm_FROBENIUS(), Anorm=%16.14e, sAnorm=%16.14e bs=%D\n",r1,r2,bs);
135: }
136: MatNorm(A,NORM_INFINITY,&r1);
137: MatNorm(sA,NORM_INFINITY,&r2);
138: rnorm = PetscAbsReal(r1-r2)/r2;
139: if (rnorm > tol && !rank) {
140: PetscPrintf(PETSC_COMM_WORLD,"Error: MatNorm_INFINITY(), Anorm=%16.14e, sAnorm=%16.14e bs=%D\n",r1,r2,bs);
141: }
142: MatNorm(A,NORM_1,&r1);
143: MatNorm(sA,NORM_1,&r2);
144: rnorm = PetscAbsReal(r1-r2)/r2;
145: if (rnorm > tol && !rank) {
146: PetscPrintf(PETSC_COMM_WORLD,"Error: MatNorm_1(), Anorm=%16.14e, sAnorm=%16.14e bs=%D\n",r1,r2,bs);
147: }
149: /* Test MatGetOwnershipRange() */
150: MatGetOwnershipRange(sA,&rstart,&rend);
151: MatGetOwnershipRange(A,&i2,&j2);
152: i2 -= rstart; j2 -= rend;
153: if (i2 || j2) {
154: PetscSynchronizedPrintf(PETSC_COMM_WORLD,"[%d], Error: MaGetOwnershipRange()\n",rank);
155: PetscSynchronizedFlush(PETSC_COMM_WORLD,PETSC_STDOUT);
156: }
158: /* Test MatDiagonalScale() */
159: MatDiagonalScale(A,x,x);
160: MatDiagonalScale(sA,x,x);
161: MatMultEqual(A,sA,10,&flg);
162: if (!flg) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_NOTSAMETYPE,"Error in MatDiagonalScale");
164: /* Test MatGetDiagonal(), MatScale() */
165: MatGetDiagonal(A,s1);
166: MatGetDiagonal(sA,s2);
167: VecNorm(s1,NORM_1,&r1);
168: VecNorm(s2,NORM_1,&r2);
169: r1 -= r2;
170: if (r1<-tol || r1>tol) {
171: PetscSynchronizedPrintf(PETSC_COMM_WORLD,"[%d], Error: MatDiagonalScale() or MatGetDiagonal(), r1=%g \n",rank,(double)r1);
172: PetscSynchronizedFlush(PETSC_COMM_WORLD,PETSC_STDOUT);
173: }
175: MatScale(A,alpha);
176: MatScale(sA,alpha);
178: /* Test MatGetRowMaxAbs() */
179: MatGetRowMaxAbs(A,s1,NULL);
180: MatGetRowMaxAbs(sA,s2,NULL);
182: VecNorm(s1,NORM_1,&r1);
183: VecNorm(s2,NORM_1,&r2);
184: r1 -= r2;
185: if (r1<-tol || r1>tol) {
186: PetscPrintf(PETSC_COMM_SELF,"Error: MatGetRowMaxAbs() \n");
187: }
189: /* Test MatMult(), MatMultAdd() */
190: MatMultEqual(A,sA,10,&flg);
191: if (!flg) {
192: PetscSynchronizedPrintf(PETSC_COMM_WORLD,"[%d], Error: MatMult() or MatScale()\n",rank);
193: PetscSynchronizedFlush(PETSC_COMM_WORLD,PETSC_STDOUT);
194: }
196: MatMultAddEqual(A,sA,10,&flg);
197: if (!flg) {
198: PetscSynchronizedPrintf(PETSC_COMM_WORLD,"[%d], Error: MatMultAdd()\n",rank);
199: PetscSynchronizedFlush(PETSC_COMM_WORLD,PETSC_STDOUT);
200: }
202: /* Test MatMultTranspose(), MatMultTransposeAdd() */
203: for (i=0; i<10; i++) {
204: VecSetRandom(x,rctx);
205: MatMultTranspose(A,x,s1);
206: MatMultTranspose(sA,x,s2);
207: VecNorm(s1,NORM_1,&r1);
208: VecNorm(s2,NORM_1,&r2);
209: r1 -= r2;
210: if (r1<-tol || r1>tol) {
211: PetscSynchronizedPrintf(PETSC_COMM_WORLD,"[%d], Error: MatMult() or MatScale(), err=%g\n",rank,(double)r1);
212: PetscSynchronizedFlush(PETSC_COMM_WORLD,PETSC_STDOUT);
213: }
214: }
215: for (i=0; i<10; i++) {
216: VecSetRandom(x,rctx);
217: VecSetRandom(y,rctx);
218: MatMultTransposeAdd(A,x,y,s1);
219: MatMultTransposeAdd(sA,x,y,s2);
220: VecNorm(s1,NORM_1,&r1);
221: VecNorm(s2,NORM_1,&r2);
222: r1 -= r2;
223: if (r1<-tol || r1>tol) {
224: PetscSynchronizedPrintf(PETSC_COMM_WORLD,"[%d], Error: MatMultAdd(), err=%g \n",rank,(double)r1);
225: PetscSynchronizedFlush(PETSC_COMM_WORLD,PETSC_STDOUT);
226: }
227: }
229: /* Test MatDuplicate() */
230: MatDuplicate(sA,MAT_COPY_VALUES,&sB);
231: MatEqual(sA,sB,&flg);
232: if (!flg) {
233: PetscPrintf(PETSC_COMM_WORLD," Error in MatDuplicate(), sA != sB \n");
234: }
235: MatMultEqual(sA,sB,5,&flg);
236: if (!flg) {
237: PetscSynchronizedPrintf(PETSC_COMM_WORLD,"[%d], Error: MatDuplicate() or MatMult()\n",rank);
238: PetscSynchronizedFlush(PETSC_COMM_WORLD,PETSC_STDOUT);
239: }
240: MatMultAddEqual(sA,sB,5,&flg);
241: if (!flg) {
242: PetscSynchronizedPrintf(PETSC_COMM_WORLD,"[%d], Error: MatDuplicate() or MatMultAdd(()\n",rank);
243: PetscSynchronizedFlush(PETSC_COMM_WORLD,PETSC_STDOUT);
244: }
245: MatDestroy(&sB);
246: VecDestroy(&u);
247: VecDestroy(&x);
248: VecDestroy(&y);
249: VecDestroy(&s1);
250: VecDestroy(&s2);
251: MatDestroy(&sA);
252: MatDestroy(&A);
253: PetscRandomDestroy(&rctx);
254: PetscFinalize();
255: return ierr;
256: }
258: /*TEST
260: test:
261: nsize: {{1 3}}
262: args: -bs {{1 2 3 5 7 8}} -mat_ignore_lower_triangular -prob {{1 2}}
264: TEST*/