Actual source code: ex129.c

  2: /*
  3:   Laplacian in 3D. Use for testing MatSolve routines. 
  4:   Modeled by the partial differential equation

  6:    - Laplacian u = 1,0 < x,y,z < 1,

  8:    with boundary conditions
  9:    u = 1 for x = 0, x = 1, y = 0, y = 1, z = 0, z = 1.
 10: */

 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";

 15: #include <petscdmda.h>


 23: int main(int argc,char **args)
 24: {
 25:   PetscErrorCode    ierr;
 26:   PetscMPIInt       size;
 27:   Vec               x,b,y,b1;
 28:   DM                da;
 29:   Mat               A,F,RHS,X,C1;
 30:   MatFactorInfo     info;
 31:   IS                perm,iperm;
 32:   PetscInt          dof=1,M=-8,m,n,nrhs;
 33:   PetscScalar       one = 1.0;
 34:   PetscReal         norm;
 35:   PetscBool         InplaceLU=PETSC_FALSE;

 37:   PetscInitialize(&argc,&args,(char *)0,help);
 38:   MPI_Comm_size(PETSC_COMM_WORLD,&size);
 39:   if(size != 1) SETERRQ(PETSC_COMM_WORLD,1,"This is a uniprocessor example only\n");
 40:   PetscOptionsGetInt(PETSC_NULL,"-dof",&dof,PETSC_NULL);
 41:   PetscOptionsGetInt(PETSC_NULL,"-M",&M,PETSC_NULL);

 43:   DMDACreate(PETSC_COMM_WORLD,&da);
 44:   DMDASetDim(da,3);
 45:   DMDASetBoundaryType(da,DMDA_BOUNDARY_NONE,DMDA_BOUNDARY_NONE,DMDA_BOUNDARY_NONE);
 46:   DMDASetStencilType(da,DMDA_STENCIL_STAR);
 47:   DMDASetSizes(da,M,M,M);
 48:   DMDASetNumProcs(da,PETSC_DECIDE,PETSC_DECIDE,PETSC_DECIDE);
 49:   DMDASetDof(da,dof);
 50:   DMDASetStencilWidth(da,1);
 51:   DMDASetOwnershipRanges(da,PETSC_NULL,PETSC_NULL,PETSC_NULL);
 52:   DMSetFromOptions(da);
 53:   DMSetUp(da);

 55:   DMCreateGlobalVector(da,&x);
 56:   DMCreateGlobalVector(da,&b);
 57:   VecDuplicate(b,&y);
 58:   ComputeRHS(da,b);
 59:   VecSet(y,one);
 60:   DMGetMatrix(da,MATBAIJ,&A);
 61:   ComputeMatrix(da,A);
 62:   MatGetSize(A,&m,&n);
 63:   nrhs = 2;
 64:   PetscOptionsGetInt(PETSC_NULL,"-nrhs",&nrhs,PETSC_NULL);
 65:   ComputeRHSMatrix(m,nrhs,&RHS);
 66:   MatDuplicate(RHS,MAT_DO_NOT_COPY_VALUES,&X);

 68:   MatGetOrdering(A,MATORDERINGND,&perm,&iperm);
 69: 
 70: 
 71:   PetscOptionsGetBool(PETSC_NULL,"-inplacelu",&InplaceLU,PETSC_NULL);
 72:   MatFactorInfoInitialize(&info);
 73:   if (!InplaceLU){
 74:     MatGetFactor(A,MATSOLVERPETSC,MAT_FACTOR_LU,&F);
 75:     info.fill = 5.0;
 76:     MatLUFactorSymbolic(F,A,perm,iperm,&info);
 77:     MatLUFactorNumeric(F,A,&info);
 78:   } else { /* Test inplace factorization */
 79:     MatDuplicate(A,MAT_COPY_VALUES,&F);
 80:     /* or create F without DMDA 
 81:     const MatType     type;
 82:     PetscInt          i,ncols;
 83:     const PetscInt    *cols;
 84:     const PetscScalar *vals;
 85:     MatGetSize(A,&m,&n);
 86:     MatGetType(A,&type);
 87:     MatCreate(((PetscObject)A)->comm,&F);
 88:     MatSetSizes(F,PETSC_DECIDE,PETSC_DECIDE,m,n);
 89:     MatSetType(F,type);
 90:     MatSetFromOptions(F);
 91:     for (i=0; i<m; i++) {
 92:       MatGetRow(A,i,&ncols,&cols,&vals);
 93:       MatSetValues(F,1,&i,ncols,cols,vals,INSERT_VALUES);
 94:     }
 95:     MatAssemblyBegin(F,MAT_FINAL_ASSEMBLY);
 96:     MatAssemblyEnd(F,MAT_FINAL_ASSEMBLY);
 97:     */
 98:     MatLUFactor(F,perm,iperm,&info);
 99:   }

101:   VecDuplicate(y,&b1);
102: 
103:   /* MatSolve */
104:   MatSolve(F,b,x);
105:   MatMult(A,x,b1);
106:   VecAXPY(b1,-1.0,b);
107:   VecNorm(b1,NORM_2,&norm);
108:   PetscPrintf(PETSC_COMM_WORLD,"MatSolve              : Error of norm %A\n",norm);
109: 
110:   /* MatSolveTranspose */
111:   MatSolveTranspose(F,b,x);
112:   MatMultTranspose(A,x,b1);
113:   VecAXPY(b1,-1.0,b);
114:   VecNorm(b1,NORM_2,&norm);
115:   PetscPrintf(PETSC_COMM_WORLD,"MatSolveTranspose     : Error of norm %A\n",norm);
116: 
117:   /* MatSolveAdd */
118:   MatSolveAdd(F,b,y,x);
119:   MatMult(A,y,b1);
120:   VecScale(b1,-1.0);
121:   MatMultAdd(A,x,b1,b1);
122:   VecAXPY(b1,-1.0,b);
123:   VecNorm(b1,NORM_2,&norm);
124:   PetscPrintf(PETSC_COMM_WORLD,"MatSolveAdd           : Error of norm %A\n",norm);
125: 
126:   /* MatSolveTransposeAdd */
127:   MatSolveTransposeAdd(F,b,y,x);
128:   MatMultTranspose(A,y,b1);
129:   VecScale(b1,-1.0);
130:   MatMultTransposeAdd(A,x,b1,b1);
131:   VecAXPY(b1,-1.0,b);
132:   VecNorm(b1,NORM_2,&norm);
133:   PetscPrintf(PETSC_COMM_WORLD,"MatSolveTransposeAdd  : Error of norm %A\n",norm);
134: 
135:   /* MatMatSolve */
136:   MatMatSolve(F,RHS,X);
137:   MatMatMult(A,X,MAT_INITIAL_MATRIX,2.0,&C1);
138:   MatAXPY(C1,-1.0,RHS,SAME_NONZERO_PATTERN);
139:   MatNorm(C1,NORM_FROBENIUS,&norm);
140:   PetscPrintf(PETSC_COMM_WORLD,"MatMatSolve           : Error of norm %A\n",norm);


143:   VecDestroy(&x);
144:   VecDestroy(&b);
145:   VecDestroy(&b1);
146:   VecDestroy(&y);
147:   MatDestroy(&A);
148:   MatDestroy(&F);
149:   MatDestroy(&RHS);
150:   MatDestroy(&C1);
151:   MatDestroy(&X);
152:   ISDestroy(&perm);
153:   ISDestroy(&iperm);
154:   DMDestroy(&da);
155:   PetscFinalize();
156:   return 0;
157: }

