Actual source code: ex77.c
petsc-3.5.4 2015-05-23
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: PetscBool 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(PETSC_COMM_WORLD,1,"This is a uniprocessor example only!");
27: PetscOptionsGetInt(NULL,"-bs",&bs,NULL);
28: PetscOptionsGetInt(NULL,"-mbs",&mbs,NULL);
30: n = mbs*bs;
31: ierr=MatCreateSeqBAIJ(PETSC_COMM_WORLD,bs,n,n,nz,NULL, &A);
32: ierr=MatCreateSeqSBAIJ(PETSC_COMM_WORLD,bs,n,n,nz,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(NULL,"-test_problem",&prob,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;
53: value[0]= 0.1; value[1]=-1; value[2]=2;
54: MatSetValues(A,1,&i,3,col,value,INSERT_VALUES);
55: MatSetValues(sA,1,&i,3,col,value,INSERT_VALUES);
57: i = 0; col[0] = 0; col[1] = 1; col[2]=n-1;
59: value[0] = 2.0; value[1] = -1.0; value[2]=0.1;
60: MatSetValues(A,1,&i,3,col,value,INSERT_VALUES);
61: MatSetValues(sA,1,&i,3,col,value,INSERT_VALUES);
62: } else if (prob ==2) { /* matrix for the five point stencil */
63: n1 = (PetscInt) (PetscSqrtReal((PetscReal)n) + 0.001);
64: if (n1*n1 - n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"sqrt(n) must be a positive interger!");
65: for (i=0; i<n1; i++) {
66: for (j=0; j<n1; j++) {
67: Ii = j + n1*i;
68: if (i>0) {
69: J = Ii - n1;
70: MatSetValues(A,1,&Ii,1,&J,&neg_one,INSERT_VALUES);
71: MatSetValues(sA,1,&Ii,1,&J,&neg_one,INSERT_VALUES);
72: }
73: if (i<n1-1) {
74: J = Ii + n1;
75: MatSetValues(A,1,&Ii,1,&J,&neg_one,INSERT_VALUES);
76: MatSetValues(sA,1,&Ii,1,&J,&neg_one,INSERT_VALUES);
77: }
78: if (j>0) {
79: J = Ii - 1;
80: MatSetValues(A,1,&Ii,1,&J,&neg_one,INSERT_VALUES);
81: MatSetValues(sA,1,&Ii,1,&J,&neg_one,INSERT_VALUES);
82: }
83: if (j<n1-1) {
84: J = Ii + 1;
85: MatSetValues(A,1,&Ii,1,&J,&neg_one,INSERT_VALUES);
86: MatSetValues(sA,1,&Ii,1,&J,&neg_one,INSERT_VALUES);
87: }
88: /*
89: MatSetValues(A,1,&I,1,&I,&four,INSERT_VALUES);
90: MatSetValues(sA,1,&I,1,&I,&four,INSERT_VALUES);
91: */
92: }
93: }
94: }
95: } else { /* bs > 1 */
96: #if defined(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;
107: value[0]=-1.0; value[1]=4.0;
108: MatSetValues(A,1,&i,2,col,value,INSERT_VALUES);
109: MatSetValues(sA,1,&i,2,col,value,INSERT_VALUES);
111: i = 0+block*bs; col[0] = 0+block*bs; col[1] = 1+block*bs;
113: value[0]=4.0; value[1] = -1.0;
114: MatSetValues(A,1,&i,2,col,value,INSERT_VALUES);
115: MatSetValues(sA,1,&i,2,col,value,INSERT_VALUES);
116: }
117: #endif
118: /* off-diagonal blocks */
119: value[0]=-1.0;
120: for (i=0; i<(n/bs-1)*bs; i++) {
121: col[0]=i+bs;
122: MatSetValues(A,1,&i,1,col,value,INSERT_VALUES);
123: MatSetValues(sA,1,&i,1,col,value,INSERT_VALUES);
124: col[0]=i; row=i+bs;
125: MatSetValues(A,1,&row,1,col,value,INSERT_VALUES);
126: MatSetValues(sA,1,&row,1,col,value,INSERT_VALUES);
127: }
128: }
129: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
130: MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
131: /* PetscPrintf(PETSC_COMM_SELF,"\n The Matrix: \n");
132: MatView(A, VIEWER_DRAW_WORLD);
133: MatView(A, VIEWER_STDOUT_WORLD); */
135: MatAssemblyBegin(sA,MAT_FINAL_ASSEMBLY);
136: MatAssemblyEnd(sA,MAT_FINAL_ASSEMBLY);
137: /* PetscPrintf(PETSC_COMM_SELF,"\n Symmetric Part of Matrix: \n");
138: MatView(sA, VIEWER_DRAW_WORLD);
139: MatView(sA, VIEWER_STDOUT_WORLD);
140: */
142: /* Test MatNorm() */
143: MatNorm(A,NORM_FROBENIUS,&norm1);
144: MatNorm(sA,NORM_FROBENIUS,&norm2);
145: norm1 -= norm2;
146: if (norm1<-tol || norm1>tol) {
147: PetscPrintf(PETSC_COMM_SELF,"Error: MatNorm(), fnorm1-fnorm2=%16.14e\n",norm1);
148: }
149: MatNorm(A,NORM_INFINITY,&norm1);
150: MatNorm(sA,NORM_INFINITY,&norm2);
151: norm1 -= norm2;
152: if (norm1<-tol || norm1>tol) {
153: PetscPrintf(PETSC_COMM_SELF,"Error: MatNorm(), inf_norm1-inf_norm2=%16.14e\n",norm1);
154: }
156: /* Test MatGetInfo(), MatGetSize(), MatGetBlockSize() */
157: MatGetInfo(A,MAT_LOCAL,&minfo1);
158: MatGetInfo(sA,MAT_LOCAL,&minfo2);
159: /*
160: printf("matrix nonzeros (BAIJ format) = %d, allocated nonzeros= %d\n", (int)minfo1.nz_used,(int)minfo1.nz_allocated);
161: printf("matrix nonzeros(SBAIJ format) = %d, allocated nonzeros= %d\n", (int)minfo2.nz_used,(int)minfo2.nz_allocated);
162: */
163: i = (int) (minfo1.nz_used - minfo2.nz_used);
164: j = (int) (minfo1.nz_allocated - minfo2.nz_allocated);
165: if (i<0 || j<0) {
166: PetscPrintf(PETSC_COMM_SELF,"Error: MatGetInfo()\n");
167: }
169: MatGetSize(A,&Ii,&J);
170: MatGetSize(sA,&i,&j);
171: if (i-Ii || j-J) {
172: PetscPrintf(PETSC_COMM_SELF,"Error: MatGetSize()\n");
173: }
175: MatGetBlockSize(A, &Ii);
176: MatGetBlockSize(sA, &i);
177: if (i-Ii) {
178: PetscPrintf(PETSC_COMM_SELF,"Error: MatGetBlockSize()\n");
179: }
181: /* Test MatDiagonalScale(), MatGetDiagonal(), MatScale() */
182: PetscRandomCreate(PETSC_COMM_SELF,&rdm);
183: PetscRandomSetFromOptions(rdm);
184: VecCreateSeq(PETSC_COMM_SELF,n,&x);
185: VecDuplicate(x,&s1);
186: VecDuplicate(x,&s2);
187: VecDuplicate(x,&y);
188: VecDuplicate(x,&b);
190: VecSetRandom(x,rdm);
192: MatDiagonalScale(A,x,x);
193: MatDiagonalScale(sA,x,x);
195: MatGetDiagonal(A,s1);
196: MatGetDiagonal(sA,s2);
197: VecNorm(s1,NORM_1,&norm1);
198: VecNorm(s2,NORM_1,&norm2);
199: norm1 -= norm2;
200: if (norm1<-tol || norm1>tol) {
201: PetscPrintf(PETSC_COMM_SELF,"Error:MatGetDiagonal() \n");
202: }
204: MatScale(A,alpha);
205: MatScale(sA,alpha);
207: /* Test MatMult(), MatMultAdd() */
208: for (i=0; i<40; i++) {
209: VecSetRandom(x,rdm);
210: MatMult(A,x,s1);
211: MatMult(sA,x,s2);
212: VecNorm(s1,NORM_1,&norm1);
213: VecNorm(s2,NORM_1,&norm2);
214: norm1 -= norm2;
215: if (norm1<-tol || norm1>tol) {
216: PetscPrintf(PETSC_COMM_SELF,"Error: MatMult(), MatDiagonalScale() or MatScale()\n");
217: }
218: }
220: for (i=0; i<40; i++) {
221: VecSetRandom(x,rdm);
222: VecSetRandom(y,rdm);
223: MatMultAdd(A,x,y,s1);
224: MatMultAdd(sA,x,y,s2);
225: VecNorm(s1,NORM_1,&norm1);
226: VecNorm(s2,NORM_1,&norm2);
227: norm1 -= norm2;
228: if (norm1<-tol || norm1>tol) {
229: PetscPrintf(PETSC_COMM_SELF,"Error:MatMultAdd(), MatDiagonalScale() or MatScale() \n");
230: }
231: }
233: /* Test MatReordering() */
234: MatGetOrdering(A,MATORDERINGNATURAL,&isrow,&iscol);
235: ip = isrow;
237: if (reorder) {
238: IS nip;
239: PetscInt *nip_ptr;
240: PetscMalloc1(mbs,&nip_ptr);
241: ISGetIndices(ip,&ip_ptr);
242: PetscMemcpy(nip_ptr,ip_ptr,mbs*sizeof(PetscInt));
243: i = nip_ptr[1]; nip_ptr[1] = nip_ptr[mbs-2]; nip_ptr[mbs-2] = i;
244: i = nip_ptr[0]; nip_ptr[0] = nip_ptr[mbs-1]; nip_ptr[mbs-1] = i;
245: ISRestoreIndices(ip,&ip_ptr);
246: ISCreateGeneral(PETSC_COMM_SELF,mbs,nip_ptr,PETSC_COPY_VALUES,&nip);
247: PetscFree(nip_ptr);
249: MatReorderingSeqSBAIJ(sA, ip);
250: ISDestroy(&nip);
251: /* ISView(ip, VIEWER_STDOUT_SELF);
252: MatView(sA,VIEWER_DRAW_SELF); */
253: }
255: ISDestroy(&iscol);
256: /* ISDestroy(&isrow);*/
258: ISDestroy(&isrow);
259: MatDestroy(&A);
260: MatDestroy(&sA);
261: VecDestroy(&x);
262: VecDestroy(&y);
263: VecDestroy(&s1);
264: VecDestroy(&s2);
265: VecDestroy(&b);
266: PetscRandomDestroy(&rdm);
268: PetscFinalize();
269: return 0;
270: }