Actual source code: ex77.c
petsc-3.8.4 2018-03-24
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>
6: int main(int argc,char **args)
7: {
8: Vec x,y,b,s1,s2;
9: Mat A; /* linear system matrix */
10: Mat sA; /* symmetric part of the matrices */
11: PetscInt n,mbs=16,bs=1,nz=3,prob=2,i,j,col[3],row,Ii,J,n1;
12: const PetscInt *ip_ptr;
13: PetscScalar neg_one = -1.0,value[3],alpha=0.1;
14: PetscMPIInt size;
16: IS ip, isrow, iscol;
17: PetscRandom rdm;
18: PetscBool reorder=PETSC_FALSE;
19: MatInfo minfo1,minfo2;
20: PetscReal norm1,norm2,tol=1.e-10;
22: PetscInitialize(&argc,&args,(char*)0,help);if (ierr) return ierr;
23: MPI_Comm_size(PETSC_COMM_WORLD,&size);
24: if (size != 1) SETERRQ(PETSC_COMM_WORLD,1,"This is a uniprocessor example only!");
25: PetscOptionsGetInt(NULL,NULL,"-bs",&bs,NULL);
26: PetscOptionsGetInt(NULL,NULL,"-mbs",&mbs,NULL);
28: n = mbs*bs;
29: ierr=MatCreateSeqBAIJ(PETSC_COMM_WORLD,bs,n,n,nz,NULL, &A);
30: ierr=MatCreateSeqSBAIJ(PETSC_COMM_WORLD,bs,n,n,nz,NULL, &sA);
32: /* Test MatGetOwnershipRange() */
33: MatGetOwnershipRange(A,&Ii,&J);
34: MatGetOwnershipRange(sA,&i,&j);
35: if (i-Ii || j-J) {
36: PetscPrintf(PETSC_COMM_SELF,"Error: MatGetOwnershipRange() in MatSBAIJ format\n");
37: }
39: /* Assemble matrix */
40: if (bs == 1) {
41: PetscOptionsGetInt(NULL,NULL,"-test_problem",&prob,NULL);
42: if (prob == 1) { /* tridiagonal matrix */
43: value[0] = -1.0; value[1] = 2.0; value[2] = -1.0;
44: for (i=1; i<n-1; i++) {
45: col[0] = i-1; col[1] = i; col[2] = i+1;
46: MatSetValues(A,1,&i,3,col,value,INSERT_VALUES);
47: MatSetValues(sA,1,&i,3,col,value,INSERT_VALUES);
48: }
49: 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;
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: } else if (prob ==2) { /* matrix for the five point stencil */
61: n1 = (PetscInt) (PetscSqrtReal((PetscReal)n) + 0.001);
62: if (n1*n1 - n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"sqrt(n) must be a positive integer!");
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: } else { /* bs > 1 */
94: #if defined(DIAGB)
95: for (block=0; block<n/bs; block++) {
96: /* diagonal blocks */
97: value[0] = -1.0; value[1] = 4.0; value[2] = -1.0;
98: for (i=1+block*bs; i<bs-1+block*bs; i++) {
99: col[0] = i-1; col[1] = i; col[2] = i+1;
100: MatSetValues(A,1,&i,3,col,value,INSERT_VALUES);
101: MatSetValues(sA,1,&i,3,col,value,INSERT_VALUES);
102: }
103: 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;
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: }
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: }
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,MATORDERINGNATURAL,&isrow,&iscol);
233: ip = isrow;
235: if (reorder) {
236: IS nip;
237: PetscInt *nip_ptr;
238: PetscMalloc1(mbs,&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,PETSC_COPY_VALUES,&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: }
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 ierr;
268: }