1 2 /* 3 Laplacian in 3D. Use for testing MatSolve routines. 4 Modeled by the partial differential equation 5 6 - Laplacian u = 1,0 < x,y,z < 1, 7 8 with boundary conditions 9 u = 1 for x = 0, x = 1, y = 0, y = 1, z = 0, z = 1. 10 */ 11 12 static char help[] = "This example is for testing different MatSolve routines :MatSolve(), MatSolveAdd(), MatSolveTranspose(), MatSolveTransposeAdd(), and MatMatSolve().\n\ 13 Example usage: ./ex129 -mat_type aij -dof 2\n\n"; 14 15 #include <petscdm.h> 16 #include <petscdmda.h> 17 18 extern PetscErrorCode ComputeMatrix(DM,Mat); 19 extern PetscErrorCode ComputeRHS(DM,Vec); 20 extern PetscErrorCode ComputeRHSMatrix(PetscInt,PetscInt,Mat*); 21 22 int main(int argc,char **args) 23 { 24 PetscMPIInt size; 25 Vec x,b,y,b1; 26 DM da; 27 Mat A,F,RHS,X,C1; 28 MatFactorInfo info; 29 IS perm,iperm; 30 PetscInt dof =1,M=8,m,n,nrhs; 31 PetscScalar one = 1.0; 32 PetscReal norm,tol = 1000*PETSC_MACHINE_EPSILON; 33 PetscBool InplaceLU=PETSC_FALSE; 34 35 PetscCall(PetscInitialize(&argc,&args,(char*)0,help)); 36 PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD,&size)); 37 PetscCheckFalse(size != 1,PETSC_COMM_WORLD,PETSC_ERR_WRONG_MPI_SIZE,"This is a uniprocessor example only"); 38 PetscCall(PetscOptionsGetInt(NULL,NULL,"-dof",&dof,NULL)); 39 PetscCall(PetscOptionsGetInt(NULL,NULL,"-M",&M,NULL)); 40 41 PetscCall(DMDACreate(PETSC_COMM_WORLD,&da)); 42 PetscCall(DMSetDimension(da,3)); 43 PetscCall(DMDASetBoundaryType(da,DM_BOUNDARY_NONE,DM_BOUNDARY_NONE,DM_BOUNDARY_NONE)); 44 PetscCall(DMDASetStencilType(da,DMDA_STENCIL_STAR)); 45 PetscCall(DMDASetSizes(da,M,M,M)); 46 PetscCall(DMDASetNumProcs(da,PETSC_DECIDE,PETSC_DECIDE,PETSC_DECIDE)); 47 PetscCall(DMDASetDof(da,dof)); 48 PetscCall(DMDASetStencilWidth(da,1)); 49 PetscCall(DMDASetOwnershipRanges(da,NULL,NULL,NULL)); 50 PetscCall(DMSetMatType(da,MATBAIJ)); 51 PetscCall(DMSetFromOptions(da)); 52 PetscCall(DMSetUp(da)); 53 54 PetscCall(DMCreateGlobalVector(da,&x)); 55 PetscCall(DMCreateGlobalVector(da,&b)); 56 PetscCall(VecDuplicate(b,&y)); 57 PetscCall(ComputeRHS(da,b)); 58 PetscCall(VecSet(y,one)); 59 PetscCall(DMCreateMatrix(da,&A)); 60 PetscCall(ComputeMatrix(da,A)); 61 PetscCall(MatGetSize(A,&m,&n)); 62 nrhs = 2; 63 PetscCall(PetscOptionsGetInt(NULL,NULL,"-nrhs",&nrhs,NULL)); 64 PetscCall(ComputeRHSMatrix(m,nrhs,&RHS)); 65 PetscCall(MatDuplicate(RHS,MAT_DO_NOT_COPY_VALUES,&X)); 66 67 PetscCall(MatGetOrdering(A,MATORDERINGND,&perm,&iperm)); 68 69 PetscCall(PetscOptionsGetBool(NULL,NULL,"-inplacelu",&InplaceLU,NULL)); 70 PetscCall(MatFactorInfoInitialize(&info)); 71 if (!InplaceLU) { 72 PetscCall(MatGetFactor(A,MATSOLVERPETSC,MAT_FACTOR_LU,&F)); 73 info.fill = 5.0; 74 PetscCall(MatLUFactorSymbolic(F,A,perm,iperm,&info)); 75 PetscCall(MatLUFactorNumeric(F,A,&info)); 76 } else { /* Test inplace factorization */ 77 PetscCall(MatDuplicate(A,MAT_COPY_VALUES,&F)); 78 PetscCall(MatLUFactor(F,perm,iperm,&info)); 79 } 80 81 PetscCall(VecDuplicate(y,&b1)); 82 83 /* MatSolve */ 84 PetscCall(MatSolve(F,b,x)); 85 PetscCall(MatMult(A,x,b1)); 86 PetscCall(VecAXPY(b1,-1.0,b)); 87 PetscCall(VecNorm(b1,NORM_2,&norm)); 88 if (norm > tol) { 89 PetscCall(PetscPrintf(PETSC_COMM_WORLD,"MatSolve : Error of norm %g\n",(double)norm)); 90 } 91 92 /* MatSolveTranspose */ 93 PetscCall(MatSolveTranspose(F,b,x)); 94 PetscCall(MatMultTranspose(A,x,b1)); 95 PetscCall(VecAXPY(b1,-1.0,b)); 96 PetscCall(VecNorm(b1,NORM_2,&norm)); 97 if (norm > tol) { 98 PetscCall(PetscPrintf(PETSC_COMM_WORLD,"MatSolveTranspose : Error of norm %g\n",(double)norm)); 99 } 100 101 /* MatSolveAdd */ 102 PetscCall(MatSolveAdd(F,b,y,x)); 103 PetscCall(MatMult(A,y,b1)); 104 PetscCall(VecScale(b1,-1.0)); 105 PetscCall(MatMultAdd(A,x,b1,b1)); 106 PetscCall(VecAXPY(b1,-1.0,b)); 107 PetscCall(VecNorm(b1,NORM_2,&norm)); 108 if (norm > tol) { 109 PetscCall(PetscPrintf(PETSC_COMM_WORLD,"MatSolveAdd : Error of norm %g\n",(double)norm)); 110 } 111 112 /* MatSolveTransposeAdd */ 113 PetscCall(MatSolveTransposeAdd(F,b,y,x)); 114 PetscCall(MatMultTranspose(A,y,b1)); 115 PetscCall(VecScale(b1,-1.0)); 116 PetscCall(MatMultTransposeAdd(A,x,b1,b1)); 117 PetscCall(VecAXPY(b1,-1.0,b)); 118 PetscCall(VecNorm(b1,NORM_2,&norm)); 119 if (norm > tol) { 120 PetscCall(PetscPrintf(PETSC_COMM_WORLD,"MatSolveTransposeAdd : Error of norm %g\n",(double)norm)); 121 } 122 123 /* MatMatSolve */ 124 PetscCall(MatMatSolve(F,RHS,X)); 125 PetscCall(MatMatMult(A,X,MAT_INITIAL_MATRIX,2.0,&C1)); 126 PetscCall(MatAXPY(C1,-1.0,RHS,SAME_NONZERO_PATTERN)); 127 PetscCall(MatNorm(C1,NORM_FROBENIUS,&norm)); 128 if (norm > tol) { 129 PetscCall(PetscPrintf(PETSC_COMM_WORLD,"MatMatSolve : Error of norm %g\n",(double)norm)); 130 } 131 132 PetscCall(VecDestroy(&x)); 133 PetscCall(VecDestroy(&b)); 134 PetscCall(VecDestroy(&b1)); 135 PetscCall(VecDestroy(&y)); 136 PetscCall(MatDestroy(&A)); 137 PetscCall(MatDestroy(&F)); 138 PetscCall(MatDestroy(&RHS)); 139 PetscCall(MatDestroy(&C1)); 140 PetscCall(MatDestroy(&X)); 141 PetscCall(ISDestroy(&perm)); 142 PetscCall(ISDestroy(&iperm)); 143 PetscCall(DMDestroy(&da)); 144 PetscCall(PetscFinalize()); 145 return 0; 146 } 147 148 PetscErrorCode ComputeRHS(DM da,Vec b) 149 { 150 PetscInt mx,my,mz; 151 PetscScalar h; 152 153 PetscFunctionBegin; 154 PetscCall(DMDAGetInfo(da,0,&mx,&my,&mz,0,0,0,0,0,0,0,0,0)); 155 h = 1.0/((mx-1)*(my-1)*(mz-1)); 156 PetscCall(VecSet(b,h)); 157 PetscFunctionReturn(0); 158 } 159 160 PetscErrorCode ComputeRHSMatrix(PetscInt m,PetscInt nrhs,Mat *C) 161 { 162 PetscRandom rand; 163 Mat RHS; 164 PetscScalar *array,rval; 165 PetscInt i,k; 166 167 PetscFunctionBegin; 168 PetscCall(MatCreate(PETSC_COMM_WORLD,&RHS)); 169 PetscCall(MatSetSizes(RHS,m,PETSC_DECIDE,PETSC_DECIDE,nrhs)); 170 PetscCall(MatSetType(RHS,MATSEQDENSE)); 171 PetscCall(MatSetUp(RHS)); 172 173 PetscCall(PetscRandomCreate(PETSC_COMM_WORLD,&rand)); 174 PetscCall(PetscRandomSetFromOptions(rand)); 175 PetscCall(MatDenseGetArray(RHS,&array)); 176 for (i=0; i<m; i++) { 177 PetscCall(PetscRandomGetValue(rand,&rval)); 178 array[i] = rval; 179 } 180 if (nrhs > 1) { 181 for (k=1; k<nrhs; k++) { 182 for (i=0; i<m; i++) { 183 array[m*k+i] = array[i]; 184 } 185 } 186 } 187 PetscCall(MatDenseRestoreArray(RHS,&array)); 188 PetscCall(MatAssemblyBegin(RHS,MAT_FINAL_ASSEMBLY)); 189 PetscCall(MatAssemblyEnd(RHS,MAT_FINAL_ASSEMBLY)); 190 *C = RHS; 191 PetscCall(PetscRandomDestroy(&rand)); 192 PetscFunctionReturn(0); 193 } 194 195 PetscErrorCode ComputeMatrix(DM da,Mat B) 196 { 197 PetscInt i,j,k,mx,my,mz,xm,ym,zm,xs,ys,zs,dof,k1,k2,k3; 198 PetscScalar *v,*v_neighbor,Hx,Hy,Hz,HxHydHz,HyHzdHx,HxHzdHy,r1,r2; 199 MatStencil row,col; 200 PetscRandom rand; 201 202 PetscFunctionBegin; 203 PetscCall(PetscRandomCreate(PETSC_COMM_WORLD,&rand)); 204 PetscCall(PetscRandomSetSeed(rand,1)); 205 PetscCall(PetscRandomSetInterval(rand,-.001,.001)); 206 PetscCall(PetscRandomSetFromOptions(rand)); 207 208 PetscCall(DMDAGetInfo(da,0,&mx,&my,&mz,0,0,0,&dof,0,0,0,0,0)); 209 /* For simplicity, this example only works on mx=my=mz */ 210 PetscCheckFalse(mx != my || mx != mz,PETSC_COMM_SELF,PETSC_ERR_SUP,"This example only works with mx %" PetscInt_FMT " = my %" PetscInt_FMT " = mz %" PetscInt_FMT,mx,my,mz); 211 212 Hx = 1.0 / (PetscReal)(mx-1); Hy = 1.0 / (PetscReal)(my-1); Hz = 1.0 / (PetscReal)(mz-1); 213 HxHydHz = Hx*Hy/Hz; HxHzdHy = Hx*Hz/Hy; HyHzdHx = Hy*Hz/Hx; 214 215 PetscCall(PetscMalloc1(2*dof*dof+1,&v)); 216 v_neighbor = v + dof*dof; 217 PetscCall(PetscArrayzero(v,2*dof*dof+1)); 218 k3 = 0; 219 for (k1=0; k1<dof; k1++) { 220 for (k2=0; k2<dof; k2++) { 221 if (k1 == k2) { 222 v[k3] = 2.0*(HxHydHz + HxHzdHy + HyHzdHx); 223 v_neighbor[k3] = -HxHydHz; 224 } else { 225 PetscCall(PetscRandomGetValue(rand,&r1)); 226 PetscCall(PetscRandomGetValue(rand,&r2)); 227 228 v[k3] = r1; 229 v_neighbor[k3] = r2; 230 } 231 k3++; 232 } 233 } 234 PetscCall(DMDAGetCorners(da,&xs,&ys,&zs,&xm,&ym,&zm)); 235 236 for (k=zs; k<zs+zm; k++) { 237 for (j=ys; j<ys+ym; j++) { 238 for (i=xs; i<xs+xm; i++) { 239 row.i = i; row.j = j; row.k = k; 240 if (i==0 || j==0 || k==0 || i==mx-1 || j==my-1 || k==mz-1) { /* boundary points */ 241 PetscCall(MatSetValuesBlockedStencil(B,1,&row,1,&row,v,INSERT_VALUES)); 242 } else { /* interior points */ 243 /* center */ 244 col.i = i; col.j = j; col.k = k; 245 PetscCall(MatSetValuesBlockedStencil(B,1,&row,1,&col,v,INSERT_VALUES)); 246 247 /* x neighbors */ 248 col.i = i-1; col.j = j; col.k = k; 249 PetscCall(MatSetValuesBlockedStencil(B,1,&row,1,&col,v_neighbor,INSERT_VALUES)); 250 col.i = i+1; col.j = j; col.k = k; 251 PetscCall(MatSetValuesBlockedStencil(B,1,&row,1,&col,v_neighbor,INSERT_VALUES)); 252 253 /* y neighbors */ 254 col.i = i; col.j = j-1; col.k = k; 255 PetscCall(MatSetValuesBlockedStencil(B,1,&row,1,&col,v_neighbor,INSERT_VALUES)); 256 col.i = i; col.j = j+1; col.k = k; 257 PetscCall(MatSetValuesBlockedStencil(B,1,&row,1,&col,v_neighbor,INSERT_VALUES)); 258 259 /* z neighbors */ 260 col.i = i; col.j = j; col.k = k-1; 261 PetscCall(MatSetValuesBlockedStencil(B,1,&row,1,&col,v_neighbor,INSERT_VALUES)); 262 col.i = i; col.j = j; col.k = k+1; 263 PetscCall(MatSetValuesBlockedStencil(B,1,&row,1,&col,v_neighbor,INSERT_VALUES)); 264 } 265 } 266 } 267 } 268 PetscCall(MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY)); 269 PetscCall(MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY)); 270 PetscCall(PetscFree(v)); 271 PetscCall(PetscRandomDestroy(&rand)); 272 PetscFunctionReturn(0); 273 } 274 275 /*TEST 276 277 test: 278 args: -dm_mat_type aij -dof 1 279 output_file: output/ex129.out 280 281 test: 282 suffix: 2 283 args: -dm_mat_type aij -dof 1 -inplacelu 284 output_file: output/ex129.out 285 286 TEST*/ 287