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