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