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 PetscScalar v,v0,v1,v2,a0=0.1,a,rhsval, *boundary_values,diag = 1.0; 14 PetscReal norm; 15 char convname[64]; 16 PetscBool upwind = PETSC_FALSE, nonlocalBC = PETSC_FALSE, zerorhs = PETSC_TRUE, convert = PETSC_FALSE; 17 18 PetscCall(PetscInitialize(&argc,&args,(char*)0,help)); 19 PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD,&rank)); 20 PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD,&size)); 21 n = nlocal*size; 22 23 PetscCall(PetscOptionsGetInt(NULL,NULL, "-bs", &bs, NULL)); 24 PetscCall(PetscOptionsGetBool(NULL,NULL, "-nonlocal_bc", &nonlocalBC, NULL)); 25 PetscCall(PetscOptionsGetScalar(NULL,NULL, "-diag", &diag, NULL)); 26 PetscCall(PetscOptionsGetString(NULL,NULL,"-convname",convname,sizeof(convname),&convert)); 27 PetscCall(PetscOptionsGetBool(NULL,NULL, "-zerorhs", &zerorhs, NULL)); 28 29 PetscCall(MatCreate(PETSC_COMM_WORLD,&A)); 30 PetscCall(MatSetSizes(A,PETSC_DECIDE,PETSC_DECIDE,m*n*bs,m*n*bs)); 31 PetscCall(MatSetFromOptions(A)); 32 PetscCall(MatSetUp(A)); 33 34 PetscCall(MatCreateVecs(A, NULL, &rhs)); 35 PetscCall(VecSetFromOptions(rhs)); 36 PetscCall(VecSetUp(rhs)); 37 38 rhsval = 0.0; 39 for (i=0; i<m; i++) { 40 for (j=nlocal*rank; j<nlocal*(rank+1); j++) { 41 a = a0; 42 for (b=0; b<bs; b++) { 43 /* let's start with a 5-point stencil diffusion term */ 44 v = -1.0; Ii = (j + n*i)*bs + b; 45 if (i>0) {J = Ii - n*bs; PetscCall(MatSetValues(A,1,&Ii,1,&J,&v,ADD_VALUES));} 46 if (i<m-1) {J = Ii + n*bs; PetscCall(MatSetValues(A,1,&Ii,1,&J,&v,ADD_VALUES));} 47 if (j>0) {J = Ii - 1*bs; PetscCall(MatSetValues(A,1,&Ii,1,&J,&v,ADD_VALUES));} 48 if (j<n-1) {J = Ii + 1*bs; PetscCall(MatSetValues(A,1,&Ii,1,&J,&v,ADD_VALUES));} 49 v = 4.0; PetscCall(MatSetValues(A,1,&Ii,1,&Ii,&v,ADD_VALUES)); 50 if (upwind) { 51 /* now add a 2nd order upwind advection term to add a little asymmetry */ 52 if (j>2) { 53 J = Ii-2*bs; v2 = 0.5*a; v1 = -2.0*a; v0 = 1.5*a; 54 PetscCall(MatSetValues(A,1,&Ii,1,&J,&v2,ADD_VALUES)); 55 } else { 56 /* fall back to 1st order upwind */ 57 v1 = -1.0*a; v0 = 1.0*a; 58 }; 59 if (j>1) {J = Ii-1*bs; PetscCall(MatSetValues(A,1,&Ii,1,&J,&v1,ADD_VALUES));} 60 PetscCall(MatSetValues(A,1,&Ii,1,&Ii,&v0,ADD_VALUES)); 61 a /= 10.; /* use a different velocity for the next component */ 62 /* add a coupling to the previous and next components */ 63 v = 0.5; 64 if (b>0) {J = Ii - 1; PetscCall(MatSetValues(A,1,&Ii,1,&J,&v,ADD_VALUES));} 65 if (b<bs-1) {J = Ii + 1; PetscCall(MatSetValues(A,1,&Ii,1,&J,&v,ADD_VALUES));} 66 } 67 /* make up some rhs */ 68 PetscCall(VecSetValue(rhs, Ii, rhsval, INSERT_VALUES)); 69 rhsval += 1.0; 70 } 71 } 72 } 73 PetscCall(MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY)); 74 PetscCall(MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY)); 75 76 if (convert) { /* Test different Mat implementations */ 77 Mat B; 78 79 PetscCall(MatConvert(A,convname,MAT_INITIAL_MATRIX,&B)); 80 PetscCall(MatDestroy(&A)); 81 A = B; 82 } 83 84 PetscCall(VecAssemblyBegin(rhs)); 85 PetscCall(VecAssemblyEnd(rhs)); 86 /* set rhs to zero to simplify */ 87 if (zerorhs) PetscCall(VecZeroEntries(rhs)); 88 89 if (nonlocalBC) { 90 /*version where boundary conditions are set by processes that don't necessarily own the nodes */ 91 if (rank == 0) { 92 nboundary_nodes = size>m ? nlocal : m-size+nlocal; 93 PetscCall(PetscMalloc1(nboundary_nodes,&boundary_nodes)); 94 k = 0; 95 for (i=size; i<m; i++,k++) {boundary_nodes[k] = n*i;}; 96 } else if (rank < m) { 97 nboundary_nodes = nlocal+1; 98 PetscCall(PetscMalloc1(nboundary_nodes,&boundary_nodes)); 99 boundary_nodes[0] = rank*n; 100 k = 1; 101 } else { 102 nboundary_nodes = nlocal; 103 PetscCall(PetscMalloc1(nboundary_nodes,&boundary_nodes)); 104 k = 0; 105 } 106 for (j=nlocal*rank; j<nlocal*(rank+1); j++,k++) {boundary_nodes[k] = j;}; 107 } else { 108 /*version where boundary conditions are set by the node owners only */ 109 PetscCall(PetscMalloc1(m*n,&boundary_nodes)); 110 k=0; 111 for (j=0; j<n; j++) { 112 Ii = j; 113 if (Ii>=rank*m*nlocal && Ii<(rank+1)*m*nlocal) boundary_nodes[k++] = Ii; 114 } 115 for (i=1; i<m; i++) { 116 Ii = n*i; 117 if (Ii>=rank*m*nlocal && Ii<(rank+1)*m*nlocal) boundary_nodes[k++] = Ii; 118 } 119 nboundary_nodes = k; 120 } 121 122 PetscCall(VecDuplicate(rhs, &x)); 123 PetscCall(VecZeroEntries(x)); 124 PetscCall(PetscMalloc2(nboundary_nodes*bs,&boundary_indices,nboundary_nodes*bs,&boundary_values)); 125 for (k=0; k<nboundary_nodes; k++) { 126 Ii = boundary_nodes[k]*bs; 127 v = 1.0*boundary_nodes[k]; 128 for (b=0; b<bs; b++, Ii++) { 129 boundary_indices[k*bs+b] = Ii; 130 boundary_values[k*bs+b] = v; 131 PetscCall(PetscSynchronizedPrintf(PETSC_COMM_WORLD, "%d %" PetscInt_FMT " %f\n", rank, Ii, (double)PetscRealPart(v))); 132 v += 0.1; 133 } 134 } 135 PetscCall(PetscSynchronizedFlush(PETSC_COMM_WORLD, NULL)); 136 PetscCall(VecSetValues(x, nboundary_nodes*bs, boundary_indices, boundary_values, INSERT_VALUES)); 137 PetscCall(VecAssemblyBegin(x)); 138 PetscCall(VecAssemblyEnd(x)); 139 140 /* We can check the rhs returned by MatZeroColumns by computing y=rhs-A*x and overwriting the boundary entries with boundary values */ 141 PetscCall(VecDuplicate(x, &y)); 142 PetscCall(MatMult(A, x, y)); 143 PetscCall(VecAYPX(y, -1.0, rhs)); 144 for (k=0; k<nboundary_nodes*bs; k++) boundary_values[k] *= diag; 145 PetscCall(VecSetValues(y, nboundary_nodes*bs, boundary_indices, boundary_values, INSERT_VALUES)); 146 PetscCall(VecAssemblyBegin(y)); 147 PetscCall(VecAssemblyEnd(y)); 148 149 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "*** Matrix A and vector x:\n")); 150 PetscCall(MatView(A, PETSC_VIEWER_STDOUT_WORLD)); 151 PetscCall(VecView(x,PETSC_VIEWER_STDOUT_WORLD)); 152 153 PetscCall(MatZeroRowsColumns(A, nboundary_nodes*bs, boundary_indices, diag, x, rhs)); 154 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "*** Vector rhs returned by MatZeroRowsColumns\n")); 155 PetscCall(VecView(rhs,PETSC_VIEWER_STDOUT_WORLD)); 156 PetscCall(VecAXPY(y, -1.0, rhs)); 157 PetscCall(VecNorm(y, NORM_INFINITY, &norm)); 158 if (norm > 1.0e-10) { 159 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "*** Difference between rhs and y, inf-norm: %f\n", (double)norm)); 160 PetscCall(VecView(y,PETSC_VIEWER_STDOUT_WORLD)); 161 SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_PLIB, "Bug in MatZeroRowsColumns"); 162 } 163 164 PetscCall(PetscFree(boundary_nodes)); 165 PetscCall(PetscFree2(boundary_indices,boundary_values)); 166 PetscCall(VecDestroy(&x)); 167 PetscCall(VecDestroy(&y)); 168 PetscCall(VecDestroy(&rhs)); 169 PetscCall(MatDestroy(&A)); 170 171 PetscCall(PetscFinalize()); 172 return 0; 173 } 174 175 /*TEST 176 177 test: 178 diff_args: -j 179 suffix: 0 180 181 test: 182 diff_args: -j 183 suffix: 1 184 nsize: 2 185 186 test: 187 diff_args: -j 188 suffix: 10 189 nsize: 2 190 args: -bs 2 -nonlocal_bc 191 192 test: 193 diff_args: -j 194 suffix: 11 195 nsize: 7 196 args: -bs 2 -nonlocal_bc 197 198 test: 199 diff_args: -j 200 suffix: 12 201 args: -bs 2 -nonlocal_bc -mat_type baij 202 203 test: 204 diff_args: -j 205 suffix: 13 206 nsize: 2 207 args: -bs 2 -nonlocal_bc -mat_type baij 208 209 test: 210 diff_args: -j 211 suffix: 14 212 nsize: 7 213 args: -bs 2 -nonlocal_bc -mat_type baij 214 215 test: 216 diff_args: -j 217 suffix: 2 218 nsize: 7 219 220 test: 221 diff_args: -j 222 suffix: 3 223 args: -mat_type baij 224 225 test: 226 diff_args: -j 227 suffix: 4 228 nsize: 2 229 args: -mat_type baij 230 231 test: 232 diff_args: -j 233 suffix: 5 234 nsize: 7 235 args: -mat_type baij 236 237 test: 238 diff_args: -j 239 suffix: 6 240 args: -bs 2 -mat_type baij 241 242 test: 243 diff_args: -j 244 suffix: 7 245 nsize: 2 246 args: -bs 2 -mat_type baij 247 248 test: 249 diff_args: -j 250 suffix: 8 251 nsize: 7 252 args: -bs 2 -mat_type baij 253 254 test: 255 diff_args: -j 256 suffix: 9 257 args: -bs 2 -nonlocal_bc 258 259 test: 260 diff_args: -j 261 suffix: 15 262 args: -bs 2 -nonlocal_bc -convname shell 263 264 test: 265 diff_args: -j 266 suffix: 16 267 nsize: 2 268 args: -bs 2 -nonlocal_bc -convname shell 269 270 test: 271 diff_args: -j 272 suffix: 17 273 args: -bs 2 -nonlocal_bc -convname dense 274 275 testset: 276 diff_args: -j 277 suffix: full 278 nsize: {{1 3}separate output} 279 args: -diag {{0.12 -0.13}separate output} -convname {{aij shell baij}separate output} -zerorhs 0 280 281 test: 282 diff_args: -j 283 requires: cuda 284 suffix: cusparse_1 285 nsize: 1 286 args: -diag {{0.12 -0.13}separate output} -convname {{seqaijcusparse mpiaijcusparse}separate output} -zerorhs 0 -mat_type {{seqaijcusparse mpiaijcusparse}separate output} 287 288 test: 289 diff_args: -j 290 requires: cuda 291 suffix: cusparse_3 292 nsize: 3 293 args: -diag {{0.12 -0.13}separate output} -convname mpiaijcusparse -zerorhs 0 -mat_type mpiaijcusparse 294 295 TEST*/ 296