Actual source code: ex77.c
2: static char help[] = "Tests the various sequential routines in MatSBAIJ format. Same as ex74.c except diagonal entries of the matrices are zeros.\n";
4: #include petscmat.h
8: int main(int argc,char **args)
9: {
10: Vec x,y,b,s1,s2;
11: Mat A; /* linear system matrix */
12: Mat sA; /* symmetric part of the matrices */
13: PetscInt n,mbs=16,bs=1,nz=3,prob=2,i,j,col[3],row,Ii,J,n1,*ip_ptr;
14: PetscScalar neg_one = -1.0,value[3],alpha=0.1;
15: PetscMPIInt size;
17: IS ip, isrow, iscol;
18: PetscRandom rdm;
19: PetscTruth reorder=PETSC_FALSE;
20: MatInfo minfo1,minfo2;
21: PetscReal norm1,norm2,tol=1.e-10;
23: PetscInitialize(&argc,&args,(char *)0,help);
24: MPI_Comm_size(PETSC_COMM_WORLD,&size);
25: if (size != 1) SETERRQ(1,"This is a uniprocessor example only!");
26: PetscOptionsGetInt(PETSC_NULL,"-bs",&bs,PETSC_NULL);
27: PetscOptionsGetInt(PETSC_NULL,"-mbs",&mbs,PETSC_NULL);
29: n = mbs*bs;
30: ierr=MatCreateSeqBAIJ(PETSC_COMM_WORLD,bs,n,n,nz,PETSC_NULL, &A);
31: ierr=MatCreateSeqSBAIJ(PETSC_COMM_WORLD,bs,n,n,nz,PETSC_NULL, &sA);
33: /* Test MatGetOwnershipRange() */
34: MatGetOwnershipRange(A,&Ii,&J);
35: MatGetOwnershipRange(sA,&i,&j);
36: if (i-Ii || j-J){
37: PetscPrintf(PETSC_COMM_SELF,"Error: MatGetOwnershipRange() in MatSBAIJ format\n");
38: }
40: /* Assemble matrix */
41: if (bs == 1){
42: PetscOptionsGetInt(PETSC_NULL,"-test_problem",&prob,PETSC_NULL);
43: if (prob == 1){ /* tridiagonal matrix */
44: value[0] = -1.0; value[1] = 2.0; value[2] = -1.0;
45: for (i=1; i<n-1; i++) {
46: col[0] = i-1; col[1] = i; col[2] = i+1;
47: MatSetValues(A,1,&i,3,col,value,INSERT_VALUES);
48: MatSetValues(sA,1,&i,3,col,value,INSERT_VALUES);
49: }
50: i = n - 1; col[0]=0; col[1] = n - 2; col[2] = n - 1;
51: value[0]= 0.1; value[1]=-1; value[2]=2;
52: MatSetValues(A,1,&i,3,col,value,INSERT_VALUES);
53: MatSetValues(sA,1,&i,3,col,value,INSERT_VALUES);
55: i = 0; col[0] = 0; col[1] = 1; col[2]=n-1;
56: value[0] = 2.0; value[1] = -1.0; value[2]=0.1;
57: MatSetValues(A,1,&i,3,col,value,INSERT_VALUES);
58: MatSetValues(sA,1,&i,3,col,value,INSERT_VALUES);
59: }
60: else if (prob ==2){ /* matrix for the five point stencil */
61: n1 = (int) (sqrt((PetscReal)n) + 0.001);
62: if (n1*n1 - n) SETERRQ(PETSC_ERR_ARG_WRONG,"sqrt(n) must be a positive interger!");
63: for (i=0; i<n1; i++) {
64: for (j=0; j<n1; j++) {
65: Ii = j + n1*i;
66: if (i>0) {
67: J = Ii - n1;
68: MatSetValues(A,1,&Ii,1,&J,&neg_one,INSERT_VALUES);
69: MatSetValues(sA,1,&Ii,1,&J,&neg_one,INSERT_VALUES);
70: }
71: if (i<n1-1) {
72: J = Ii + n1;
73: MatSetValues(A,1,&Ii,1,&J,&neg_one,INSERT_VALUES);
74: MatSetValues(sA,1,&Ii,1,&J,&neg_one,INSERT_VALUES);
75: }
76: if (j>0) {
77: J = Ii - 1;
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 (j<n1-1) {
82: J = Ii + 1;
83: MatSetValues(A,1,&Ii,1,&J,&neg_one,INSERT_VALUES);
84: MatSetValues(sA,1,&Ii,1,&J,&neg_one,INSERT_VALUES);
85: }
86: /*
87: MatSetValues(A,1,&I,1,&I,&four,INSERT_VALUES);
88: MatSetValues(sA,1,&I,1,&I,&four,INSERT_VALUES);
89: */
90: }
91: }
92: }
93: }
94: else { /* bs > 1 */
95: #ifdef DIAGB
96: for (block=0; block<n/bs; block++){
97: /* diagonal blocks */
98: value[0] = -1.0; value[1] = 4.0; value[2] = -1.0;
99: for (i=1+block*bs; i<bs-1+block*bs; i++) {
100: col[0] = i-1; col[1] = i; col[2] = i+1;
101: MatSetValues(A,1,&i,3,col,value,INSERT_VALUES);
102: MatSetValues(sA,1,&i,3,col,value,INSERT_VALUES);
103: }
104: i = bs - 1+block*bs; col[0] = bs - 2+block*bs; col[1] = bs - 1+block*bs;
105: value[0]=-1.0; value[1]=4.0;
106: MatSetValues(A,1,&i,2,col,value,INSERT_VALUES);
107: MatSetValues(sA,1,&i,2,col,value,INSERT_VALUES);
109: i = 0+block*bs; col[0] = 0+block*bs; col[1] = 1+block*bs;
110: value[0]=4.0; value[1] = -1.0;
111: MatSetValues(A,1,&i,2,col,value,INSERT_VALUES);
112: MatSetValues(sA,1,&i,2,col,value,INSERT_VALUES);
113: }
114: #endif
115: /* off-diagonal blocks */
116: value[0]=-1.0;
117: for (i=0; i<(n/bs-1)*bs; i++){
118: col[0]=i+bs;
119: MatSetValues(A,1,&i,1,col,value,INSERT_VALUES);
120: MatSetValues(sA,1,&i,1,col,value,INSERT_VALUES);
121: col[0]=i; row=i+bs;
122: MatSetValues(A,1,&row,1,col,value,INSERT_VALUES);
123: MatSetValues(sA,1,&row,1,col,value,INSERT_VALUES);
124: }
125: }
126: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
127: MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
128: /* PetscPrintf(PETSC_COMM_SELF,"\n The Matrix: \n");
129: MatView(A, VIEWER_DRAW_WORLD);
130: MatView(A, VIEWER_STDOUT_WORLD); */
132: MatAssemblyBegin(sA,MAT_FINAL_ASSEMBLY);
133: MatAssemblyEnd(sA,MAT_FINAL_ASSEMBLY);
134: /* PetscPrintf(PETSC_COMM_SELF,"\n Symmetric Part of Matrix: \n");
135: MatView(sA, VIEWER_DRAW_WORLD);
136: MatView(sA, VIEWER_STDOUT_WORLD);
137: */
139: /* Test MatNorm() */
140: MatNorm(A,NORM_FROBENIUS,&norm1);
141: MatNorm(sA,NORM_FROBENIUS,&norm2);
142: norm1 -= norm2;
143: if (norm1<-tol || norm1>tol){
144: PetscPrintf(PETSC_COMM_SELF,"Error: MatNorm(), fnorm1-fnorm2=%16.14e\n",norm1);
145: }
146: MatNorm(A,NORM_INFINITY,&norm1);
147: MatNorm(sA,NORM_INFINITY,&norm2);
148: norm1 -= norm2;
149: if (norm1<-tol || norm1>tol){
150: PetscPrintf(PETSC_COMM_SELF,"Error: MatNorm(), inf_norm1-inf_norm2=%16.14e\n",norm1);
151: }
152:
153: /* Test MatGetInfo(), MatGetSize(), MatGetBlockSize() */
154: MatGetInfo(A,MAT_LOCAL,&minfo1);
155: MatGetInfo(sA,MAT_LOCAL,&minfo2);
156: /*
157: printf("matrix nonzeros (BAIJ format) = %d, allocated nonzeros= %d\n", (int)minfo1.nz_used,(int)minfo1.nz_allocated);
158: printf("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: if (i<0 || j<0) {
163: PetscPrintf(PETSC_COMM_SELF,"Error: MatGetInfo()\n");
164: }
166: MatGetSize(A,&Ii,&J);
167: MatGetSize(sA,&i,&j);
168: if (i-Ii || j-J) {
169: PetscPrintf(PETSC_COMM_SELF,"Error: MatGetSize()\n");
170: }
171:
172: MatGetBlockSize(A, &Ii);
173: MatGetBlockSize(sA, &i);
174: if (i-Ii){
175: PetscPrintf(PETSC_COMM_SELF,"Error: MatGetBlockSize()\n");
176: }
178: /* Test MatDiagonalScale(), MatGetDiagonal(), MatScale() */
179: PetscRandomCreate(PETSC_COMM_SELF,&rdm);
180: PetscRandomSetFromOptions(rdm);
181: VecCreateSeq(PETSC_COMM_SELF,n,&x);
182: VecDuplicate(x,&s1);
183: VecDuplicate(x,&s2);
184: VecDuplicate(x,&y);
185: VecDuplicate(x,&b);
187: VecSetRandom(x,rdm);
189: MatDiagonalScale(A,x,x);
190: MatDiagonalScale(sA,x,x);
192: MatGetDiagonal(A,s1);
193: MatGetDiagonal(sA,s2);
194: VecNorm(s1,NORM_1,&norm1);
195: VecNorm(s2,NORM_1,&norm2);
196: norm1 -= norm2;
197: if (norm1<-tol || norm1>tol) {
198: PetscPrintf(PETSC_COMM_SELF,"Error:MatGetDiagonal() \n");
199: }
201: MatScale(A,alpha);
202: MatScale(sA,alpha);
204: /* Test MatMult(), MatMultAdd() */
205: for (i=0; i<40; i++) {
206: VecSetRandom(x,rdm);
207: MatMult(A,x,s1);
208: MatMult(sA,x,s2);
209: VecNorm(s1,NORM_1,&norm1);
210: VecNorm(s2,NORM_1,&norm2);
211: norm1 -= norm2;
212: if (norm1<-tol || norm1>tol) {
213: PetscPrintf(PETSC_COMM_SELF,"Error: MatMult(), MatDiagonalScale() or MatScale()\n");
214: }
215: }
217: for (i=0; i<40; i++) {
218: VecSetRandom(x,rdm);
219: VecSetRandom(y,rdm);
220: MatMultAdd(A,x,y,s1);
221: MatMultAdd(sA,x,y,s2);
222: VecNorm(s1,NORM_1,&norm1);
223: VecNorm(s2,NORM_1,&norm2);
224: norm1 -= norm2;
225: if (norm1<-tol || norm1>tol) {
226: PetscPrintf(PETSC_COMM_SELF,"Error:MatMultAdd(), MatDiagonalScale() or MatScale() \n");
227: }
228: }
230: /* Test MatReordering() */
231: MatGetOrdering(A,MATORDERING_NATURAL,&isrow,&iscol);
232: ip = isrow;
234: if (reorder){
235: ISGetIndices(ip,&ip_ptr);
236: i = ip_ptr[1]; ip_ptr[1] = ip_ptr[mbs-2]; ip_ptr[mbs-2] = i;
237: i = ip_ptr[0]; ip_ptr[0] = ip_ptr[mbs-1]; ip_ptr[mbs-1] = i;
238: ierr= ISRestoreIndices(ip,&ip_ptr);
240: MatReorderingSeqSBAIJ(sA, ip);
241: /* ISView(ip, VIEWER_STDOUT_SELF);
242: MatView(sA,VIEWER_DRAW_SELF); */
243: }
244:
245: ISDestroy(iscol);
246: /* ISDestroy(isrow);*/
248: ISDestroy(isrow);
249: MatDestroy(A);
250: MatDestroy(sA);
251: VecDestroy(x);
252: VecDestroy(y);
253: VecDestroy(s1);
254: VecDestroy(s2);
255: VecDestroy(b);
256: PetscRandomDestroy(rdm);
258: PetscFinalize();
259: return 0;
260: }