161: PetscErrorCode ComputeRHS(DM da,Vec b)
162: {
164:   PetscInt       mx,my,mz;
165:   PetscScalar    h;

168:   DMDAGetInfo(da,0,&mx,&my,&mz,0,0,0,0,0,0,0,0,0);
169:   h    = 1.0/((mx-1)*(my-1)*(mz-1));
170:   VecSet(b,h);
171:   return(0);
172: }

176: PetscErrorCode ComputeRHSMatrix(PetscInt m,PetscInt nrhs,Mat* C)
177: {
179:   PetscRandom    rand;
180:   Mat            RHS;
181:   PetscScalar    *array,rval;
182:   PetscInt       i,k;

185:   MatCreate(PETSC_COMM_WORLD,&RHS);
186:   MatSetSizes(RHS,m,PETSC_DECIDE,PETSC_DECIDE,nrhs);
187:   MatSetType(RHS,MATSEQDENSE);
188: 
189:   PetscRandomCreate(PETSC_COMM_WORLD,&rand);
190:   PetscRandomSetFromOptions(rand);
191:   MatGetArray(RHS,&array);
192:   for (i=0; i<m; i++){
193:     PetscRandomGetValue(rand,&rval);
194:     array[i] = rval;
195:   }
196:   if (nrhs > 1){
197:     for (k=1; k<nrhs; k++){
198:       for (i=0; i<m; i++){
199:         array[m*k+i] = array[i];
200:       }
201:     }
202:   }
203:   MatRestoreArray(RHS,&array);
204:   MatAssemblyBegin(RHS,MAT_FINAL_ASSEMBLY);
205:   MatAssemblyEnd(RHS,MAT_FINAL_ASSEMBLY);
206:   *C = RHS;
207:   PetscRandomDestroy(&rand);
208:   return(0);
209: }

