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 PetscErrorCode ierr; 25 PetscMPIInt size; 26 Vec x,b,y,b1; 27 DM da; 28 Mat A,F,RHS,X,C1; 29 MatFactorInfo info; 30 IS perm,iperm; 31 PetscInt dof =1,M=8,m,n,nrhs; 32 PetscScalar one = 1.0; 33 PetscReal norm,tol = 1000*PETSC_MACHINE_EPSILON; 34 PetscBool InplaceLU=PETSC_FALSE; 35 36 ierr = PetscInitialize(&argc,&args,(char*)0,help);if (ierr) return ierr; 37 CHKERRMPI(MPI_Comm_size(PETSC_COMM_WORLD,&size)); 38 PetscCheckFalse(size != 1,PETSC_COMM_WORLD,PETSC_ERR_WRONG_MPI_SIZE,"This is a uniprocessor example only"); 39 CHKERRQ(PetscOptionsGetInt(NULL,NULL,"-dof",&dof,NULL)); 40 CHKERRQ(PetscOptionsGetInt(NULL,NULL,"-M",&M,NULL)); 41 42 CHKERRQ(DMDACreate(PETSC_COMM_WORLD,&da)); 43 CHKERRQ(DMSetDimension(da,3)); 44 CHKERRQ(DMDASetBoundaryType(da,DM_BOUNDARY_NONE,DM_BOUNDARY_NONE,DM_BOUNDARY_NONE)); 45 CHKERRQ(DMDASetStencilType(da,DMDA_STENCIL_STAR)); 46 CHKERRQ(DMDASetSizes(da,M,M,M)); 47 CHKERRQ(DMDASetNumProcs(da,PETSC_DECIDE,PETSC_DECIDE,PETSC_DECIDE)); 48 CHKERRQ(DMDASetDof(da,dof)); 49 CHKERRQ(DMDASetStencilWidth(da,1)); 50 CHKERRQ(DMDASetOwnershipRanges(da,NULL,NULL,NULL)); 51 CHKERRQ(DMSetMatType(da,MATBAIJ)); 52 CHKERRQ(DMSetFromOptions(da)); 53 CHKERRQ(DMSetUp(da)); 54 55 CHKERRQ(DMCreateGlobalVector(da,&x)); 56 CHKERRQ(DMCreateGlobalVector(da,&b)); 57 CHKERRQ(VecDuplicate(b,&y)); 58 CHKERRQ(ComputeRHS(da,b)); 59 CHKERRQ(VecSet(y,one)); 60 CHKERRQ(DMCreateMatrix(da,&A)); 61 CHKERRQ(ComputeMatrix(da,A)); 62 CHKERRQ(MatGetSize(A,&m,&n)); 63 nrhs = 2; 64 CHKERRQ(PetscOptionsGetInt(NULL,NULL,"-nrhs",&nrhs,NULL)); 65 CHKERRQ(ComputeRHSMatrix(m,nrhs,&RHS)); 66 CHKERRQ(MatDuplicate(RHS,MAT_DO_NOT_COPY_VALUES,&X)); 67 68 CHKERRQ(MatGetOrdering(A,MATORDERINGND,&perm,&iperm)); 69 70 CHKERRQ(PetscOptionsGetBool(NULL,NULL,"-inplacelu",&InplaceLU,NULL)); 71 CHKERRQ(MatFactorInfoInitialize(&info)); 72 if (!InplaceLU) { 73 CHKERRQ(MatGetFactor(A,MATSOLVERPETSC,MAT_FACTOR_LU,&F)); 74 info.fill = 5.0; 75 CHKERRQ(MatLUFactorSymbolic(F,A,perm,iperm,&info)); 76 CHKERRQ(MatLUFactorNumeric(F,A,&info)); 77 } else { /* Test inplace factorization */ 78 CHKERRQ(MatDuplicate(A,MAT_COPY_VALUES,&F)); 79 CHKERRQ(MatLUFactor(F,perm,iperm,&info)); 80 } 81 82 CHKERRQ(VecDuplicate(y,&b1)); 83 84 /* MatSolve */ 85 CHKERRQ(MatSolve(F,b,x)); 86 CHKERRQ(MatMult(A,x,b1)); 87 CHKERRQ(VecAXPY(b1,-1.0,b)); 88 CHKERRQ(VecNorm(b1,NORM_2,&norm)); 89 if (norm > tol) { 90 CHKERRQ(PetscPrintf(PETSC_COMM_WORLD,"MatSolve : Error of norm %g\n",(double)norm)); 91 } 92 93 /* MatSolveTranspose */ 94 CHKERRQ(MatSolveTranspose(F,b,x)); 95 CHKERRQ(MatMultTranspose(A,x,b1)); 96 CHKERRQ(VecAXPY(b1,-1.0,b)); 97 CHKERRQ(VecNorm(b1,NORM_2,&norm)); 98 if (norm > tol) { 99 CHKERRQ(PetscPrintf(PETSC_COMM_WORLD,"MatSolveTranspose : Error of norm %g\n",(double)norm)); 100 } 101 102 /* MatSolveAdd */ 103 CHKERRQ(MatSolveAdd(F,b,y,x)); 104 CHKERRQ(MatMult(A,y,b1)); 105 CHKERRQ(VecScale(b1,-1.0)); 106 CHKERRQ(MatMultAdd(A,x,b1,b1)); 107 CHKERRQ(VecAXPY(b1,-1.0,b)); 108 CHKERRQ(VecNorm(b1,NORM_2,&norm)); 109 if (norm > tol) { 110 CHKERRQ(PetscPrintf(PETSC_COMM_WORLD,"MatSolveAdd : Error of norm %g\n",(double)norm)); 111 } 112 113 /* MatSolveTransposeAdd */ 114 CHKERRQ(MatSolveTransposeAdd(F,b,y,x)); 115 CHKERRQ(MatMultTranspose(A,y,b1)); 116 CHKERRQ(VecScale(b1,-1.0)); 117 CHKERRQ(MatMultTransposeAdd(A,x,b1,b1)); 118 CHKERRQ(VecAXPY(b1,-1.0,b)); 119 CHKERRQ(VecNorm(b1,NORM_2,&norm)); 120 if (norm > tol) { 121 CHKERRQ(PetscPrintf(PETSC_COMM_WORLD,"MatSolveTransposeAdd : Error of norm %g\n",(double)norm)); 122 } 123 124 /* MatMatSolve */ 125 CHKERRQ(MatMatSolve(F,RHS,X)); 126 CHKERRQ(MatMatMult(A,X,MAT_INITIAL_MATRIX,2.0,&C1)); 127 CHKERRQ(MatAXPY(C1,-1.0,RHS,SAME_NONZERO_PATTERN)); 128 CHKERRQ(MatNorm(C1,NORM_FROBENIUS,&norm)); 129 if (norm > tol) { 130 CHKERRQ(PetscPrintf(PETSC_COMM_WORLD,"MatMatSolve : Error of norm %g\n",(double)norm)); 131 } 132 133 CHKERRQ(VecDestroy(&x)); 134 CHKERRQ(VecDestroy(&b)); 135 CHKERRQ(VecDestroy(&b1)); 136 CHKERRQ(VecDestroy(&y)); 137 CHKERRQ(MatDestroy(&A)); 138 CHKERRQ(MatDestroy(&F)); 139 CHKERRQ(MatDestroy(&RHS)); 140 CHKERRQ(MatDestroy(&C1)); 141 CHKERRQ(MatDestroy(&X)); 142 CHKERRQ(ISDestroy(&perm)); 143 CHKERRQ(ISDestroy(&iperm)); 144 CHKERRQ(DMDestroy(&da)); 145 ierr = PetscFinalize(); 146 return ierr; 147 } 148 149 PetscErrorCode ComputeRHS(DM da,Vec b) 150 { 151 PetscInt mx,my,mz; 152 PetscScalar h; 153 154 PetscFunctionBegin; 155 CHKERRQ(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 CHKERRQ(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 CHKERRQ(MatCreate(PETSC_COMM_WORLD,&RHS)); 170 CHKERRQ(MatSetSizes(RHS,m,PETSC_DECIDE,PETSC_DECIDE,nrhs)); 171 CHKERRQ(MatSetType(RHS,MATSEQDENSE)); 172 CHKERRQ(MatSetUp(RHS)); 173 174 CHKERRQ(PetscRandomCreate(PETSC_COMM_WORLD,&rand)); 175 CHKERRQ(PetscRandomSetFromOptions(rand)); 176 CHKERRQ(MatDenseGetArray(RHS,&array)); 177 for (i=0; i<m; i++) { 178 CHKERRQ(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 CHKERRQ(MatDenseRestoreArray(RHS,&array)); 189 CHKERRQ(MatAssemblyBegin(RHS,MAT_FINAL_ASSEMBLY)); 190 CHKERRQ(MatAssemblyEnd(RHS,MAT_FINAL_ASSEMBLY)); 191 *C = RHS; 192 CHKERRQ(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 CHKERRQ(PetscRandomCreate(PETSC_COMM_WORLD,&rand)); 205 CHKERRQ(PetscRandomSetSeed(rand,1)); 206 CHKERRQ(PetscRandomSetInterval(rand,-.001,.001)); 207 CHKERRQ(PetscRandomSetFromOptions(rand)); 208 209 CHKERRQ(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 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); 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 CHKERRQ(PetscMalloc1(2*dof*dof+1,&v)); 217 v_neighbor = v + dof*dof; 218 CHKERRQ(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 CHKERRQ(PetscRandomGetValue(rand,&r1)); 227 CHKERRQ(PetscRandomGetValue(rand,&r2)); 228 229 v[k3] = r1; 230 v_neighbor[k3] = r2; 231 } 232 k3++; 233 } 234 } 235 CHKERRQ(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 CHKERRQ(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 CHKERRQ(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 CHKERRQ(MatSetValuesBlockedStencil(B,1,&row,1,&col,v_neighbor,INSERT_VALUES)); 251 col.i = i+1; col.j = j; col.k = k; 252 CHKERRQ(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 CHKERRQ(MatSetValuesBlockedStencil(B,1,&row,1,&col,v_neighbor,INSERT_VALUES)); 257 col.i = i; col.j = j+1; col.k = k; 258 CHKERRQ(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 CHKERRQ(MatSetValuesBlockedStencil(B,1,&row,1,&col,v_neighbor,INSERT_VALUES)); 263 col.i = i; col.j = j; col.k = k+1; 264 CHKERRQ(MatSetValuesBlockedStencil(B,1,&row,1,&col,v_neighbor,INSERT_VALUES)); 265 } 266 } 267 } 268 } 269 CHKERRQ(MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY)); 270 CHKERRQ(MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY)); 271 CHKERRQ(PetscFree(v)); 272 CHKERRQ(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