Actual source code: ex74.c
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,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],alpha=0.1;
18: IS perm, iscol;
19: PetscRandom rdm;
20: PetscTruth 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_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: ierr=MatCreateSeqBAIJ(PETSC_COMM_SELF,bs,n,n,nz,PETSC_NULL, &A);
33: MatCreate(PETSC_COMM_SELF,&sA);
34: MatSetSizes(sA,n,n,PETSC_DETERMINE,PETSC_DETERMINE);
35: MatSetType(sA,MATSEQSBAIJ);
36: /* -mat_type <seqsbaij_derived type>, e.g., seqsbaijspooles, sbaijmumps */
37: MatSetFromOptions(sA);
38: MatGetType(sA,&type);
39: PetscTypeCompare((PetscObject)sA,MATSEQSBAIJ,&doIcc);
40: MatSeqSBAIJSetPreallocation(sA,bs,nz,PETSC_NULL);
41: MatSetOption(sA,MAT_IGNORE_LOWER_TRIANGULAR);
43: /* Test MatGetOwnershipRange() */
44: MatGetOwnershipRange(A,&Ii,&J);
45: MatGetOwnershipRange(sA,&i,&j);
46: if (i-Ii || j-J){
47: PetscPrintf(PETSC_COMM_SELF,"Error: MatGetOwnershipRange() in MatSBAIJ format\n");
48: }
50: /* Assemble matrix */
51: if (bs == 1){
52: PetscOptionsGetInt(PETSC_NULL,"-test_problem",&prob,PETSC_NULL);
53: if (prob == 1){ /* tridiagonal matrix */
54: value[0] = -1.0; value[1] = 2.0; value[2] = -1.0;
55: for (i=1; i<n-1; i++) {
56: col[0] = i-1; col[1] = i; col[2] = i+1;
57: MatSetValues(A,1,&i,3,col,value,INSERT_VALUES);
58: MatSetValues(sA,1,&i,3,col,value,INSERT_VALUES);
59: }
60: i = n - 1; col[0]=0; col[1] = n - 2; col[2] = n - 1;
61: value[0]= 0.1; value[1]=-1; value[2]=2;
62: MatSetValues(A,1,&i,3,col,value,INSERT_VALUES);
63: MatSetValues(sA,1,&i,3,col,value,INSERT_VALUES);
65: i = 0;
66: col[0] = n-1; col[1] = 1; col[2]=0;
67: value[0] = 0.1; value[1] = -1.0; value[2]=2;
68: MatSetValues(A,1,&i,3,col,value,INSERT_VALUES);
69: MatSetValues(sA,1,&i,3,col,value,INSERT_VALUES);
70: }
71: else if (prob ==2){ /* matrix for the five point stencil */
72: n1 = (int) (sqrt((PetscReal)n) + 0.001);
73: if (n1*n1 - n) SETERRQ(PETSC_ERR_ARG_WRONG,"sqrt(n) must be a positive interger!");
74: for (i=0; i<n1; i++) {
75: for (j=0; j<n1; j++) {
76: Ii = j + n1*i;
77: if (i>0) {
78: J = Ii - n1;
79: MatSetValues(A,1,&Ii,1,&J,&neg_one,INSERT_VALUES);
80: MatSetValues(sA,1,&Ii,1,&J,&neg_one,INSERT_VALUES);
81: }
82: if (i<n1-1) {
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 (j>0) {
88: J = Ii - 1;
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<n1-1) {
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: MatSetValues(A,1,&Ii,1,&Ii,&four,INSERT_VALUES);
98: MatSetValues(sA,1,&Ii,1,&Ii,&four,INSERT_VALUES);
99: }
100: }
101: }
102: }
103: else { /* bs > 1 */
104: for (block=0; block<n/bs; block++){
105: /* diagonal blocks */
106: value[0] = -1.0; value[1] = 4.0; value[2] = -1.0;
107: for (i=1+block*bs; i<bs-1+block*bs; i++) {
108: col[0] = i-1; col[1] = i; col[2] = i+1;
109: MatSetValues(A,1,&i,3,col,value,INSERT_VALUES);
110: MatSetValues(sA,1,&i,3,col,value,INSERT_VALUES);
111: }
112: i = bs - 1+block*bs; col[0] = bs - 2+block*bs; col[1] = bs - 1+block*bs;
113: value[0]=-1.0; value[1]=4.0;
114: MatSetValues(A,1,&i,2,col,value,INSERT_VALUES);
115: MatSetValues(sA,1,&i,2,col,value,INSERT_VALUES);
117: i = 0+block*bs; col[0] = 0+block*bs; col[1] = 1+block*bs;
118: value[0]=4.0; value[1] = -1.0;
119: MatSetValues(A,1,&i,2,col,value,INSERT_VALUES);
120: MatSetValues(sA,1,&i,2,col,value,INSERT_VALUES);
121: }
122: /* off-diagonal blocks */
123: value[0]=-1.0;
124: for (i=0; i<(n/bs-1)*bs; i++){
125: col[0]=i+bs;
126: MatSetValues(A,1,&i,1,col,value,INSERT_VALUES);
127: MatSetValues(sA,1,&i,1,col,value,INSERT_VALUES);
128: col[0]=i; row=i+bs;
129: MatSetValues(A,1,&row,1,col,value,INSERT_VALUES);
130: MatSetValues(sA,1,&row,1,col,value,INSERT_VALUES);
131: }
132: }
133: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
134: MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
135: /* PetscPrintf(PETSC_COMM_SELF,"\n The Matrix: \n");
136: MatView(A, PETSC_VIEWER_DRAW_WORLD);
137: MatView(A, PETSC_VIEWER_STDOUT_WORLD); */
139: MatAssemblyBegin(sA,MAT_FINAL_ASSEMBLY);
140: MatAssemblyEnd(sA,MAT_FINAL_ASSEMBLY);
141: /* PetscPrintf(PETSC_COMM_SELF,"\n Symmetric Part of Matrix: \n");
142: MatView(sA, PETSC_VIEWER_DRAW_WORLD);
143: MatView(sA, PETSC_VIEWER_STDOUT_WORLD);
144: */
146: /* Test MatDuplicate() */
147: MatNorm(A,NORM_FROBENIUS,&norm1);
148: MatDuplicate(sA,MAT_COPY_VALUES,&sB);
149: MatEqual(sA,sB,&equal);
150: if (!equal) SETERRQ(PETSC_ERR_ARG_NOTSAMETYPE,"Error in MatDuplicate()");
152: /* Test MatNorm() */
153: MatNorm(A,NORM_FROBENIUS,&norm1);
154: MatNorm(sB,NORM_FROBENIUS,&norm2);
155: rnorm = PetscAbsScalar(norm1-norm2)/norm2;
156: if (rnorm > tol){
157: PetscPrintf(PETSC_COMM_SELF,"Error: MatNorm_FROBENIUS, NormA=%16.14e NormsB=%16.14e\n",norm1,norm2);
158: }
159: MatNorm(A,NORM_INFINITY,&norm1);
160: MatNorm(sB,NORM_INFINITY,&norm2);
161: rnorm = PetscAbsScalar(norm1-norm2)/norm2;
162: if (rnorm > tol){
163: PetscPrintf(PETSC_COMM_SELF,"Error: MatNorm_INFINITY(), NormA=%16.14e NormsB=%16.14e\n",norm1,norm2);
164: }
165: MatNorm(A,NORM_1,&norm1);
166: MatNorm(sB,NORM_1,&norm2);
167: rnorm = PetscAbsScalar(norm1-norm2)/norm2;
168: if (rnorm > tol){
169: PetscPrintf(PETSC_COMM_SELF,"Error: MatNorm_INFINITY(), NormA=%16.14e NormsB=%16.14e\n",norm1,norm2);
170: }
172: /* Test MatGetInfo(), MatGetSize(), MatGetBlockSize() */
173: MatGetInfo(A,MAT_LOCAL,&minfo1);
174: MatGetInfo(sB,MAT_LOCAL,&minfo2);
175: /*
176: printf("matrix nonzeros (BAIJ format) = %d, allocated nonzeros= %d\n", (int)minfo1.nz_used,(int)minfo1.nz_allocated);
177: printf("matrix nonzeros(SBAIJ format) = %d, allocated nonzeros= %d\n", (int)minfo2.nz_used,(int)minfo2.nz_allocated);
178: */
179: i = (int) (minfo1.nz_used - minfo2.nz_used);
180: j = (int) (minfo2.nz_allocated - minfo2.nz_used);
181: if (i<0 || j<0) {
182: PetscPrintf(PETSC_COMM_SELF,"Error: MatGetInfo()\n");
183: }
185: MatGetSize(A,&Ii,&J);
186: MatGetSize(sB,&i,&j);
187: if (i-Ii || j-J) {
188: PetscPrintf(PETSC_COMM_SELF,"Error: MatGetSize()\n");
189: }
190:
191: MatGetBlockSize(A, &Ii);
192: MatGetBlockSize(sB, &i);
193: if (i-Ii){
194: PetscPrintf(PETSC_COMM_SELF,"Error: MatGetBlockSize()\n");
195: }
197: PetscRandomCreate(PETSC_COMM_SELF,&rdm);
198: PetscRandomSetFromOptions(rdm);
199: VecCreateSeq(PETSC_COMM_SELF,n,&x);
200: VecDuplicate(x,&s1);
201: VecDuplicate(x,&s2);
202: VecDuplicate(x,&y);
203: VecDuplicate(x,&b);
204: VecSetRandom(x,rdm);
206: /* Test MatDiagonalScale(), MatGetDiagonal(), MatScale() */
207: MatDiagonalScale(A,x,x);
208: MatDiagonalScale(sB,x,x);
209: MatMultEqual(A,sB,10,&equal);
210: if (!equal) SETERRQ(PETSC_ERR_ARG_NOTSAMETYPE,"Error in MatDiagonalScale");
212: MatGetDiagonal(A,s1);
213: MatGetDiagonal(sB,s2);
214: VecAXPY(s2,neg_one,s1);
215: VecNorm(s2,NORM_1,&norm1);
216: if ( norm1>tol) {
217: PetscPrintf(PETSC_COMM_SELF,"Error:MatGetDiagonal(), ||s1-s2||=%G\n",norm1);
218: }
220: MatScale(A,alpha);
221: MatScale(sB,alpha);
223: /* Test MatGetRowMax() */
224: MatGetRowMax(A,s1);
225: MatGetRowMax(sB,s2);
226: VecNorm(s1,NORM_1,&norm1);
227: VecNorm(s2,NORM_1,&norm2);
228: norm1 -= norm2;
229: if (norm1<-tol || norm1>tol) {
230: PetscPrintf(PETSC_COMM_SELF,"Error:MatGetRowMax() \n");
231: }
233: /* Test MatMult() */
234: for (i=0; i<40; i++) {
235: VecSetRandom(x,rdm);
236: MatMult(A,x,s1);
237: MatMult(sB,x,s2);
238: VecNorm(s1,NORM_1,&norm1);
239: VecNorm(s2,NORM_1,&norm2);
240: norm1 -= norm2;
241: if (norm1<-tol || norm1>tol) {
242: PetscPrintf(PETSC_COMM_SELF,"Error: MatMult(), norm1-norm2: %G\n",norm1);
243: }
244: }
246: /* MatMultAdd() */
247: for (i=0; i<40; i++) {
248: VecSetRandom(x,rdm);
249: VecSetRandom(y,rdm);
250: MatMultAdd(A,x,y,s1);
251: MatMultAdd(sB,x,y,s2);
252: VecNorm(s1,NORM_1,&norm1);
253: VecNorm(s2,NORM_1,&norm2);
254: norm1 -= norm2;
255: if (norm1<-tol || norm1>tol) {
256: PetscPrintf(PETSC_COMM_SELF,"Error:MatMultAdd(), norm1-norm2: %G\n",norm1);
257: }
258: }
260: /* Test MatCholeskyFactor(), MatICCFactor() with natural ordering */
261: MatGetOrdering(A,MATORDERING_NATURAL,&perm,&iscol);
262: ISDestroy(iscol);
263: norm1 = tol;
264: inc = bs;
266: /* initialize factinfo */
267: PetscMemzero(&factinfo,sizeof(MatFactorInfo));
269: for (lf=-1; lf<10; lf += inc){
270: if (lf==-1) { /* Cholesky factor of sB (duplicate sA) */
271: factinfo.fill = 5.0;
272: MatCholeskyFactorSymbolic(sB,perm,&factinfo,&sC);
273: } else if (!doIcc){
274: break;
275: } else { /* incomplete Cholesky factor */
276: factinfo.fill = 5.0;
277: factinfo.levels = lf;
278: MatICCFactorSymbolic(sB,perm,&factinfo,&sC);
279: }
280: MatCholeskyFactorNumeric(sB,&factinfo,&sC);
281: /* MatView(sC, PETSC_VIEWER_DRAW_WORLD); */
283: /* test MatGetDiagonal on numeric factor */
284: /*
285: if (lf == -1) {
286: MatGetDiagonal(sC,s1);
287: printf(" in ex74.c, diag: \n");
288: VecView(s1,PETSC_VIEWER_STDOUT_SELF);
289: }
290: */
292: MatMult(sB,x,b);
294: /* test MatForwardSolve() and MatBackwardSolve() */
295: if (lf == -1){
296: MatForwardSolve(sC,b,s1);
297: MatBackwardSolve(sC,s1,s2);
298: VecAXPY(s2,neg_one,x);
299: VecNorm(s2,NORM_2,&norm2);
300: if (10*norm1 < norm2){
301: PetscPrintf(PETSC_COMM_SELF,"MatForwardSolve and BackwardSolve: Norm of error=%G, bs=%d\n",norm2,bs);
302: }
303: }
305: /* test MatSolve() */
306: MatSolve(sC,b,y);
307: MatDestroy(sC);
308: /* Check the error */
309: VecAXPY(y,neg_one,x);
310: VecNorm(y,NORM_2,&norm2);
311: /* printf("lf: %d, error: %G\n", lf,norm2); */
312: if (10*norm1 < norm2 && lf-inc != -1){
313: PetscPrintf(PETSC_COMM_SELF,"lf=%D, %D, Norm of error=%G, %G\n",lf-inc,lf,norm1,norm2);
314: }
315: norm1 = norm2;
316: if (norm2 < tol && lf != -1) break;
317: }
319: ISDestroy(perm);
321: MatDestroy(A);
322: MatDestroy(sB);
323: MatDestroy(sA);
324: VecDestroy(x);
325: VecDestroy(y);
326: VecDestroy(s1);
327: VecDestroy(s2);
328: VecDestroy(b);
329: PetscRandomDestroy(rdm);
331: PetscFinalize();
332: return 0;
333: }