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