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