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 petscda.h
 16:  #include petscmg.h


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

 37:   PetscInitialize(&argc,&args,(char *)0,help);
 38:   MPI_Comm_size(PETSC_COMM_WORLD,&size);
 39:   if(size != 1) SETERRQ(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:   DACreate(PETSC_COMM_WORLD,&da);
 44:   DASetDim(da,3);
 45:   DASetPeriodicity(da,DA_NONPERIODIC);
 46:   DASetStencilType(da,DA_STENCIL_STAR);
 47:   DASetSizes(da,M,M,M);
 48:   DASetNumProcs(da,PETSC_DECIDE,PETSC_DECIDE,PETSC_DECIDE);
 49:   DASetDof(da,dof);
 50:   DASetStencilWidth(da,1);
 51:   DASetVertexDivision(da,PETSC_NULL,PETSC_NULL,PETSC_NULL);
 52:   DASetFromOptions(da);

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

 68:   MatGetOrdering(A,MATORDERING_ND,&perm,&iperm);
 69:   MatGetFactor(A,MAT_SOLVER_PETSC,MAT_FACTOR_LU,&F);
 70:   MatFactorInfoInitialize(&info);
 71:   info.fill = 5.0;
 72:   MatLUFactorSymbolic(F,A,perm,iperm,&info);
 73:   MatLUFactorNumeric(F,A,&info);
 74:   VecDuplicate(y,&b1);
 75: 
 76:   /* MatSolve */
 77:   MatSolve(F,b,x);
 78:   MatMult(A,x,b1);
 79:   VecAXPY(b1,-1.0,b);
 80:   VecNorm(b1,NORM_2,&norm);
 81:   PetscPrintf(PETSC_COMM_WORLD,"MatSolve              : Error of norm %A\n",norm);
 82: 
 83:   /* MatSolveTranspose */
 84:   MatSolveTranspose(F,b,x);
 85:   MatMultTranspose(A,x,b1);
 86:   VecAXPY(b1,-1.0,b);
 87:   VecNorm(b1,NORM_2,&norm);
 88:   PetscPrintf(PETSC_COMM_WORLD,"MatSolveTranspose     : Error of norm %A\n",norm);
 89: 
 90:   /* MatSolveAdd */
 91:   MatSolveAdd(F,b,y,x);
 92:   MatMult(A,y,b1);
 93:   VecScale(b1,-1.0);
 94:   MatMultAdd(A,x,b1,b1);
 95:   VecAXPY(b1,-1.0,b);
 96:   VecNorm(b1,NORM_2,&norm);
 97:   PetscPrintf(PETSC_COMM_WORLD,"MatSolveAdd           : Error of norm %A\n",norm);
 98: 
 99:   /* MatSolveTransposeAdd */
100:   MatSolveTransposeAdd(F,b,y,x);
101:   MatMultTranspose(A,y,b1);
102:   VecScale(b1,-1.0);
103:   MatMultTransposeAdd(A,x,b1,b1);
104:   VecAXPY(b1,-1.0,b);
105:   VecNorm(b1,NORM_2,&norm);
106:   PetscPrintf(PETSC_COMM_WORLD,"MatSolveTransposeAdd  : Error of norm %A\n",norm);
107: 
108:   /* MatMatSolve */
109:   MatMatSolve(F,C,X);
110:   MatMatMult(A,X,MAT_INITIAL_MATRIX,2.0,&C1);
111:   MatAXPY(C1,-1.0,C,SAME_NONZERO_PATTERN);
112:   MatNorm(C1,NORM_FROBENIUS,&norm);
113:   PetscPrintf(PETSC_COMM_WORLD,"MatMatSolve           : Error of norm %A\n",norm);

115:   VecDestroy(x);
116:   VecDestroy(b);
117:   VecDestroy(b1);
118:   VecDestroy(y);
119:   MatDestroy(A);
120:   MatDestroy(F);
121:   MatDestroy(C);
122:   MatDestroy(C1);
123:   MatDestroy(X);
124:   ISDestroy(perm);
125:   ISDestroy(iperm);
126:   DADestroy(da);
127:   PetscFinalize();
128:   return 0;
129: }

133: PetscErrorCode ComputeRHS(DA da,Vec b)
134: {
136:   PetscInt       mx,my,mz;
137:   PetscScalar    h;

140:   DAGetInfo(da,0,&mx,&my,&mz,0,0,0,0,0,0,0);
141:   h    = 1.0/((mx-1)*(my-1)*(mz-1));
142:   VecSet(b,h);
143:   return(0);
144: }