211: 
214: PetscErrorCode ComputeMatrix(DM da,Mat B)
215: {
217:   PetscInt       i,j,k,mx,my,mz,xm,ym,zm,xs,ys,zs,dof,k1,k2,k3;
218:   PetscScalar    *v,*v_neighbor,Hx,Hy,Hz,HxHydHz,HyHzdHx,HxHzdHy,r1,r2;
219:   MatStencil     row,col;
220:   PetscRandom    rand;

223:   PetscRandomCreate(PETSC_COMM_WORLD,&rand);
224:   PetscRandomSetType(rand,PETSCRAND);
225:   PetscRandomSetSeed(rand,1);
226:   PetscRandomSetInterval(rand,-.001,.001);
227:   PetscRandomSetFromOptions(rand);

229:   DMDAGetInfo(da,0,&mx,&my,&mz,0,0,0,&dof,0,0,0,0,0);
230:   /* For simplicity, this example only works on mx=my=mz */
231:   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);

233:   Hx = 1.0 / (PetscReal)(mx-1); Hy = 1.0 / (PetscReal)(my-1); Hz = 1.0 / (PetscReal)(mz-1);
234:   HxHydHz = Hx*Hy/Hz; HxHzdHy = Hx*Hz/Hy; HyHzdHx = Hy*Hz/Hx;

236:   PetscMalloc((2*dof*dof+1)*sizeof(PetscScalar),&v);
237:   v_neighbor = v + dof*dof;
238:   PetscMemzero(v,(2*dof*dof+1)*sizeof(PetscScalar));
239:   k3 = 0;
240:   for (k1=0; k1<dof; k1++){
241:     for (k2=0; k2<dof; k2++){
242:       if (k1 == k2){
243:         v[k3]          = 2.0*(HxHydHz + HxHzdHy + HyHzdHx);
244:         v_neighbor[k3] = -HxHydHz;
245:       } else {
246:         PetscRandomGetValue(rand,&r1);
247:         PetscRandomGetValue(rand,&r2);
248:         v[k3] = r1;
249:         v_neighbor[k3] = r2;
250:       }
251:       k3++;
252:     }
253:   }
254:   DMDAGetCorners(da,&xs,&ys,&zs,&xm,&ym,&zm);
255: 
256:   for (k=zs; k<zs+zm; k++){
257:     for (j=ys; j<ys+ym; j++){
258:       for(i=xs; i<xs+xm; i++){
259:         row.i = i; row.j = j; row.k = k;
260:         if (i==0 || j==0 || k==0 || i==mx-1 || j==my-1 || k==mz-1){ /* boudary points */
261:           MatSetValuesBlockedStencil(B,1,&row,1,&row,v,INSERT_VALUES);
262:         } else { /* interior points */
263:           /* center */
264:           col.i = i; col.j = j; col.k = k;
265:           MatSetValuesBlockedStencil(B,1,&row,1,&col,v,INSERT_VALUES);
266: 
267:           /* x neighbors */
268:           col.i = i-1; col.j = j; col.k = k;
269:           MatSetValuesBlockedStencil(B,1,&row,1,&col,v_neighbor,INSERT_VALUES);
270:           col.i = i+1; col.j = j; col.k = k;
271:           MatSetValuesBlockedStencil(B,1,&row,1,&col,v_neighbor,INSERT_VALUES);
272: 
273:           /* y neighbors */
274:           col.i = i; col.j = j-1; col.k = k;
275:           MatSetValuesBlockedStencil(B,1,&row,1,&col,v_neighbor,INSERT_VALUES);
276:           col.i = i; col.j = j+1; col.k = k;
277:           MatSetValuesBlockedStencil(B,1,&row,1,&col,v_neighbor,INSERT_VALUES);
278: 
279:           /* z neighbors */
280:           col.i = i; col.j = j; col.k = k-1;
281:           MatSetValuesBlockedStencil(B,1,&row,1,&col,v_neighbor,INSERT_VALUES);
282:           col.i = i; col.j = j; col.k = k+1;
283:           MatSetValuesBlockedStencil(B,1,&row,1,&col,v_neighbor,INSERT_VALUES);
284:         }
285:       }
286:     }
287:   }
288:   MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);
289:   MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);
290:   PetscFree(v);
291:   PetscRandomDestroy(&rand);
292:   return(0);
293: }