Actual source code: ex74.c
petsc-3.13.6 2020-09-29
2: static char help[] = "Tests the various sequential routines in MATSEQSBAIJ format.\n";
4: #include <petscmat.h>
6: int main(int argc,char **args)
7: {
8: PetscMPIInt size;
10: Vec x,y,b,s1,s2;
11: Mat A; /* linear system matrix */
12: Mat sA,sB,sFactor,B,C; /* symmetric matrices */
13: PetscInt n,mbs=16,bs=1,nz=3,prob=1,i,j,k1,k2,col[3],lf,block, row,Ii,J,n1,inc;
14: PetscReal norm1,norm2,rnorm,tol=10*PETSC_SMALL;
15: PetscScalar neg_one=-1.0,four=4.0,value[3];
16: IS perm, iscol;
17: PetscRandom rdm;
18: PetscBool doIcc=PETSC_TRUE,equal;
19: MatInfo minfo1,minfo2;
20: MatFactorInfo factinfo;
21: MatType type;
23: PetscInitialize(&argc,&args,(char*)0,help);if (ierr) return ierr;
24: MPI_Comm_size(PETSC_COMM_WORLD,&size);
25: if (size != 1) SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_SUP,"This is a uniprocessor example only!");
26: PetscOptionsGetInt(NULL,NULL,"-bs",&bs,NULL);
27: PetscOptionsGetInt(NULL,NULL,"-mbs",&mbs,NULL);
29: n = mbs*bs;
30: MatCreate(PETSC_COMM_SELF,&A);
31: MatSetSizes(A,n,n,PETSC_DETERMINE,PETSC_DETERMINE);
32: MatSetType(A,MATSEQBAIJ);
33: MatSetFromOptions(A);
34: MatSeqBAIJSetPreallocation(A,bs,nz,NULL);
36: MatCreate(PETSC_COMM_SELF,&sA);
37: MatSetSizes(sA,n,n,PETSC_DETERMINE,PETSC_DETERMINE);
38: MatSetType(sA,MATSEQSBAIJ);
39: MatSetFromOptions(sA);
40: MatGetType(sA,&type);
41: PetscObjectTypeCompare((PetscObject)sA,MATSEQSBAIJ,&doIcc);
42: MatSeqSBAIJSetPreallocation(sA,bs,nz,NULL);
43: MatSetOption(sA,MAT_IGNORE_LOWER_TRIANGULAR,PETSC_TRUE);
45: /* Test MatGetOwnershipRange() */
46: MatGetOwnershipRange(A,&Ii,&J);
47: MatGetOwnershipRange(sA,&i,&j);
48: if (i-Ii || j-J) {
49: PetscPrintf(PETSC_COMM_SELF,"Error: MatGetOwnershipRange() in MatSBAIJ format\n");
50: }
52: /* Assemble matrix */
53: if (bs == 1) {
54: PetscOptionsGetInt(NULL,NULL,"-test_problem",&prob,NULL);
55: if (prob == 1) { /* tridiagonal matrix */
56: value[0] = -1.0; value[1] = 2.0; value[2] = -1.0;
57: for (i=1; i<n-1; i++) {
58: col[0] = i-1; col[1] = i; col[2] = i+1;
59: MatSetValues(A,1,&i,3,col,value,INSERT_VALUES);
60: MatSetValues(sA,1,&i,3,col,value,INSERT_VALUES);
61: }
62: i = n - 1; col[0]=0; col[1] = n - 2; col[2] = n - 1;
64: value[0]= 0.1; value[1]=-1; value[2]=2;
66: MatSetValues(A,1,&i,3,col,value,INSERT_VALUES);
67: MatSetValues(sA,1,&i,3,col,value,INSERT_VALUES);
69: i = 0;
70: col[0] = n-1; col[1] = 1; col[2] = 0;
71: value[0] = 0.1; value[1] = -1.0; value[2] = 2;
73: MatSetValues(A,1,&i,3,col,value,INSERT_VALUES);
74: MatSetValues(sA,1,&i,3,col,value,INSERT_VALUES);
76: } else if (prob == 2) { /* matrix for the five point stencil */
77: n1 = (PetscInt) (PetscSqrtReal((PetscReal)n) + 0.001);
78: if (n1*n1 - n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"sqrt(n) must be a positive integer!");
79: for (i=0; i<n1; i++) {
80: for (j=0; j<n1; j++) {
81: Ii = j + n1*i;
82: if (i>0) {
83: J = Ii - n1;
84: MatSetValues(A,1,&Ii,1,&J,&neg_one,INSERT_VALUES);
85: MatSetValues(sA,1,&Ii,1,&J,&neg_one,INSERT_VALUES);
86: }
87: if (i<n1-1) {
88: J = Ii + n1;
89: MatSetValues(A,1,&Ii,1,&J,&neg_one,INSERT_VALUES);
90: MatSetValues(sA,1,&Ii,1,&J,&neg_one,INSERT_VALUES);
91: }
92: if (j>0) {
93: J = Ii - 1;
94: MatSetValues(A,1,&Ii,1,&J,&neg_one,INSERT_VALUES);
95: MatSetValues(sA,1,&Ii,1,&J,&neg_one,INSERT_VALUES);
96: }
97: if (j<n1-1) {
98: J = Ii + 1;
99: MatSetValues(A,1,&Ii,1,&J,&neg_one,INSERT_VALUES);
100: MatSetValues(sA,1,&Ii,1,&J,&neg_one,INSERT_VALUES);
101: }
102: MatSetValues(A,1,&Ii,1,&Ii,&four,INSERT_VALUES);
103: MatSetValues(sA,1,&Ii,1,&Ii,&four,INSERT_VALUES);
104: }
105: }
106: }
108: } else { /* bs > 1 */
109: for (block=0; block<n/bs; block++) {
110: /* diagonal blocks */
111: value[0] = -1.0; value[1] = 4.0; value[2] = -1.0;
112: for (i=1+block*bs; i<bs-1+block*bs; i++) {
113: col[0] = i-1; col[1] = i; col[2] = i+1;
114: MatSetValues(A,1,&i,3,col,value,INSERT_VALUES);
115: MatSetValues(sA,1,&i,3,col,value,INSERT_VALUES);
116: }
117: i = bs - 1+block*bs; col[0] = bs - 2+block*bs; col[1] = bs - 1+block*bs;
119: value[0]=-1.0; value[1]=4.0;
121: MatSetValues(A,1,&i,2,col,value,INSERT_VALUES);
122: MatSetValues(sA,1,&i,2,col,value,INSERT_VALUES);
124: i = 0+block*bs; col[0] = 0+block*bs; col[1] = 1+block*bs;
126: value[0]=4.0; value[1] = -1.0;
128: MatSetValues(A,1,&i,2,col,value,INSERT_VALUES);
129: MatSetValues(sA,1,&i,2,col,value,INSERT_VALUES);
130: }
131: /* off-diagonal blocks */
132: value[0]=-1.0;
133: for (i=0; i<(n/bs-1)*bs; i++) {
134: col[0]=i+bs;
136: MatSetValues(A,1,&i,1,col,value,INSERT_VALUES);
137: MatSetValues(sA,1,&i,1,col,value,INSERT_VALUES);
139: col[0]=i; row=i+bs;
141: MatSetValues(A,1,&row,1,col,value,INSERT_VALUES);
142: MatSetValues(sA,1,&row,1,col,value,INSERT_VALUES);
143: }
144: }
145: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
146: MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
148: MatAssemblyBegin(sA,MAT_FINAL_ASSEMBLY);
149: MatAssemblyEnd(sA,MAT_FINAL_ASSEMBLY);
151: /* Test MatGetInfo() of A and sA */
152: MatGetInfo(A,MAT_LOCAL,&minfo1);
153: MatGetInfo(sA,MAT_LOCAL,&minfo2);
154: i = (int) (minfo1.nz_used - minfo2.nz_used);
155: j = (int) (minfo1.nz_allocated - minfo2.nz_allocated);
156: k1 = (int) (minfo1.nz_allocated - minfo1.nz_used);
157: k2 = (int) (minfo2.nz_allocated - minfo2.nz_used);
158: if (i < 0 || j < 0 || k1 < 0 || k2 < 0) {
159: PetscPrintf(PETSC_COMM_SELF,"Error (compare A and sA): MatGetInfo()\n");
160: }
162: /* Test MatDuplicate() */
163: MatNorm(A,NORM_FROBENIUS,&norm1);
164: MatDuplicate(sA,MAT_COPY_VALUES,&sB);
165: MatEqual(sA,sB,&equal);
166: if (!equal) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_NOTSAMETYPE,"Error in MatDuplicate()");
168: /* Test MatNorm() */
169: MatNorm(A,NORM_FROBENIUS,&norm1);
170: MatNorm(sB,NORM_FROBENIUS,&norm2);
171: rnorm = PetscAbsReal(norm1-norm2)/norm2;
172: if (rnorm > tol) {
173: PetscPrintf(PETSC_COMM_SELF,"Error: MatNorm_FROBENIUS, NormA=%16.14e NormsB=%16.14e\n",norm1,norm2);
174: }
175: MatNorm(A,NORM_INFINITY,&norm1);
176: MatNorm(sB,NORM_INFINITY,&norm2);
177: rnorm = PetscAbsReal(norm1-norm2)/norm2;
178: if (rnorm > tol) {
179: PetscPrintf(PETSC_COMM_SELF,"Error: MatNorm_INFINITY(), NormA=%16.14e NormsB=%16.14e\n",norm1,norm2);
180: }
181: MatNorm(A,NORM_1,&norm1);
182: MatNorm(sB,NORM_1,&norm2);
183: rnorm = PetscAbsReal(norm1-norm2)/norm2;
184: if (rnorm > tol) {
185: PetscPrintf(PETSC_COMM_SELF,"Error: MatNorm_INFINITY(), NormA=%16.14e NormsB=%16.14e\n",norm1,norm2);
186: }
188: /* Test MatGetInfo(), MatGetSize(), MatGetBlockSize() */
189: MatGetInfo(A,MAT_LOCAL,&minfo1);
190: MatGetInfo(sB,MAT_LOCAL,&minfo2);
191: i = (int) (minfo1.nz_used - minfo2.nz_used);
192: j = (int) (minfo1.nz_allocated - minfo2.nz_allocated);
193: k1 = (int) (minfo1.nz_allocated - minfo1.nz_used);
194: k2 = (int) (minfo2.nz_allocated - minfo2.nz_used);
195: if (i < 0 || j < 0 || k1 < 0 || k2 < 0) {
196: PetscPrintf(PETSC_COMM_SELF,"Error(compare A and sB): MatGetInfo()\n");
197: }
199: MatGetSize(A,&Ii,&J);
200: MatGetSize(sB,&i,&j);
201: if (i-Ii || j-J) {
202: PetscPrintf(PETSC_COMM_SELF,"Error: MatGetSize()\n");
203: }
205: MatGetBlockSize(A, &Ii);
206: MatGetBlockSize(sB, &i);
207: if (i-Ii) {
208: PetscPrintf(PETSC_COMM_SELF,"Error: MatGetBlockSize()\n");
209: }
211: PetscRandomCreate(PETSC_COMM_SELF,&rdm);
212: PetscRandomSetFromOptions(rdm);
213: VecCreateSeq(PETSC_COMM_SELF,n,&x);
214: VecDuplicate(x,&s1);
215: VecDuplicate(x,&s2);
216: VecDuplicate(x,&y);
217: VecDuplicate(x,&b);
218: VecSetRandom(x,rdm);
220: /* Test MatDiagonalScale(), MatGetDiagonal(), MatScale() */
221: #if !defined(PETSC_USE_COMPLEX)
222: /* Scaling matrix with complex numbers results non-spd matrix,
223: causing crash of MatForwardSolve() and MatBackwardSolve() */
224: MatDiagonalScale(A,x,x);
225: MatDiagonalScale(sB,x,x);
226: MatMultEqual(A,sB,10,&equal);
227: if (!equal) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_NOTSAMETYPE,"Error in MatDiagonalScale");
229: MatGetDiagonal(A,s1);
230: MatGetDiagonal(sB,s2);
231: VecAXPY(s2,neg_one,s1);
232: VecNorm(s2,NORM_1,&norm1);
233: if (norm1>tol) {
234: PetscPrintf(PETSC_COMM_SELF,"Error:MatGetDiagonal(), ||s1-s2||=%g\n",(double)norm1);
235: }
237: {
238: PetscScalar alpha=0.1;
239: MatScale(A,alpha);
240: MatScale(sB,alpha);
241: }
242: #endif
244: /* Test MatGetRowMaxAbs() */
245: MatGetRowMaxAbs(A,s1,NULL);
246: MatGetRowMaxAbs(sB,s2,NULL);
247: VecNorm(s1,NORM_1,&norm1);
248: VecNorm(s2,NORM_1,&norm2);
249: norm1 -= norm2;
250: if (norm1<-tol || norm1>tol) {
251: PetscPrintf(PETSC_COMM_SELF,"Error:MatGetRowMaxAbs() \n");
252: }
254: /* Test MatMult() */
255: for (i=0; i<40; i++) {
256: VecSetRandom(x,rdm);
257: MatMult(A,x,s1);
258: MatMult(sB,x,s2);
259: VecNorm(s1,NORM_1,&norm1);
260: VecNorm(s2,NORM_1,&norm2);
261: norm1 -= norm2;
262: if (norm1<-tol || norm1>tol) {
263: PetscPrintf(PETSC_COMM_SELF,"Error: MatMult(), norm1-norm2: %g\n",(double)norm1);
264: }
265: }
267: /* MatMultAdd() */
268: for (i=0; i<40; i++) {
269: VecSetRandom(x,rdm);
270: VecSetRandom(y,rdm);
271: MatMultAdd(A,x,y,s1);
272: MatMultAdd(sB,x,y,s2);
273: VecNorm(s1,NORM_1,&norm1);
274: VecNorm(s2,NORM_1,&norm2);
275: norm1 -= norm2;
276: if (norm1<-tol || norm1>tol) {
277: PetscPrintf(PETSC_COMM_SELF,"Error:MatMultAdd(), norm1-norm2: %g\n",(double)norm1);
278: }
279: }
281: /* Test MatMatMult() for sbaij and dense matrices */
282: MatCreateSeqDense(PETSC_COMM_SELF,n,5*n,NULL,&B);
283: MatSetRandom(B,rdm);
284: MatMatMult(sA,B,MAT_INITIAL_MATRIX,PETSC_DEFAULT,&C);
285: MatMatMultEqual(sA,B,C,5*n,&equal);
286: if (!equal) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Error: MatMatMult()");
287: MatDestroy(&C);
288: MatDestroy(&B);
290: /* Test MatCholeskyFactor(), MatICCFactor() with natural ordering */
291: MatGetOrdering(A,MATORDERINGNATURAL,&perm,&iscol);
292: ISDestroy(&iscol);
293: norm1 = tol;
294: inc = bs;
296: /* initialize factinfo */
297: PetscMemzero(&factinfo,sizeof(MatFactorInfo));
299: for (lf=-1; lf<10; lf += inc) {
300: if (lf==-1) { /* Cholesky factor of sB (duplicate sA) */
301: factinfo.fill = 5.0;
303: MatGetFactor(sB,MATSOLVERPETSC,MAT_FACTOR_CHOLESKY,&sFactor);
304: MatCholeskyFactorSymbolic(sFactor,sB,perm,&factinfo);
305: } else if (!doIcc) break;
306: else { /* incomplete Cholesky factor */
307: factinfo.fill = 5.0;
308: factinfo.levels = lf;
310: MatGetFactor(sB,MATSOLVERPETSC,MAT_FACTOR_ICC,&sFactor);
311: MatICCFactorSymbolic(sFactor,sB,perm,&factinfo);
312: }
313: MatCholeskyFactorNumeric(sFactor,sB,&factinfo);
314: /* MatView(sFactor, PETSC_VIEWER_DRAW_WORLD); */
316: /* test MatGetDiagonal on numeric factor */
317: /*
318: if (lf == -1) {
319: MatGetDiagonal(sFactor,s1);
320: printf(" in ex74.c, diag: \n");
321: VecView(s1,PETSC_VIEWER_STDOUT_SELF);
322: }
323: */
325: MatMult(sB,x,b);
327: /* test MatForwardSolve() and MatBackwardSolve() */
328: if (lf == -1) {
329: MatForwardSolve(sFactor,b,s1);
330: MatBackwardSolve(sFactor,s1,s2);
331: VecAXPY(s2,neg_one,x);
332: VecNorm(s2,NORM_2,&norm2);
333: if (10*norm1 < norm2) {
334: PetscPrintf(PETSC_COMM_SELF,"MatForwardSolve and BackwardSolve: Norm of error=%g, bs=%D\n",(double)norm2,bs);
335: }
336: }
338: /* test MatSolve() */
339: MatSolve(sFactor,b,y);
340: MatDestroy(&sFactor);
341: /* Check the error */
342: VecAXPY(y,neg_one,x);
343: VecNorm(y,NORM_2,&norm2);
344: if (10*norm1 < norm2 && lf-inc != -1) {
345: PetscPrintf(PETSC_COMM_SELF,"lf=%D, %D, Norm of error=%g, %g\n",lf-inc,lf,(double)norm1,(double)norm2);
346: }
347: norm1 = norm2;
348: if (norm2 < tol && lf != -1) break;
349: }
351: #if defined(PETSC_HAVE_MUMPS)
352: MatGetFactor(sA,MATSOLVERMUMPS,MAT_FACTOR_CHOLESKY,&sFactor);
353: MatCholeskyFactorSymbolic(sFactor,sA,NULL,NULL);
354: MatCholeskyFactorNumeric(sFactor,sA,NULL);
355: for (i=0; i<10; i++) {
356: VecSetRandom(b,rdm);
357: MatSolve(sFactor,b,y);
358: /* Check the error */
359: MatMult(sA,y,x);
360: VecAXPY(x,neg_one,b);
361: VecNorm(x,NORM_2,&norm2);
362: if (norm2>tol) {
363: PetscPrintf(PETSC_COMM_SELF,"Error:MatSolve(), norm2: %g\n",(double)norm2);
364: }
365: }
366: MatDestroy(&sFactor);
367: #endif
369: ISDestroy(&perm);
371: MatDestroy(&A);
372: MatDestroy(&sB);
373: MatDestroy(&sA);
374: VecDestroy(&x);
375: VecDestroy(&y);
376: VecDestroy(&s1);
377: VecDestroy(&s2);
378: VecDestroy(&b);
379: PetscRandomDestroy(&rdm);
381: PetscFinalize();
382: return ierr;
383: }
386: /*TEST
388: test:
389: args: -bs {{1 2 3 4 5 6 7 8}}
391: TEST*/