Actual source code: ex74.c
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;
9: Vec x,y,b,s1,s2;
10: Mat A; /* linear system matrix */
11: Mat sA,sB,sFactor,B,C; /* symmetric matrices */
12: PetscInt n,mbs=16,bs=1,nz=3,prob=1,i,j,k1,k2,col[3],lf,block, row,Ii,J,n1,inc;
13: PetscReal norm1,norm2,rnorm,tol=10*PETSC_SMALL;
14: PetscScalar neg_one=-1.0,four=4.0,value[3];
15: IS perm, iscol;
16: PetscRandom rdm;
17: PetscBool doIcc=PETSC_TRUE,equal;
18: MatInfo minfo1,minfo2;
19: MatFactorInfo factinfo;
20: MatType type;
22: PetscInitialize(&argc,&args,(char*)0,help);
23: MPI_Comm_size(PETSC_COMM_WORLD,&size);
25: PetscOptionsGetInt(NULL,NULL,"-bs",&bs,NULL);
26: PetscOptionsGetInt(NULL,NULL,"-mbs",&mbs,NULL);
28: n = mbs*bs;
29: MatCreate(PETSC_COMM_SELF,&A);
30: MatSetSizes(A,n,n,PETSC_DETERMINE,PETSC_DETERMINE);
31: MatSetType(A,MATSEQBAIJ);
32: MatSetFromOptions(A);
33: MatSeqBAIJSetPreallocation(A,bs,nz,NULL);
35: MatCreate(PETSC_COMM_SELF,&sA);
36: MatSetSizes(sA,n,n,PETSC_DETERMINE,PETSC_DETERMINE);
37: MatSetType(sA,MATSEQSBAIJ);
38: MatSetFromOptions(sA);
39: MatGetType(sA,&type);
40: PetscObjectTypeCompare((PetscObject)sA,MATSEQSBAIJ,&doIcc);
41: MatSeqSBAIJSetPreallocation(sA,bs,nz,NULL);
42: MatSetOption(sA,MAT_IGNORE_LOWER_TRIANGULAR,PETSC_TRUE);
44: /* Test MatGetOwnershipRange() */
45: MatGetOwnershipRange(A,&Ii,&J);
46: MatGetOwnershipRange(sA,&i,&j);
47: if (i-Ii || j-J) {
48: PetscPrintf(PETSC_COMM_SELF,"Error: MatGetOwnershipRange() in MatSBAIJ format\n");
49: }
51: /* Assemble matrix */
52: if (bs == 1) {
53: PetscOptionsGetInt(NULL,NULL,"-test_problem",&prob,NULL);
54: if (prob == 1) { /* tridiagonal matrix */
55: value[0] = -1.0; value[1] = 2.0; value[2] = -1.0;
56: for (i=1; i<n-1; i++) {
57: col[0] = i-1; col[1] = i; col[2] = i+1;
58: MatSetValues(A,1,&i,3,col,value,INSERT_VALUES);
59: MatSetValues(sA,1,&i,3,col,value,INSERT_VALUES);
60: }
61: i = n - 1; col[0]=0; col[1] = n - 2; col[2] = n - 1;
63: value[0]= 0.1; value[1]=-1; value[2]=2;
65: MatSetValues(A,1,&i,3,col,value,INSERT_VALUES);
66: MatSetValues(sA,1,&i,3,col,value,INSERT_VALUES);
68: i = 0;
69: col[0] = n-1; col[1] = 1; col[2] = 0;
70: value[0] = 0.1; value[1] = -1.0; value[2] = 2;
72: MatSetValues(A,1,&i,3,col,value,INSERT_VALUES);
73: MatSetValues(sA,1,&i,3,col,value,INSERT_VALUES);
75: } else if (prob == 2) { /* matrix for the five point stencil */
76: n1 = (PetscInt) (PetscSqrtReal((PetscReal)n) + 0.001);
78: for (i=0; i<n1; i++) {
79: for (j=0; j<n1; j++) {
80: Ii = j + n1*i;
81: if (i>0) {
82: J = Ii - n1;
83: MatSetValues(A,1,&Ii,1,&J,&neg_one,INSERT_VALUES);
84: MatSetValues(sA,1,&Ii,1,&J,&neg_one,INSERT_VALUES);
85: }
86: if (i<n1-1) {
87: J = Ii + n1;
88: MatSetValues(A,1,&Ii,1,&J,&neg_one,INSERT_VALUES);
89: MatSetValues(sA,1,&Ii,1,&J,&neg_one,INSERT_VALUES);
90: }
91: if (j>0) {
92: J = Ii - 1;
93: MatSetValues(A,1,&Ii,1,&J,&neg_one,INSERT_VALUES);
94: MatSetValues(sA,1,&Ii,1,&J,&neg_one,INSERT_VALUES);
95: }
96: if (j<n1-1) {
97: J = Ii + 1;
98: MatSetValues(A,1,&Ii,1,&J,&neg_one,INSERT_VALUES);
99: MatSetValues(sA,1,&Ii,1,&J,&neg_one,INSERT_VALUES);
100: }
101: MatSetValues(A,1,&Ii,1,&Ii,&four,INSERT_VALUES);
102: MatSetValues(sA,1,&Ii,1,&Ii,&four,INSERT_VALUES);
103: }
104: }
105: }
107: } else { /* bs > 1 */
108: for (block=0; block<n/bs; block++) {
109: /* diagonal blocks */
110: value[0] = -1.0; value[1] = 4.0; value[2] = -1.0;
111: for (i=1+block*bs; i<bs-1+block*bs; i++) {
112: col[0] = i-1; col[1] = i; col[2] = i+1;
113: MatSetValues(A,1,&i,3,col,value,INSERT_VALUES);
114: MatSetValues(sA,1,&i,3,col,value,INSERT_VALUES);
115: }
116: i = bs - 1+block*bs; col[0] = bs - 2+block*bs; col[1] = bs - 1+block*bs;
118: value[0]=-1.0; value[1]=4.0;
120: MatSetValues(A,1,&i,2,col,value,INSERT_VALUES);
121: MatSetValues(sA,1,&i,2,col,value,INSERT_VALUES);
123: i = 0+block*bs; col[0] = 0+block*bs; col[1] = 1+block*bs;
125: value[0]=4.0; value[1] = -1.0;
127: MatSetValues(A,1,&i,2,col,value,INSERT_VALUES);
128: MatSetValues(sA,1,&i,2,col,value,INSERT_VALUES);
129: }
130: /* off-diagonal blocks */
131: value[0]=-1.0;
132: for (i=0; i<(n/bs-1)*bs; i++) {
133: col[0]=i+bs;
135: MatSetValues(A,1,&i,1,col,value,INSERT_VALUES);
136: MatSetValues(sA,1,&i,1,col,value,INSERT_VALUES);
138: col[0]=i; row=i+bs;
140: MatSetValues(A,1,&row,1,col,value,INSERT_VALUES);
141: MatSetValues(sA,1,&row,1,col,value,INSERT_VALUES);
142: }
143: }
144: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
145: MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
147: MatAssemblyBegin(sA,MAT_FINAL_ASSEMBLY);
148: MatAssemblyEnd(sA,MAT_FINAL_ASSEMBLY);
150: /* Test MatGetInfo() of A and sA */
151: MatGetInfo(A,MAT_LOCAL,&minfo1);
152: MatGetInfo(sA,MAT_LOCAL,&minfo2);
153: i = (int) (minfo1.nz_used - minfo2.nz_used);
154: j = (int) (minfo1.nz_allocated - minfo2.nz_allocated);
155: k1 = (int) (minfo1.nz_allocated - minfo1.nz_used);
156: k2 = (int) (minfo2.nz_allocated - minfo2.nz_used);
157: if (i < 0 || j < 0 || k1 < 0 || k2 < 0) {
158: PetscPrintf(PETSC_COMM_SELF,"Error (compare A and sA): MatGetInfo()\n");
159: }
161: /* Test MatDuplicate() */
162: MatNorm(A,NORM_FROBENIUS,&norm1);
163: MatDuplicate(sA,MAT_COPY_VALUES,&sB);
164: MatEqual(sA,sB,&equal);
167: /* Test MatNorm() */
168: MatNorm(A,NORM_FROBENIUS,&norm1);
169: MatNorm(sB,NORM_FROBENIUS,&norm2);
170: rnorm = PetscAbsReal(norm1-norm2)/norm2;
171: if (rnorm > tol) {
172: PetscPrintf(PETSC_COMM_SELF,"Error: MatNorm_FROBENIUS, NormA=%16.14e NormsB=%16.14e\n",(double)norm1,(double)norm2);
173: }
174: MatNorm(A,NORM_INFINITY,&norm1);
175: MatNorm(sB,NORM_INFINITY,&norm2);
176: rnorm = PetscAbsReal(norm1-norm2)/norm2;
177: if (rnorm > tol) {
178: PetscPrintf(PETSC_COMM_SELF,"Error: MatNorm_INFINITY(), NormA=%16.14e NormsB=%16.14e\n",(double)norm1,(double)norm2);
179: }
180: MatNorm(A,NORM_1,&norm1);
181: MatNorm(sB,NORM_1,&norm2);
182: rnorm = PetscAbsReal(norm1-norm2)/norm2;
183: if (rnorm > tol) {
184: PetscPrintf(PETSC_COMM_SELF,"Error: MatNorm_INFINITY(), NormA=%16.14e NormsB=%16.14e\n",(double)norm1,(double)norm2);
185: }
187: /* Test MatGetInfo(), MatGetSize(), MatGetBlockSize() */
188: MatGetInfo(A,MAT_LOCAL,&minfo1);
189: MatGetInfo(sB,MAT_LOCAL,&minfo2);
190: i = (int) (minfo1.nz_used - minfo2.nz_used);
191: j = (int) (minfo1.nz_allocated - minfo2.nz_allocated);
192: k1 = (int) (minfo1.nz_allocated - minfo1.nz_used);
193: k2 = (int) (minfo2.nz_allocated - minfo2.nz_used);
194: if (i < 0 || j < 0 || k1 < 0 || k2 < 0) {
195: PetscPrintf(PETSC_COMM_SELF,"Error(compare A and sB): MatGetInfo()\n");
196: }
198: MatGetSize(A,&Ii,&J);
199: MatGetSize(sB,&i,&j);
200: if (i-Ii || j-J) {
201: PetscPrintf(PETSC_COMM_SELF,"Error: MatGetSize()\n");
202: }
204: MatGetBlockSize(A, &Ii);
205: MatGetBlockSize(sB, &i);
206: if (i-Ii) {
207: PetscPrintf(PETSC_COMM_SELF,"Error: MatGetBlockSize()\n");
208: }
210: PetscRandomCreate(PETSC_COMM_SELF,&rdm);
211: PetscRandomSetFromOptions(rdm);
212: VecCreateSeq(PETSC_COMM_SELF,n,&x);
213: VecDuplicate(x,&s1);
214: VecDuplicate(x,&s2);
215: VecDuplicate(x,&y);
216: VecDuplicate(x,&b);
217: VecSetRandom(x,rdm);
219: /* Test MatDiagonalScale(), MatGetDiagonal(), MatScale() */
220: #if !defined(PETSC_USE_COMPLEX)
221: /* Scaling matrix with complex numbers results non-spd matrix,
222: causing crash of MatForwardSolve() and MatBackwardSolve() */
223: MatDiagonalScale(A,x,x);
224: MatDiagonalScale(sB,x,x);
225: MatMultEqual(A,sB,10,&equal);
228: MatGetDiagonal(A,s1);
229: MatGetDiagonal(sB,s2);
230: VecAXPY(s2,neg_one,s1);
231: VecNorm(s2,NORM_1,&norm1);
232: if (norm1>tol) {
233: PetscPrintf(PETSC_COMM_SELF,"Error:MatGetDiagonal(), ||s1-s2||=%g\n",(double)norm1);
234: }
236: {
237: PetscScalar alpha=0.1;
238: MatScale(A,alpha);
239: MatScale(sB,alpha);
240: }
241: #endif
243: /* Test MatGetRowMaxAbs() */
244: MatGetRowMaxAbs(A,s1,NULL);
245: MatGetRowMaxAbs(sB,s2,NULL);
246: VecNorm(s1,NORM_1,&norm1);
247: VecNorm(s2,NORM_1,&norm2);
248: norm1 -= norm2;
249: if (norm1<-tol || norm1>tol) {
250: PetscPrintf(PETSC_COMM_SELF,"Error:MatGetRowMaxAbs() \n");
251: }
253: /* Test MatMult() */
254: for (i=0; i<40; i++) {
255: VecSetRandom(x,rdm);
256: MatMult(A,x,s1);
257: MatMult(sB,x,s2);
258: VecNorm(s1,NORM_1,&norm1);
259: VecNorm(s2,NORM_1,&norm2);
260: norm1 -= norm2;
261: if (norm1<-tol || norm1>tol) {
262: PetscPrintf(PETSC_COMM_SELF,"Error: MatMult(), norm1-norm2: %g\n",(double)norm1);
263: }
264: }
266: /* MatMultAdd() */
267: for (i=0; i<40; i++) {
268: VecSetRandom(x,rdm);
269: VecSetRandom(y,rdm);
270: MatMultAdd(A,x,y,s1);
271: MatMultAdd(sB,x,y,s2);
272: VecNorm(s1,NORM_1,&norm1);
273: VecNorm(s2,NORM_1,&norm2);
274: norm1 -= norm2;
275: if (norm1<-tol || norm1>tol) {
276: PetscPrintf(PETSC_COMM_SELF,"Error:MatMultAdd(), norm1-norm2: %g\n",(double)norm1);
277: }
278: }
280: /* Test MatMatMult() for sbaij and dense matrices */
281: MatCreateSeqDense(PETSC_COMM_SELF,n,5*n,NULL,&B);
282: MatSetRandom(B,rdm);
283: MatMatMult(sA,B,MAT_INITIAL_MATRIX,PETSC_DEFAULT,&C);
284: MatMatMultEqual(sA,B,C,5*n,&equal);
286: MatDestroy(&C);
287: MatDestroy(&B);
289: /* Test MatCholeskyFactor(), MatICCFactor() with natural ordering */
290: MatGetOrdering(A,MATORDERINGNATURAL,&perm,&iscol);
291: ISDestroy(&iscol);
292: norm1 = tol;
293: inc = bs;
295: /* initialize factinfo */
296: PetscMemzero(&factinfo,sizeof(MatFactorInfo));
298: for (lf=-1; lf<10; lf += inc) {
299: if (lf==-1) { /* Cholesky factor of sB (duplicate sA) */
300: factinfo.fill = 5.0;
302: MatGetFactor(sB,MATSOLVERPETSC,MAT_FACTOR_CHOLESKY,&sFactor);
303: MatCholeskyFactorSymbolic(sFactor,sB,perm,&factinfo);
304: } else if (!doIcc) break;
305: else { /* incomplete Cholesky factor */
306: factinfo.fill = 5.0;
307: factinfo.levels = lf;
309: MatGetFactor(sB,MATSOLVERPETSC,MAT_FACTOR_ICC,&sFactor);
310: MatICCFactorSymbolic(sFactor,sB,perm,&factinfo);
311: }
312: MatCholeskyFactorNumeric(sFactor,sB,&factinfo);
313: /* MatView(sFactor, PETSC_VIEWER_DRAW_WORLD); */
315: /* test MatGetDiagonal on numeric factor */
316: /*
317: if (lf == -1) {
318: MatGetDiagonal(sFactor,s1);
319: printf(" in ex74.c, diag: \n");
320: VecView(s1,PETSC_VIEWER_STDOUT_SELF);
321: }
322: */
324: MatMult(sB,x,b);
326: /* test MatForwardSolve() and MatBackwardSolve() */
327: if (lf == -1) {
328: MatForwardSolve(sFactor,b,s1);
329: MatBackwardSolve(sFactor,s1,s2);
330: VecAXPY(s2,neg_one,x);
331: VecNorm(s2,NORM_2,&norm2);
332: if (10*norm1 < norm2) {
333: PetscPrintf(PETSC_COMM_SELF,"MatForwardSolve and BackwardSolve: Norm of error=%g, bs=%" PetscInt_FMT "\n",(double)norm2,bs);
334: }
335: }
337: /* test MatSolve() */
338: MatSolve(sFactor,b,y);
339: MatDestroy(&sFactor);
340: /* Check the error */
341: VecAXPY(y,neg_one,x);
342: VecNorm(y,NORM_2,&norm2);
343: if (10*norm1 < norm2 && lf-inc != -1) {
344: PetscPrintf(PETSC_COMM_SELF,"lf=%" PetscInt_FMT ", %" PetscInt_FMT ", Norm of error=%g, %g\n",lf-inc,lf,(double)norm1,(double)norm2);
345: }
346: norm1 = norm2;
347: if (norm2 < tol && lf != -1) break;
348: }
350: #if defined(PETSC_HAVE_MUMPS)
351: MatGetFactor(sA,MATSOLVERMUMPS,MAT_FACTOR_CHOLESKY,&sFactor);
352: MatCholeskyFactorSymbolic(sFactor,sA,NULL,NULL);
353: MatCholeskyFactorNumeric(sFactor,sA,NULL);
354: for (i=0; i<10; i++) {
355: VecSetRandom(b,rdm);
356: MatSolve(sFactor,b,y);
357: /* Check the error */
358: MatMult(sA,y,x);
359: VecAXPY(x,neg_one,b);
360: VecNorm(x,NORM_2,&norm2);
361: if (norm2>tol) {
362: PetscPrintf(PETSC_COMM_SELF,"Error:MatSolve(), norm2: %g\n",(double)norm2);
363: }
364: }
365: MatDestroy(&sFactor);
366: #endif
368: ISDestroy(&perm);
370: MatDestroy(&A);
371: MatDestroy(&sB);
372: MatDestroy(&sA);
373: VecDestroy(&x);
374: VecDestroy(&y);
375: VecDestroy(&s1);
376: VecDestroy(&s2);
377: VecDestroy(&b);
378: PetscRandomDestroy(&rdm);
380: PetscFinalize();
381: return 0;
382: }
384: /*TEST
386: test:
387: args: -bs {{1 2 3 4 5 6 7 8}}
389: TEST*/