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"; 3 4 #include <petscmat.h> 5 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; 13 PetscErrorCode ierr; 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; 18 19 ierr = PetscInitialize(&argc,&args,(char*)0,help);if (ierr) return ierr; 20 ierr = MPI_Comm_rank(PETSC_COMM_WORLD,&rank);CHKERRQ(ierr); 21 ierr = MPI_Comm_size(PETSC_COMM_WORLD,&size);CHKERRQ(ierr); 22 n = nlocal*size; 23 24 ierr = PetscOptionsGetInt(NULL,NULL, "-bs", &bs, NULL);CHKERRQ(ierr); 25 ierr = PetscOptionsGetBool(NULL,NULL, "-nonlocal_bc", &nonlocalBC, NULL);CHKERRQ(ierr); 26 ierr = PetscOptionsGetScalar(NULL,NULL, "-diag", &diag, NULL);CHKERRQ(ierr); 27 ierr = PetscOptionsGetString(NULL,NULL,"-convname",convname,sizeof(convname),&convert);CHKERRQ(ierr); 28 ierr = PetscOptionsGetBool(NULL,NULL, "-zerorhs", &zerorhs, NULL);CHKERRQ(ierr); 29 30 ierr = MatCreate(PETSC_COMM_WORLD,&A);CHKERRQ(ierr); 31 ierr = MatSetSizes(A,PETSC_DECIDE,PETSC_DECIDE,m*n*bs,m*n*bs);CHKERRQ(ierr); 32 ierr = MatSetFromOptions(A);CHKERRQ(ierr); 33 ierr = MatSetUp(A);CHKERRQ(ierr); 34 35 ierr = VecCreate(PETSC_COMM_WORLD, &rhs);CHKERRQ(ierr); 36 ierr = VecSetSizes(rhs, PETSC_DECIDE, m*n*bs);CHKERRQ(ierr); 37 ierr = VecSetFromOptions(rhs);CHKERRQ(ierr); 38 ierr = VecSetUp(rhs);CHKERRQ(ierr); 39 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; ierr = MatSetValues(A,1,&Ii,1,&J,&v,ADD_VALUES);CHKERRQ(ierr);} 48 if (i<m-1) {J = Ii + n*bs; ierr = MatSetValues(A,1,&Ii,1,&J,&v,ADD_VALUES);CHKERRQ(ierr);} 49 if (j>0) {J = Ii - 1*bs; ierr = MatSetValues(A,1,&Ii,1,&J,&v,ADD_VALUES);CHKERRQ(ierr);} 50 if (j<n-1) {J = Ii + 1*bs; ierr = MatSetValues(A,1,&Ii,1,&J,&v,ADD_VALUES);CHKERRQ(ierr);} 51 v = 4.0; ierr = MatSetValues(A,1,&Ii,1,&Ii,&v,ADD_VALUES);CHKERRQ(ierr); 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 ierr = MatSetValues(A,1,&Ii,1,&J,&v2,ADD_VALUES);CHKERRQ(ierr); 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; ierr = MatSetValues(A,1,&Ii,1,&J,&v1,ADD_VALUES);CHKERRQ(ierr);} 62 ierr = MatSetValues(A,1,&Ii,1,&Ii,&v0,ADD_VALUES);CHKERRQ(ierr); 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; ierr = MatSetValues(A,1,&Ii,1,&J,&v,ADD_VALUES);CHKERRQ(ierr);} 67 if (b<bs-1) {J = Ii + 1; ierr = MatSetValues(A,1,&Ii,1,&J,&v,ADD_VALUES);CHKERRQ(ierr);} 68 } 69 /* make up some rhs */ 70 ierr = VecSetValue(rhs, Ii, rhsval, INSERT_VALUES);CHKERRQ(ierr); 71 rhsval += 1.0; 72 } 73 } 74 } 75 ierr = MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 76 ierr = MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 77 78 if (convert) { /* Test different Mat implementations */ 79 Mat B; 80 81 ierr = MatConvert(A,convname,MAT_INITIAL_MATRIX,&B);CHKERRQ(ierr); 82 ierr = MatDestroy(&A);CHKERRQ(ierr); 83 A = B; 84 } 85 86 ierr = VecAssemblyBegin(rhs);CHKERRQ(ierr); 87 ierr = VecAssemblyEnd(rhs);CHKERRQ(ierr); 88 /* set rhs to zero to simplify */ 89 if (zerorhs) { 90 ierr = VecZeroEntries(rhs);CHKERRQ(ierr); 91 } 92 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 ierr = PetscMalloc1(nboundary_nodes,&boundary_nodes);CHKERRQ(ierr); 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 ierr = PetscMalloc1(nboundary_nodes,&boundary_nodes);CHKERRQ(ierr); 103 boundary_nodes[0] = rank*n; 104 k = 1; 105 } else { 106 nboundary_nodes = nlocal; 107 ierr = PetscMalloc1(nboundary_nodes,&boundary_nodes);CHKERRQ(ierr); 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 ierr = PetscMalloc1(m*n,&boundary_nodes);CHKERRQ(ierr); 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 } 125 126 ierr = VecDuplicate(rhs, &x);CHKERRQ(ierr); 127 ierr = VecZeroEntries(x);CHKERRQ(ierr); 128 ierr = PetscMalloc2(nboundary_nodes*bs,&boundary_indices,nboundary_nodes*bs,&boundary_values);CHKERRQ(ierr); 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 ierr = PetscSynchronizedPrintf(PETSC_COMM_WORLD, "%d %D %f\n", rank, Ii, (double)PetscRealPart(v));CHKERRQ(ierr); 136 v += 0.1; 137 } 138 } 139 ierr = PetscSynchronizedFlush(PETSC_COMM_WORLD, NULL);CHKERRQ(ierr); 140 ierr = VecSetValues(x, nboundary_nodes*bs, boundary_indices, boundary_values, INSERT_VALUES);CHKERRQ(ierr); 141 ierr = VecAssemblyBegin(x);CHKERRQ(ierr); 142 ierr = VecAssemblyEnd(x);CHKERRQ(ierr); 143 144 /* We can check the rhs returned by MatZeroColumns by computing y=rhs-A*x and overwriting the boundary entries with boundary values */ 145 ierr = VecDuplicate(x, &y);CHKERRQ(ierr); 146 ierr = MatMult(A, x, y);CHKERRQ(ierr); 147 ierr = VecAYPX(y, -1.0, rhs);CHKERRQ(ierr); 148 for (k=0; k<nboundary_nodes*bs; k++) boundary_values[k] *= diag; 149 ierr = VecSetValues(y, nboundary_nodes*bs, boundary_indices, boundary_values, INSERT_VALUES);CHKERRQ(ierr); 150 ierr = VecAssemblyBegin(y);CHKERRQ(ierr); 151 ierr = VecAssemblyEnd(y);CHKERRQ(ierr); 152 153 ierr = PetscPrintf(PETSC_COMM_WORLD, "*** Matrix A and vector x:\n");CHKERRQ(ierr); 154 ierr = MatView(A, PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr); 155 ierr = VecView(x,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr); 156 157 ierr = MatZeroRowsColumns(A, nboundary_nodes*bs, boundary_indices, diag, x, rhs);CHKERRQ(ierr); 158 ierr = PetscPrintf(PETSC_COMM_WORLD, "*** Vector rhs returned by MatZeroRowsColumns\n");CHKERRQ(ierr); 159 ierr = VecView(rhs,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr); 160 ierr = VecAXPY(y, -1.0, rhs);CHKERRQ(ierr); 161 ierr = VecNorm(y, NORM_INFINITY, &norm);CHKERRQ(ierr); 162 if (norm > 1.0e-10) { 163 ierr = PetscPrintf(PETSC_COMM_WORLD, "*** Difference between rhs and y, inf-norm: %f\n", (double)norm);CHKERRQ(ierr); 164 ierr = VecView(y,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr); 165 SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_PLIB, "Bug in MatZeroRowsColumns"); 166 } 167 168 ierr = PetscFree(boundary_nodes);CHKERRQ(ierr); 169 ierr = PetscFree2(boundary_indices,boundary_values);CHKERRQ(ierr); 170 ierr = VecDestroy(&x);CHKERRQ(ierr); 171 ierr = VecDestroy(&y);CHKERRQ(ierr); 172 ierr = VecDestroy(&rhs);CHKERRQ(ierr); 173 ierr = MatDestroy(&A);CHKERRQ(ierr); 174 175 ierr = PetscFinalize(); 176 return ierr; 177 } 178 179 180 /*TEST 181 182 test: 183 suffix: 0 184 185 test: 186 suffix: 1 187 nsize: 2 188 189 test: 190 suffix: 10 191 nsize: 2 192 args: -bs 2 -nonlocal_bc 193 194 test: 195 suffix: 11 196 nsize: 7 197 args: -bs 2 -nonlocal_bc 198 199 test: 200 suffix: 12 201 args: -bs 2 -nonlocal_bc -mat_type baij 202 203 test: 204 suffix: 13 205 nsize: 2 206 args: -bs 2 -nonlocal_bc -mat_type baij 207 208 test: 209 suffix: 14 210 nsize: 7 211 args: -bs 2 -nonlocal_bc -mat_type baij 212 213 test: 214 suffix: 2 215 nsize: 7 216 217 test: 218 suffix: 3 219 args: -mat_type baij 220 221 test: 222 suffix: 4 223 nsize: 2 224 args: -mat_type baij 225 226 test: 227 suffix: 5 228 nsize: 7 229 args: -mat_type baij 230 231 test: 232 suffix: 6 233 args: -bs 2 -mat_type baij 234 235 test: 236 suffix: 7 237 nsize: 2 238 args: -bs 2 -mat_type baij 239 240 test: 241 suffix: 8 242 nsize: 7 243 args: -bs 2 -mat_type baij 244 245 test: 246 suffix: 9 247 args: -bs 2 -nonlocal_bc 248 249 test: 250 suffix: 15 251 args: -bs 2 -nonlocal_bc -convname shell 252 253 test: 254 suffix: 16 255 nsize: 2 256 args: -bs 2 -nonlocal_bc -convname shell 257 258 test: 259 suffix: 17 260 args: -bs 2 -nonlocal_bc -convname dense 261 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*/ 267