148: PetscErrorCode ComputeRHSMatrix(PetscInt m,PetscInt nrhs,Mat* C)
149: {
151:   PetscRandom    rand;
152:   Mat            RHS;
153:   PetscScalar    *array,rval;
154:   PetscInt       i,k;

157:   MatCreate(PETSC_COMM_WORLD,&RHS);
158:   MatSetSizes(RHS,m,PETSC_DECIDE,PETSC_DECIDE,nrhs);
159:   MatSetType(RHS,MATSEQDENSE);
160: 
161:   PetscRandomCreate(PETSC_COMM_WORLD,&rand);
162:   PetscRandomSetFromOptions(rand);
163:   MatGetArray(RHS,&array);
164:   for (i=0; i<m; i++){
165:     PetscRandomGetValue(rand,&rval);
166:     array[i] = rval;
167:   }
168:   if (nrhs > 1){
169:     for (k=1; k<nrhs; k++){
170:       for (i=0; i<m; i++){
171:         array[m*k+i] = array[i];
172:       }
173:     }
174:   }
175:   MatRestoreArray(RHS,&array);
176:   MatAssemblyBegin(RHS,MAT_FINAL_ASSEMBLY);
177:   MatAssemblyEnd(RHS,MAT_FINAL_ASSEMBLY);
178:   *C = RHS;
179:   PetscRandomDestroy(rand);
180:   return(0);
181: }

183: 
186: PetscErrorCode ComputeMatrix(DA da,Mat B)
187: {
189:   PetscInt       i,j,k,mx,my,mz,xm,ym,zm,xs,ys,zs,dof,k1,k2,k3;
190:   PetscScalar    *v,*v_neighbor,Hx,Hy,Hz,HxHydHz,HyHzdHx,HxHzdHy,r1,r2;
191:   MatStencil     row,col;
192:   PetscRandom    rand;

195:   PetscRandomCreate(PETSC_COMM_WORLD,&rand);
196:   PetscRandomSetType(rand,PETSCRAND);
197:   PetscRandomSetSeed(rand,1);
198:   PetscRandomSetInterval(rand,-.001,.001);
199:   PetscRandomSetFromOptions(rand);

201:   DAGetInfo(da,0,&mx,&my,&mz,0,0,0,&dof,0,0,0);
202:   /* For simplicity, this example only works on mx=my=mz */
203:   if ( mx != my || mx != mz) SETERRQ3(1,"This example only works with mx %d = my %d = mz %d\n",mx,my,mz);

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

208:   PetscMalloc((2*dof*dof+1)*sizeof(PetscScalar),&v);
209:   v_neighbor = v + dof*dof;
210:   PetscMemzero(v,(2*dof*dof+1)*sizeof(PetscScalar));
211:   k3 = 0;
212:   for (k1=0; k1<dof; k1++){
213:     for (k2=0; k2<dof; k2++){
214:       if (k1 == k2){
215:         v[k3]          = 2.0*(HxHydHz + HxHzdHy + HyHzdHx);
216:         v_neighbor[k3] = -HxHydHz;
217:       } else {
218:         PetscRandomGetValue(rand,&r1);
219:         PetscRandomGetValue(rand,&r2);
220:         v[k3] = r1;
221:         v_neighbor[k3] = r2;
222:       }
223:       k3++;
224:     }
225:   }
226:   DAGetCorners(da,&xs,&ys,&zs,&xm,&ym,&zm);
227: 
228:   for (k=zs; k<zs+zm; k++){
229:     for (j=ys; j<ys+ym; j++){
230:       for(i=xs; i<xs+xm; i++){
231:         row.i = i; row.j = j; row.k = k;
232:         if (i==0 || j==0 || k==0 || i==mx-1 || j==my-1 || k==mz-1){ /* boudary points */
233:           MatSetValuesBlockedStencil(B,1,&row,1,&row,v,INSERT_VALUES);
234:         } else { /* interior points */
235:           /* center */
236:           col.i = i; col.j = j; col.k = k;
237:           MatSetValuesBlockedStencil(B,1,&row,1,&col,v,INSERT_VALUES);
238: 
239:           /* x neighbors */
240:           col.i = i-1; col.j = j; col.k = k;
241:           MatSetValuesBlockedStencil(B,1,&row,1,&col,v_neighbor,INSERT_VALUES);
242:           col.i = i+1; col.j = j; col.k = k;
243:           MatSetValuesBlockedStencil(B,1,&row,1,&col,v_neighbor,INSERT_VALUES);
244: 
245:           /* y neighbors */
246:           col.i = i; col.j = j-1; col.k = k;
247:           MatSetValuesBlockedStencil(B,1,&row,1,&col,v_neighbor,INSERT_VALUES);
248:           col.i = i; col.j = j+1; col.k = k;
249:           MatSetValuesBlockedStencil(B,1,&row,1,&col,v_neighbor,INSERT_VALUES);
250: 
251:           /* z neighbors */
252:           col.i = i; col.j = j; col.k = k-1;
253:           MatSetValuesBlockedStencil(B,1,&row,1,&col,v_neighbor,INSERT_VALUES);
254:           col.i = i; col.j = j; col.k = k+1;
255:           MatSetValuesBlockedStencil(B,1,&row,1,&col,v_neighbor,INSERT_VALUES);
256:         }
257:       }
258:     }
259:   }
260:   MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);
261:   MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);
262:   PetscFree(v);
263:   PetscRandomDestroy(rand);
264:   return(0);
265: }