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