Actual source code: ex18.c
petsc-3.12.5 2020-03-29
1: static char help[] = "Tests the use of MatZeroRowsColumns() for parallel matrices.\n\
2: Contributed-by: Stephan Kramer <s.kramer@imperial.ac.uk>\n\n";
4: #include <petscmat.h>
6: int main(int argc,char **args)
7: {
8: Mat A;
9: Vec x, rhs, y;
10: PetscInt i,j,k,b,m = 3,n,nlocal=2,bs=1,Ii,J;
11: PetscInt *boundary_nodes, nboundary_nodes, *boundary_indices;
12: PetscMPIInt rank,size;
14: PetscScalar v,v0,v1,v2,a0=0.1,a,rhsval, *boundary_values,diag = 1.0;
15: PetscReal norm;
16: char convname[64];
17: PetscBool upwind = PETSC_FALSE, nonlocalBC = PETSC_FALSE, zerorhs = PETSC_TRUE, convert = PETSC_FALSE;
19: PetscInitialize(&argc,&args,(char*)0,help);if (ierr) return ierr;
20: MPI_Comm_rank(PETSC_COMM_WORLD,&rank);
21: MPI_Comm_size(PETSC_COMM_WORLD,&size);
22: n = nlocal*size;
24: PetscOptionsGetInt(NULL,NULL, "-bs", &bs, NULL);
25: PetscOptionsGetBool(NULL,NULL, "-nonlocal_bc", &nonlocalBC, NULL);
26: PetscOptionsGetScalar(NULL,NULL, "-diag", &diag, NULL);
27: PetscOptionsGetString(NULL,NULL,"-convname",convname,sizeof(convname),&convert);
28: PetscOptionsGetBool(NULL,NULL, "-zerorhs", &zerorhs, NULL);
30: MatCreate(PETSC_COMM_WORLD,&A);
31: MatSetSizes(A,PETSC_DECIDE,PETSC_DECIDE,m*n*bs,m*n*bs);
32: MatSetFromOptions(A);
33: MatSetUp(A);
35: VecCreate(PETSC_COMM_WORLD, &rhs);
36: VecSetSizes(rhs, PETSC_DECIDE, m*n*bs);
37: VecSetFromOptions(rhs);
38: VecSetUp(rhs);
40: rhsval = 0.0;
41: for (i=0; i<m; i++) {
42: for (j=nlocal*rank; j<nlocal*(rank+1); j++) {
43: a = a0;
44: for (b=0; b<bs; b++) {
45: /* let's start with a 5-point stencil diffusion term */
46: v = -1.0; Ii = (j + n*i)*bs + b;
47: if (i>0) {J = Ii - n*bs; MatSetValues(A,1,&Ii,1,&J,&v,ADD_VALUES);}
48: if (i<m-1) {J = Ii + n*bs; MatSetValues(A,1,&Ii,1,&J,&v,ADD_VALUES);}
49: if (j>0) {J = Ii - 1*bs; MatSetValues(A,1,&Ii,1,&J,&v,ADD_VALUES);}
50: if (j<n-1) {J = Ii + 1*bs; MatSetValues(A,1,&Ii,1,&J,&v,ADD_VALUES);}
51: v = 4.0; MatSetValues(A,1,&Ii,1,&Ii,&v,ADD_VALUES);
52: if (upwind) {
53: /* now add a 2nd order upwind advection term to add a little asymmetry */
54: if (j>2) {
55: J = Ii-2*bs; v2 = 0.5*a; v1 = -2.0*a; v0 = 1.5*a;
56: MatSetValues(A,1,&Ii,1,&J,&v2,ADD_VALUES);
57: } else {
58: /* fall back to 1st order upwind */
59: v1 = -1.0*a; v0 = 1.0*a;
60: };
61: if (j>1) {J = Ii-1*bs; MatSetValues(A,1,&Ii,1,&J,&v1,ADD_VALUES);}
62: MatSetValues(A,1,&Ii,1,&Ii,&v0,ADD_VALUES);
63: a /= 10.; /* use a different velocity for the next component */
64: /* add a coupling to the previous and next components */
65: v = 0.5;
66: if (b>0) {J = Ii - 1; MatSetValues(A,1,&Ii,1,&J,&v,ADD_VALUES);}
67: if (b<bs-1) {J = Ii + 1; MatSetValues(A,1,&Ii,1,&J,&v,ADD_VALUES);}
68: }
69: /* make up some rhs */
70: VecSetValue(rhs, Ii, rhsval, INSERT_VALUES);
71: rhsval += 1.0;
72: }
73: }
74: }
75: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
76: MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
78: if (convert) { /* Test different Mat implementations */
79: Mat B;
81: MatConvert(A,convname,MAT_INITIAL_MATRIX,&B);
82: MatDestroy(&A);
83: A = B;
84: }
86: VecAssemblyBegin(rhs);
87: VecAssemblyEnd(rhs);
88: /* set rhs to zero to simplify */
89: if (zerorhs) {
90: VecZeroEntries(rhs);
91: }
93: if (nonlocalBC) {
94: /*version where boundary conditions are set by processes that don't necessarily own the nodes */
95: if (!rank) {
96: nboundary_nodes = size>m ? nlocal : m-size+nlocal;
97: PetscMalloc1(nboundary_nodes,&boundary_nodes);
98: k = 0;
99: for (i=size; i<m; i++,k++) {boundary_nodes[k] = n*i;};
100: } else if (rank < m) {
101: nboundary_nodes = nlocal+1;
102: PetscMalloc1(nboundary_nodes,&boundary_nodes);
103: boundary_nodes[0] = rank*n;
104: k = 1;
105: } else {
106: nboundary_nodes = nlocal;
107: PetscMalloc1(nboundary_nodes,&boundary_nodes);
108: k = 0;
109: }
110: for (j=nlocal*rank; j<nlocal*(rank+1); j++,k++) {boundary_nodes[k] = j;};
111: } else {
112: /*version where boundary conditions are set by the node owners only */
113: PetscMalloc1(m*n,&boundary_nodes);
114: k=0;
115: for (j=0; j<n; j++) {
116: Ii = j;
117: if (Ii>=rank*m*nlocal && Ii<(rank+1)*m*nlocal) boundary_nodes[k++] = Ii;
118: }
119: for (i=1; i<m; i++) {
120: Ii = n*i;
121: if (Ii>=rank*m*nlocal && Ii<(rank+1)*m*nlocal) boundary_nodes[k++] = Ii;
122: }
123: nboundary_nodes = k;
124: }
126: VecDuplicate(rhs, &x);
127: VecZeroEntries(x);
128: PetscMalloc2(nboundary_nodes*bs,&boundary_indices,nboundary_nodes*bs,&boundary_values);
129: for (k=0; k<nboundary_nodes; k++) {
130: Ii = boundary_nodes[k]*bs;
131: v = 1.0*boundary_nodes[k];
132: for (b=0; b<bs; b++, Ii++) {
133: boundary_indices[k*bs+b] = Ii;
134: boundary_values[k*bs+b] = v;
135: PetscSynchronizedPrintf(PETSC_COMM_WORLD, "%d %D %f\n", rank, Ii, (double)PetscRealPart(v));
136: v += 0.1;
137: }
138: }
139: PetscSynchronizedFlush(PETSC_COMM_WORLD, NULL);
140: VecSetValues(x, nboundary_nodes*bs, boundary_indices, boundary_values, INSERT_VALUES);
141: VecAssemblyBegin(x);
142: VecAssemblyEnd(x);
144: /* We can check the rhs returned by MatZeroColumns by computing y=rhs-A*x and overwriting the boundary entries with boundary values */
145: VecDuplicate(x, &y);
146: MatMult(A, x, y);
147: VecAYPX(y, -1.0, rhs);
148: for (k=0; k<nboundary_nodes*bs; k++) boundary_values[k] *= diag;
149: VecSetValues(y, nboundary_nodes*bs, boundary_indices, boundary_values, INSERT_VALUES);
150: VecAssemblyBegin(y);
151: VecAssemblyEnd(y);
153: PetscPrintf(PETSC_COMM_WORLD, "*** Matrix A and vector x:\n");
154: MatView(A, PETSC_VIEWER_STDOUT_WORLD);
155: VecView(x,PETSC_VIEWER_STDOUT_WORLD);
157: MatZeroRowsColumns(A, nboundary_nodes*bs, boundary_indices, diag, x, rhs);
158: PetscPrintf(PETSC_COMM_WORLD, "*** Vector rhs returned by MatZeroRowsColumns\n");
159: VecView(rhs,PETSC_VIEWER_STDOUT_WORLD);
160: VecAXPY(y, -1.0, rhs);
161: VecNorm(y, NORM_INFINITY, &norm);
162: if (norm > 1.0e-10) {
163: PetscPrintf(PETSC_COMM_WORLD, "*** Difference between rhs and y, inf-norm: %f\n", (double)norm);
164: VecView(y,PETSC_VIEWER_STDOUT_WORLD);
165: SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_PLIB, "Bug in MatZeroRowsColumns");
166: }
168: PetscFree(boundary_nodes);
169: PetscFree2(boundary_indices,boundary_values);
170: VecDestroy(&x);
171: VecDestroy(&y);
172: VecDestroy(&rhs);
173: MatDestroy(&A);
175: PetscFinalize();
176: return ierr;
177: }
180: /*TEST
182: test:
183: suffix: 0
185: test:
186: suffix: 1
187: nsize: 2
189: test:
190: suffix: 10
191: nsize: 2
192: args: -bs 2 -nonlocal_bc
194: test:
195: suffix: 11
196: nsize: 7
197: args: -bs 2 -nonlocal_bc
199: test:
200: suffix: 12
201: args: -bs 2 -nonlocal_bc -mat_type baij
203: test:
204: suffix: 13
205: nsize: 2
206: args: -bs 2 -nonlocal_bc -mat_type baij
208: test:
209: suffix: 14
210: nsize: 7
211: args: -bs 2 -nonlocal_bc -mat_type baij
213: test:
214: suffix: 2
215: nsize: 7
217: test:
218: suffix: 3
219: args: -mat_type baij
221: test:
222: suffix: 4
223: nsize: 2
224: args: -mat_type baij
226: test:
227: suffix: 5
228: nsize: 7
229: args: -mat_type baij
231: test:
232: suffix: 6
233: args: -bs 2 -mat_type baij
235: test:
236: suffix: 7
237: nsize: 2
238: args: -bs 2 -mat_type baij
240: test:
241: suffix: 8
242: nsize: 7
243: args: -bs 2 -mat_type baij
245: test:
246: suffix: 9
247: args: -bs 2 -nonlocal_bc
249: test:
250: suffix: 15
251: args: -bs 2 -nonlocal_bc -convname shell
253: test:
254: suffix: 16
255: nsize: 2
256: args: -bs 2 -nonlocal_bc -convname shell
258: test:
259: suffix: 17
260: args: -bs 2 -nonlocal_bc -convname dense
262: testset:
263: suffix: full
264: nsize: {{1 3}separate output}
265: args: -diag {{0.12 -0.13}separate output} -convname {{aij shell baij}separate output} -zerorhs 0
266: TEST*/