Actual source code: ex77.c
1: /*$Id: ex77.c,v 1.12 2001/08/07 21:30:08 bsmith Exp $*/
3: static char help[] = "Tests the various sequential routines in MatSBAIJ format. Same as ex74.c except diagonal entries of the matrices are zeros.\n";
5: #include petscmat.h
7: /* extern int MatReorderingSeqSBAIJ(Mat,IS); */
11: int main(int argc,char **args)
12: {
13: Vec x,y,b,s1,s2;
14: Mat A; /* linear system matrix */
15: Mat sA,sC; /* symmetric part of the matrices */
17: int n,mbs=16,bs=1,nz=3,prob=2;
18: PetscScalar neg_one = -1.0,value[3],alpha=0.1;
19: int ierr,i,j,col[3],size,row,I,J,n1,*ip_ptr;
20: IS ip, isrow, iscol;
21: PetscRandom rand;
22: PetscTruth reorder=PETSC_FALSE,getrow=PETSC_FALSE;
23: MatInfo minfo1,minfo2;
24: PetscScalar *vr1,*vr2,*vr1_wk,*vr2_wk;
25: int *cols1,*cols2;
26: PetscReal norm1,norm2,tol=1.e-10;
28: PetscInitialize(&argc,&args,(char *)0,help);
29: MPI_Comm_size(PETSC_COMM_WORLD,&size);
30: if (size != 1) SETERRQ(1,"This is a uniprocessor example only!");
31: PetscOptionsGetInt(PETSC_NULL,"-bs",&bs,PETSC_NULL);
32: PetscOptionsGetInt(PETSC_NULL,"-mbs",&mbs,PETSC_NULL);
34: n = mbs*bs;
35: ierr=MatCreateSeqBAIJ(PETSC_COMM_WORLD,bs,n,n,nz,PETSC_NULL, &A);
36: ierr=MatCreateSeqSBAIJ(PETSC_COMM_WORLD,bs,n,n,nz,PETSC_NULL, &sA);
38: /* Test MatGetOwnershipRange() */
39: MatGetOwnershipRange(A,&I,&J);
40: MatGetOwnershipRange(sA,&i,&j);
41: if (i-I || j-J){
42: PetscPrintf(PETSC_COMM_SELF,"Error: MatGetOwnershipRange() in MatSBAIJ format\n");
43: }
45: /* Assemble matrix */
46: if (bs == 1){
47: PetscOptionsGetInt(PETSC_NULL,"-test_problem",&prob,PETSC_NULL);
48: if (prob == 1){ /* tridiagonal matrix */
49: value[0] = -1.0; value[1] = 2.0; value[2] = -1.0;
50: for (i=1; i<n-1; i++) {
51: col[0] = i-1; col[1] = i; col[2] = i+1;
52: MatSetValues(A,1,&i,3,col,value,INSERT_VALUES);
53: MatSetValues(sA,1,&i,3,col,value,INSERT_VALUES);
54: }
55: i = n - 1; col[0]=0; col[1] = n - 2; col[2] = n - 1;
56: value[0]= 0.1; value[1]=-1; value[2]=2;
57: MatSetValues(A,1,&i,3,col,value,INSERT_VALUES);
58: MatSetValues(sA,1,&i,3,col,value,INSERT_VALUES);
60: i = 0; col[0] = 0; col[1] = 1; col[2]=n-1;
61: value[0] = 2.0; value[1] = -1.0; value[2]=0.1;
62: MatSetValues(A,1,&i,3,col,value,INSERT_VALUES);
63: MatSetValues(sA,1,&i,3,col,value,INSERT_VALUES);
64: }
65: else if (prob ==2){ /* matrix for the five point stencil */
66: n1 = (int) (sqrt((PetscReal)n) + 0.001);
67: if (n1*n1 - n) SETERRQ(PETSC_ERR_ARG_WRONG,"sqrt(n) must be a positive interger!");
68: for (i=0; i<n1; i++) {
69: for (j=0; j<n1; j++) {
70: I = j + n1*i;
71: if (i>0) {
72: J = I - n1;
73: MatSetValues(A,1,&I,1,&J,&neg_one,INSERT_VALUES);
74: MatSetValues(sA,1,&I,1,&J,&neg_one,INSERT_VALUES);
75: }
76: if (i<n1-1) {
77: J = I + n1;
78: MatSetValues(A,1,&I,1,&J,&neg_one,INSERT_VALUES);
79: MatSetValues(sA,1,&I,1,&J,&neg_one,INSERT_VALUES);
80: }
81: if (j>0) {
82: J = I - 1;
83: MatSetValues(A,1,&I,1,&J,&neg_one,INSERT_VALUES);
84: MatSetValues(sA,1,&I,1,&J,&neg_one,INSERT_VALUES);
85: }
86: if (j<n1-1) {
87: J = I + 1;
88: MatSetValues(A,1,&I,1,&J,&neg_one,INSERT_VALUES);
89: MatSetValues(sA,1,&I,1,&J,&neg_one,INSERT_VALUES);
90: }
91: /*
92: MatSetValues(A,1,&I,1,&I,&four,INSERT_VALUES);
93: MatSetValues(sA,1,&I,1,&I,&four,INSERT_VALUES);
94: */
95: }
96: }
97: }
98: }
99: else { /* bs > 1 */
100: #ifdef DIAGB
101: for (block=0; block<n/bs; block++){
102: /* diagonal blocks */
103: value[0] = -1.0; value[1] = 4.0; value[2] = -1.0;
104: for (i=1+block*bs; i<bs-1+block*bs; i++) {
105: col[0] = i-1; col[1] = i; col[2] = i+1;
106: MatSetValues(A,1,&i,3,col,value,INSERT_VALUES);
107: MatSetValues(sA,1,&i,3,col,value,INSERT_VALUES);
108: }
109: i = bs - 1+block*bs; col[0] = bs - 2+block*bs; col[1] = bs - 1+block*bs;
110: value[0]=-1.0; value[1]=4.0;
111: MatSetValues(A,1,&i,2,col,value,INSERT_VALUES);
112: MatSetValues(sA,1,&i,2,col,value,INSERT_VALUES);
114: i = 0+block*bs; col[0] = 0+block*bs; col[1] = 1+block*bs;
115: value[0]=4.0; value[1] = -1.0;
116: MatSetValues(A,1,&i,2,col,value,INSERT_VALUES);
117: MatSetValues(sA,1,&i,2,col,value,INSERT_VALUES);
118: }
119: #endif
120: /* off-diagonal blocks */
121: value[0]=-1.0;
122: for (i=0; i<(n/bs-1)*bs; i++){
123: col[0]=i+bs;
124: MatSetValues(A,1,&i,1,col,value,INSERT_VALUES);
125: MatSetValues(sA,1,&i,1,col,value,INSERT_VALUES);
126: col[0]=i; row=i+bs;
127: MatSetValues(A,1,&row,1,col,value,INSERT_VALUES);
128: MatSetValues(sA,1,&row,1,col,value,INSERT_VALUES);
129: }
130: }
131: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
132: MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
133: /* PetscPrintf(PETSC_COMM_SELF,"\n The Matrix: \n");
134: MatView(A, VIEWER_DRAW_WORLD);
135: MatView(A, VIEWER_STDOUT_WORLD); */
137: MatAssemblyBegin(sA,MAT_FINAL_ASSEMBLY);
138: MatAssemblyEnd(sA,MAT_FINAL_ASSEMBLY);
139: /* PetscPrintf(PETSC_COMM_SELF,"\n Symmetric Part of Matrix: \n");
140: MatView(sA, VIEWER_DRAW_WORLD);
141: MatView(sA, VIEWER_STDOUT_WORLD);
142: */
144: /* Test MatNorm() */
145: MatNorm(A,NORM_FROBENIUS,&norm1);
146: MatNorm(sA,NORM_FROBENIUS,&norm2);
147: norm1 -= norm2;
148: if (norm1<-tol || norm1>tol){
149: PetscPrintf(PETSC_COMM_SELF,"Error: MatNorm(), fnorm1-fnorm2=%16.14e\n",norm1);
150: }
151: MatNorm(A,NORM_INFINITY,&norm1);
152: MatNorm(sA,NORM_INFINITY,&norm2);
153: norm1 -= norm2;
154: if (norm1<-tol || norm1>tol){
155: PetscPrintf(PETSC_COMM_SELF,"Error: MatNorm(), inf_norm1-inf_norm2=%16.14e\n",norm1);
156: }
157:
158: /* Test MatGetInfo(), MatGetSize(), MatGetBlockSize() */
159: MatGetInfo(A,MAT_LOCAL,&minfo1);
160: MatGetInfo(sA,MAT_LOCAL,&minfo2);
161: /*
162: printf("matrix nonzeros (BAIJ format) = %d, allocated nonzeros= %d\n", (int)minfo1.nz_used,(int)minfo1.nz_allocated);
163: printf("matrix nonzeros(SBAIJ format) = %d, allocated nonzeros= %d\n", (int)minfo2.nz_used,(int)minfo2.nz_allocated);
164: */
165: i = (int) (minfo1.nz_used - minfo2.nz_used);
166: j = (int) (minfo1.nz_allocated - minfo2.nz_allocated);
167: if (i<0 || j<0) {
168: PetscPrintf(PETSC_COMM_SELF,"Error: MatGetInfo()\n");
169: }
171: MatGetSize(A,&I,&J);
172: MatGetSize(sA,&i,&j);
173: if (i-I || j-J) {
174: PetscPrintf(PETSC_COMM_SELF,"Error: MatGetSize()\n");
175: }
176:
177: MatGetBlockSize(A, &I);
178: MatGetBlockSize(sA, &i);
179: if (i-I){
180: PetscPrintf(PETSC_COMM_SELF,"Error: MatGetBlockSize()\n");
181: }
183: /* Test MatGetRow() */
184: if (getrow){
185: row = n/2;
186: PetscMalloc(n*sizeof(PetscScalar),&vr1);
187: vr1_wk = vr1;
188: PetscMalloc(n*sizeof(PetscScalar),&vr2);
189: vr2_wk = vr2;
190: MatGetRow(A,row,&J,&cols1,&vr1);
191: vr1_wk += J-1;
192: MatGetRow(sA,row,&j,&cols2,&vr2);
193: vr2_wk += j-1;
194: VecCreateSeq(PETSC_COMM_SELF,j,&x);
195:
196: for (i=j-1; i>-1; i--){
197: VecSetValue(x,i,*vr2_wk - *vr1_wk,INSERT_VALUES);
198: vr2_wk--; vr1_wk--;
199: }
200: VecNorm(x,NORM_1,&norm2);
201: if (norm2<-tol || norm2>tol) {
202: PetscPrintf(PETSC_COMM_SELF,"Error: MatGetRow()\n");
203: }
204: VecDestroy(x);
205: MatRestoreRow(A,row,&J,&cols1,&vr1);
206: MatRestoreRow(sA,row,&j,&cols2,&vr2);
207: PetscFree(vr1);
208: PetscFree(vr2);
210: /* Test GetSubMatrix() */
211: /* get a submatrix consisting of every next block row and column of the original matrix */
212: /* for symm. matrix, iscol=isrow. */
213: PetscMalloc(n*sizeof(IS),&isrow);
214: PetscMalloc(n*sizeof(int),&ip_ptr);
215: j = 0;
216: for (n1=0; n1<mbs; n1 += 2){ /* n1: block row */
217: for (i=0; i<bs; i++) ip_ptr[j++] = n1*bs + i;
218: }
219: ISCreateGeneral(PETSC_COMM_SELF, j, ip_ptr, &isrow);
220: /* ISView(isrow, VIEWER_STDOUT_SELF); */
221:
222: MatGetSubMatrix(sA,isrow,isrow,PETSC_DECIDE,MAT_INITIAL_MATRIX,&sC);
223: ISDestroy(isrow);
224: PetscFree(ip_ptr);
225: printf("sA =\n");
226: MatView(sA,PETSC_VIEWER_STDOUT_WORLD);
227: printf("submatrix of sA =\n");
228: MatView(sC,PETSC_VIEWER_STDOUT_WORLD);
229: MatDestroy(sC);
230: }
232: /* Test MatDiagonalScale(), MatGetDiagonal(), MatScale() */
233: PetscRandomCreate(PETSC_COMM_SELF,RANDOM_DEFAULT,&rand);
234: VecCreateSeq(PETSC_COMM_SELF,n,&x);
235: VecDuplicate(x,&s1);
236: VecDuplicate(x,&s2);
237: VecDuplicate(x,&y);
238: VecDuplicate(x,&b);
240: VecSetRandom(rand,x);
242: MatDiagonalScale(A,x,x);
243: MatDiagonalScale(sA,x,x);
245: MatGetDiagonal(A,s1);
246: MatGetDiagonal(sA,s2);
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:MatGetDiagonal() \n");
252: }
254: MatScale(&alpha,A);
255: MatScale(&alpha,sA);
257: /* Test MatMult(), MatMultAdd() */
258: for (i=0; i<40; i++) {
259: VecSetRandom(rand,x);
260: MatMult(A,x,s1);
261: MatMult(sA,x,s2);
262: VecNorm(s1,NORM_1,&norm1);
263: VecNorm(s2,NORM_1,&norm2);
264: norm1 -= norm2;
265: if (norm1<-tol || norm1>tol) {
266: PetscPrintf(PETSC_COMM_SELF,"Error: MatMult(), MatDiagonalScale() or MatScale()\n");
267: }
268: }
270: for (i=0; i<40; i++) {
271: VecSetRandom(rand,x);
272: VecSetRandom(rand,y);
273: MatMultAdd(A,x,y,s1);
274: MatMultAdd(sA,x,y,s2);
275: VecNorm(s1,NORM_1,&norm1);
276: VecNorm(s2,NORM_1,&norm2);
277: norm1 -= norm2;
278: if (norm1<-tol || norm1>tol) {
279: PetscPrintf(PETSC_COMM_SELF,"Error:MatMultAdd(), MatDiagonalScale() or MatScale() \n");
280: }
281: }
283: /* Test MatReordering() */
284: MatGetOrdering(A,MATORDERING_NATURAL,&isrow,&iscol);
285: ip = isrow;
287: if (reorder){
288: ISGetIndices(ip,&ip_ptr);
289: i = ip_ptr[1]; ip_ptr[1] = ip_ptr[mbs-2]; ip_ptr[mbs-2] = i;
290: i = ip_ptr[0]; ip_ptr[0] = ip_ptr[mbs-1]; ip_ptr[mbs-1] = i;
291: ierr= ISRestoreIndices(ip,&ip_ptr);
293: MatReorderingSeqSBAIJ(sA, ip);
294: /* ISView(ip, VIEWER_STDOUT_SELF);
295: MatView(sA,VIEWER_DRAW_SELF); */
296: }
297:
298: ISDestroy(iscol);
299: /* ISDestroy(isrow);*/
301: ISDestroy(isrow);
302: MatDestroy(A);
303: MatDestroy(sA);
304: VecDestroy(x);
305: VecDestroy(y);
306: VecDestroy(s1);
307: VecDestroy(s2);
308: VecDestroy(b);
309: PetscRandomDestroy(rand);
311: PetscFinalize();
312: return 0;
313: }