Actual source code: ex74.c
petsc-3.3-p6 2013-02-11
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;
11: PetscErrorCode ierr;
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=1.e-10;
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: const 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(PETSC_NULL,"-bs",&bs,PETSC_NULL);
29: PetscOptionsGetInt(PETSC_NULL,"-mbs",&mbs,PETSC_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,PETSC_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,PETSC_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(PETSC_NULL,"-test_problem",&prob,PETSC_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;
65: 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;
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 = (int) (PetscSqrtReal((PetscReal)n) + 0.001);
77: if (n1*n1 - n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"sqrt(n) must be a positive interger!");
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;
117: value[0]=-1.0; value[1]=4.0;
118: MatSetValues(A,1,&i,2,col,value,INSERT_VALUES);
119: MatSetValues(sA,1,&i,2,col,value,INSERT_VALUES);
121: i = 0+block*bs; col[0] = 0+block*bs; col[1] = 1+block*bs;
122: value[0]=4.0; value[1] = -1.0;
123: MatSetValues(A,1,&i,2,col,value,INSERT_VALUES);
124: MatSetValues(sA,1,&i,2,col,value,INSERT_VALUES);
125: }
126: /* off-diagonal blocks */
127: value[0]=-1.0;
128: for (i=0; i<(n/bs-1)*bs; i++){
129: col[0]=i+bs;
130: MatSetValues(A,1,&i,1,col,value,INSERT_VALUES);
131: MatSetValues(sA,1,&i,1,col,value,INSERT_VALUES);
132: col[0]=i; row=i+bs;
133: MatSetValues(A,1,&row,1,col,value,INSERT_VALUES);
134: MatSetValues(sA,1,&row,1,col,value,INSERT_VALUES);
135: }
136: }
137: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
138: MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
139:
140: MatAssemblyBegin(sA,MAT_FINAL_ASSEMBLY);
141: MatAssemblyEnd(sA,MAT_FINAL_ASSEMBLY);
143: /* Test MatGetInfo() of A and sA */
144: MatGetInfo(A,MAT_LOCAL,&minfo1);
145: MatGetInfo(sA,MAT_LOCAL,&minfo2);
146: /*
147: printf("A matrix nonzeros (BAIJ format) = %d, allocated nonzeros= %d\n", (int)minfo1.nz_used,(int)minfo1.nz_allocated);
148: printf("sA matrix nonzeros(SBAIJ format) = %d, allocated nonzeros= %d\n", (int)minfo2.nz_used,(int)minfo2.nz_allocated);
149: */
150: i = (int) (minfo1.nz_used - minfo2.nz_used);
151: j = (int) (minfo1.nz_allocated - minfo2.nz_allocated);
152: k1 = (int) (minfo1.nz_allocated - minfo1.nz_used);
153: k2 = (int) (minfo2.nz_allocated - minfo2.nz_used);
154: if (i < 0 || j < 0 || k1 < 0 || k2 < 0) {
155: PetscPrintf(PETSC_COMM_SELF,"Error (compare A and sA): MatGetInfo()\n");
156: }
158: /* Test MatDuplicate() */
159: MatNorm(A,NORM_FROBENIUS,&norm1);
160: MatDuplicate(sA,MAT_COPY_VALUES,&sB);
161: MatEqual(sA,sB,&equal);
162: if (!equal) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_NOTSAMETYPE,"Error in MatDuplicate()");
164: /* Test MatNorm() */
165: MatNorm(A,NORM_FROBENIUS,&norm1);
166: MatNorm(sB,NORM_FROBENIUS,&norm2);
167: rnorm = PetscAbsScalar(norm1-norm2)/norm2;
168: if (rnorm > tol){
169: PetscPrintf(PETSC_COMM_SELF,"Error: MatNorm_FROBENIUS, NormA=%16.14e NormsB=%16.14e\n",norm1,norm2);
170: }
171: MatNorm(A,NORM_INFINITY,&norm1);
172: MatNorm(sB,NORM_INFINITY,&norm2);
173: rnorm = PetscAbsScalar(norm1-norm2)/norm2;
174: if (rnorm > tol){
175: PetscPrintf(PETSC_COMM_SELF,"Error: MatNorm_INFINITY(), NormA=%16.14e NormsB=%16.14e\n",norm1,norm2);
176: }
177: MatNorm(A,NORM_1,&norm1);
178: MatNorm(sB,NORM_1,&norm2);
179: rnorm = PetscAbsScalar(norm1-norm2)/norm2;
180: if (rnorm > tol){
181: PetscPrintf(PETSC_COMM_SELF,"Error: MatNorm_INFINITY(), NormA=%16.14e NormsB=%16.14e\n",norm1,norm2);
182: }
184: /* Test MatGetInfo(), MatGetSize(), MatGetBlockSize() */
185: MatGetInfo(A,MAT_LOCAL,&minfo1);
186: MatGetInfo(sB,MAT_LOCAL,&minfo2);
187: /*
188: printf("matrix nonzeros (BAIJ format) = %d, allocated nonzeros= %d\n", (int)minfo1.nz_used,(int)minfo1.nz_allocated);
189: printf("matrix nonzeros(SBAIJ format) = %d, allocated nonzeros= %d\n", (int)minfo2.nz_used,(int)minfo2.nz_allocated);
190: */
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: }
204:
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",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,PETSC_NULL);
246: MatGetRowMaxAbs(sB,s2,PETSC_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",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",norm1);
278: }
279: }
281: /* Test MatCholeskyFactor(), MatICCFactor() with natural ordering */
282: MatGetOrdering(A,MATORDERINGNATURAL,&perm,&iscol);
283: ISDestroy(&iscol);
284: norm1 = tol;
285: inc = bs;
287: /* initialize factinfo */
288: PetscMemzero(&factinfo,sizeof(MatFactorInfo));
290: for (lf=-1; lf<10; lf += inc){
291: if (lf==-1) { /* Cholesky factor of sB (duplicate sA) */
292: factinfo.fill = 5.0;
293: MatGetFactor(sB,MATSOLVERPETSC,MAT_FACTOR_CHOLESKY,&sC);
294: MatCholeskyFactorSymbolic(sC,sB,perm,&factinfo);
295: } else if (!doIcc){
296: break;
297: } else { /* incomplete Cholesky factor */
298: factinfo.fill = 5.0;
299: factinfo.levels = lf;
300: MatGetFactor(sB,MATSOLVERPETSC,MAT_FACTOR_ICC,&sC);
301: MatICCFactorSymbolic(sC,sB,perm,&factinfo);
302: }
303: MatCholeskyFactorNumeric(sC,sB,&factinfo);
304: /* MatView(sC, PETSC_VIEWER_DRAW_WORLD); */
306: /* test MatGetDiagonal on numeric factor */
307: /*
308: if (lf == -1) {
309: MatGetDiagonal(sC,s1);
310: printf(" in ex74.c, diag: \n");
311: VecView(s1,PETSC_VIEWER_STDOUT_SELF);
312: }
313: */
315: MatMult(sB,x,b);
317: /* test MatForwardSolve() and MatBackwardSolve() */
318: if (lf == -1){
319: MatForwardSolve(sC,b,s1);
320: MatBackwardSolve(sC,s1,s2);
321: VecAXPY(s2,neg_one,x);
322: VecNorm(s2,NORM_2,&norm2);
323: if (10*norm1 < norm2){
324: PetscPrintf(PETSC_COMM_SELF,"MatForwardSolve and BackwardSolve: Norm of error=%G, bs=%d\n",norm2,bs);
325: }
326: }
328: /* test MatSolve() */
329: MatSolve(sC,b,y);
330: MatDestroy(&sC);
331: /* Check the error */
332: VecAXPY(y,neg_one,x);
333: VecNorm(y,NORM_2,&norm2);
334: /* printf("lf: %d, error: %G\n", lf,norm2); */
335: if (10*norm1 < norm2 && lf-inc != -1){
336: PetscPrintf(PETSC_COMM_SELF,"lf=%D, %D, Norm of error=%G, %G\n",lf-inc,lf,norm1,norm2);
337: }
338: norm1 = norm2;
339: if (norm2 < tol && lf != -1) break;
340: }
342: ISDestroy(&perm);
344: MatDestroy(&A);
345: MatDestroy(&sB);
346: MatDestroy(&sA);
347: VecDestroy(&x);
348: VecDestroy(&y);
349: VecDestroy(&s1);
350: VecDestroy(&s2);
351: VecDestroy(&b);
352: PetscRandomDestroy(&rdm);
354: PetscFinalize();
355: return 0;
356: }