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