Actual source code: ex77.c
petsc-3.14.6 2021-03-30
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=10*PETSC_SMALL;
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,PETSC_ERR_WRONG_MPI_SIZE,"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: MatCreateSeqBAIJ(PETSC_COMM_WORLD,bs,n,n,nz,NULL, &A);
30: MatSetOption(A,MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_FALSE);
31: MatCreateSeqSBAIJ(PETSC_COMM_WORLD,bs,n,n,nz,NULL, &sA);
32: MatSetOption(sA,MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_FALSE);
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,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 integer!");
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: }
90: }
91: } else { /* bs > 1 */
92: #if defined(DIAGB)
93: for (block=0; block<n/bs; block++) {
94: /* diagonal blocks */
95: value[0] = -1.0; value[1] = 4.0; value[2] = -1.0;
96: for (i=1+block*bs; i<bs-1+block*bs; i++) {
97: col[0] = i-1; col[1] = i; col[2] = i+1;
98: MatSetValues(A,1,&i,3,col,value,INSERT_VALUES);
99: MatSetValues(sA,1,&i,3,col,value,INSERT_VALUES);
100: }
101: i = bs - 1+block*bs; col[0] = bs - 2+block*bs; col[1] = bs - 1+block*bs;
103: value[0]=-1.0; value[1]=4.0;
104: MatSetValues(A,1,&i,2,col,value,INSERT_VALUES);
105: MatSetValues(sA,1,&i,2,col,value,INSERT_VALUES);
107: i = 0+block*bs; col[0] = 0+block*bs; col[1] = 1+block*bs;
109: value[0]=4.0; value[1] = -1.0;
110: MatSetValues(A,1,&i,2,col,value,INSERT_VALUES);
111: MatSetValues(sA,1,&i,2,col,value,INSERT_VALUES);
112: }
113: #endif
114: /* off-diagonal blocks */
115: value[0]=-1.0;
116: for (i=0; i<(n/bs-1)*bs; i++) {
117: col[0]=i+bs;
118: MatSetValues(A,1,&i,1,col,value,INSERT_VALUES);
119: MatSetValues(sA,1,&i,1,col,value,INSERT_VALUES);
120: col[0]=i; row=i+bs;
121: MatSetValues(A,1,&row,1,col,value,INSERT_VALUES);
122: MatSetValues(sA,1,&row,1,col,value,INSERT_VALUES);
123: }
124: }
125: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
126: MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
128: MatAssemblyBegin(sA,MAT_FINAL_ASSEMBLY);
129: MatAssemblyEnd(sA,MAT_FINAL_ASSEMBLY);
131: /* Test MatNorm() */
132: MatNorm(A,NORM_FROBENIUS,&norm1);
133: MatNorm(sA,NORM_FROBENIUS,&norm2);
134: norm1 -= norm2;
135: if (norm1<-tol || norm1>tol) {
136: PetscPrintf(PETSC_COMM_SELF,"Error: MatNorm(), fnorm1-fnorm2=%16.14e\n",norm1);
137: }
138: MatNorm(A,NORM_INFINITY,&norm1);
139: MatNorm(sA,NORM_INFINITY,&norm2);
140: norm1 -= norm2;
141: if (norm1<-tol || norm1>tol) {
142: PetscPrintf(PETSC_COMM_SELF,"Error: MatNorm(), inf_norm1-inf_norm2=%16.14e\n",norm1);
143: }
145: /* Test MatGetInfo(), MatGetSize(), MatGetBlockSize() */
146: MatGetInfo(A,MAT_LOCAL,&minfo1);
147: MatGetInfo(sA,MAT_LOCAL,&minfo2);
148: i = (int) (minfo1.nz_used - minfo2.nz_used);
149: j = (int) (minfo1.nz_allocated - minfo2.nz_allocated);
150: if (i<0 || j<0) {
151: PetscPrintf(PETSC_COMM_SELF,"Error: MatGetInfo()\n");
152: }
154: MatGetSize(A,&Ii,&J);
155: MatGetSize(sA,&i,&j);
156: if (i-Ii || j-J) {
157: PetscPrintf(PETSC_COMM_SELF,"Error: MatGetSize()\n");
158: }
160: MatGetBlockSize(A, &Ii);
161: MatGetBlockSize(sA, &i);
162: if (i-Ii) {
163: PetscPrintf(PETSC_COMM_SELF,"Error: MatGetBlockSize()\n");
164: }
166: /* Test MatDiagonalScale(), MatGetDiagonal(), MatScale() */
167: PetscRandomCreate(PETSC_COMM_SELF,&rdm);
168: PetscRandomSetFromOptions(rdm);
169: VecCreateSeq(PETSC_COMM_SELF,n,&x);
170: VecDuplicate(x,&s1);
171: VecDuplicate(x,&s2);
172: VecDuplicate(x,&y);
173: VecDuplicate(x,&b);
175: VecSetRandom(x,rdm);
177: MatDiagonalScale(A,x,x);
178: MatDiagonalScale(sA,x,x);
180: MatGetDiagonal(A,s1);
181: MatGetDiagonal(sA,s2);
182: VecNorm(s1,NORM_1,&norm1);
183: VecNorm(s2,NORM_1,&norm2);
184: norm1 -= norm2;
185: if (norm1<-tol || norm1>tol) {
186: PetscPrintf(PETSC_COMM_SELF,"Error:MatGetDiagonal() \n");
187: }
189: MatScale(A,alpha);
190: MatScale(sA,alpha);
192: /* Test MatMult(), MatMultAdd() */
193: for (i=0; i<40; i++) {
194: VecSetRandom(x,rdm);
195: MatMult(A,x,s1);
196: MatMult(sA,x,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: MatMult(), MatDiagonalScale() or MatScale()\n");
202: }
203: }
205: for (i=0; i<40; i++) {
206: VecSetRandom(x,rdm);
207: VecSetRandom(y,rdm);
208: MatMultAdd(A,x,y,s1);
209: MatMultAdd(sA,x,y,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:MatMultAdd(), MatDiagonalScale() or MatScale() \n");
215: }
216: }
218: /* Test MatReordering() */
219: MatGetOrdering(A,MATORDERINGNATURAL,&isrow,&iscol);
220: ip = isrow;
222: if (reorder) {
223: IS nip;
224: PetscInt *nip_ptr;
225: PetscMalloc1(mbs,&nip_ptr);
226: ISGetIndices(ip,&ip_ptr);
227: PetscArraycpy(nip_ptr,ip_ptr,mbs);
228: i = nip_ptr[1]; nip_ptr[1] = nip_ptr[mbs-2]; nip_ptr[mbs-2] = i;
229: i = nip_ptr[0]; nip_ptr[0] = nip_ptr[mbs-1]; nip_ptr[mbs-1] = i;
230: ISRestoreIndices(ip,&ip_ptr);
231: ISCreateGeneral(PETSC_COMM_SELF,mbs,nip_ptr,PETSC_COPY_VALUES,&nip);
232: PetscFree(nip_ptr);
234: MatReorderingSeqSBAIJ(sA, ip);
235: ISDestroy(&nip);
236: }
238: ISDestroy(&iscol);
239: ISDestroy(&isrow);
240: MatDestroy(&A);
241: MatDestroy(&sA);
242: VecDestroy(&x);
243: VecDestroy(&y);
244: VecDestroy(&s1);
245: VecDestroy(&s2);
246: VecDestroy(&b);
247: PetscRandomDestroy(&rdm);
249: PetscFinalize();
250: return ierr;
251: }
253: /*TEST
255: test:
256: args: -bs {{1 2 3 4 5 6 7 8}}
258: TEST*/