Actual source code: ex129.c
petsc-3.4.2 2013-07-02

      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>

     17: extern PetscErrorCode ComputeMatrix(DM,Mat);
     18: extern PetscErrorCode ComputeRHS(DM,Vec);
     19: extern PetscErrorCode ComputeRHSMatrix(PetscInt,PetscInt,Mat*);

     23: int main(int argc,char **args)
     24: {
     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,tol=1.e-13;
     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(NULL,"-dof",&dof,NULL);
     41:   PetscOptionsGetInt(NULL,"-M",&M,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,NULL,NULL,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:   DMCreateMatrix(da,MATBAIJ,&A);
     61:   ComputeMatrix(da,A);
     62:   MatGetSize(A,&m,&n);
     63:   nrhs = 2;
     64:   PetscOptionsGetInt(NULL,"-nrhs",&nrhs,NULL);
     65:   ComputeRHSMatrix(m,nrhs,&RHS);
     66:   MatDuplicate(RHS,MAT_DO_NOT_COPY_VALUES,&X);

     68:   MatGetOrdering(A,MATORDERINGND,&perm,&iperm);

     71:   PetscOptionsGetBool(NULL,"-inplacelu",&InplaceLU,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:     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(PetscObjectComm((PetscObject)A),&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);

    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:   if (norm > tol) {
    109:     PetscPrintf(PETSC_COMM_WORLD,"MatSolve              : Error of norm %G\n",norm);
    110:   }

    112:   /* MatSolveTranspose */
    113:   MatSolveTranspose(F,b,x);
    114:   MatMultTranspose(A,x,b1);
    115:   VecAXPY(b1,-1.0,b);
    116:   VecNorm(b1,NORM_2,&norm);
    117:   if (norm > tol) {
    118:     PetscPrintf(PETSC_COMM_WORLD,"MatSolveTranspose     : Error of norm %G\n",norm);
    119:   }

    121:   /* MatSolveAdd */
    122:   MatSolveAdd(F,b,y,x);
    123:   MatMult(A,y,b1);
    124:   VecScale(b1,-1.0);
    125:   MatMultAdd(A,x,b1,b1);
    126:   VecAXPY(b1,-1.0,b);
    127:   VecNorm(b1,NORM_2,&norm);
    128:   if (norm > tol) {
    129:     PetscPrintf(PETSC_COMM_WORLD,"MatSolveAdd           : Error of norm %G\n",norm);
    130:   }

    132:   /* MatSolveTransposeAdd */
    133:   MatSolveTransposeAdd(F,b,y,x);
    134:   MatMultTranspose(A,y,b1);
    135:   VecScale(b1,-1.0);
    136:   MatMultTransposeAdd(A,x,b1,b1);
    137:   VecAXPY(b1,-1.0,b);
    138:   VecNorm(b1,NORM_2,&norm);
    139:   if (norm > tol) {
    140:     PetscPrintf(PETSC_COMM_WORLD,"MatSolveTransposeAdd  : Error of norm %G\n",norm);
    141:   }

    143:   /* MatMatSolve */
    144:   MatMatSolve(F,RHS,X);
    145:   MatMatMult(A,X,MAT_INITIAL_MATRIX,2.0,&C1);
    146:   MatAXPY(C1,-1.0,RHS,SAME_NONZERO_PATTERN);
    147:   MatNorm(C1,NORM_FROBENIUS,&norm);
    148:   if (norm > tol) {
    149:     PetscPrintf(PETSC_COMM_WORLD,"MatMatSolve           : Error of norm %G\n",norm);
    150:   }

    152:   VecDestroy(&x);
    153:   VecDestroy(&b);
    154:   VecDestroy(&b1);
    155:   VecDestroy(&y);
    156:   MatDestroy(&A);
    157:   MatDestroy(&F);
    158:   MatDestroy(&RHS);
    159:   MatDestroy(&C1);
    160:   MatDestroy(&X);
    161:   ISDestroy(&perm);
    162:   ISDestroy(&iperm);
    163:   DMDestroy(&da);
    164:   PetscFinalize();
    165:   return 0;
    166: }

    170: PetscErrorCode ComputeRHS(DM da,Vec b)
    171: {
    173:   PetscInt       mx,my,mz;
    174:   PetscScalar    h;

    177:   DMDAGetInfo(da,0,&mx,&my,&mz,0,0,0,0,0,0,0,0,0);
    178:   h    = 1.0/((mx-1)*(my-1)*(mz-1));
    179:   VecSet(b,h);
    180:   return(0);
    181: }

    185: PetscErrorCode ComputeRHSMatrix(PetscInt m,PetscInt nrhs,Mat *C)
    186: {
    188:   PetscRandom    rand;
    189:   Mat            RHS;
    190:   PetscScalar    *array,rval;
    191:   PetscInt       i,k;

    194:   MatCreate(PETSC_COMM_WORLD,&RHS);
    195:   MatSetSizes(RHS,m,PETSC_DECIDE,PETSC_DECIDE,nrhs);
    196:   MatSetType(RHS,MATSEQDENSE);
    197:   MatSetUp(RHS);

    199:   PetscRandomCreate(PETSC_COMM_WORLD,&rand);
    200:   PetscRandomSetFromOptions(rand);
    201:   MatDenseGetArray(RHS,&array);
    202:   for (i=0; i<m; i++) {
    203:     PetscRandomGetValue(rand,&rval);
    204:     array[i] = rval;
    205:   }
    206:   if (nrhs > 1) {
    207:     for (k=1; k<nrhs; k++) {
    208:       for (i=0; i<m; i++) {
    209:         array[m*k+i] = array[i];
    210:       }
    211:     }
    212:   }
    213:   MatDenseRestoreArray(RHS,&array);
    214:   MatAssemblyBegin(RHS,MAT_FINAL_ASSEMBLY);
    215:   MatAssemblyEnd(RHS,MAT_FINAL_ASSEMBLY);
    216:   *C   = RHS;
    217:   PetscRandomDestroy(&rand);
    218:   return(0);
    219: }

    224: PetscErrorCode ComputeMatrix(DM da,Mat B)
    225: {
    227:   PetscInt       i,j,k,mx,my,mz,xm,ym,zm,xs,ys,zs,dof,k1,k2,k3;
    228:   PetscScalar    *v,*v_neighbor,Hx,Hy,Hz,HxHydHz,HyHzdHx,HxHzdHy,r1,r2;
    229:   MatStencil     row,col;
    230:   PetscRandom    rand;

    233:   PetscRandomCreate(PETSC_COMM_WORLD,&rand);
    234:   PetscRandomSetType(rand,PETSCRAND);
    235:   PetscRandomSetSeed(rand,1);
    236:   PetscRandomSetInterval(rand,-.001,.001);
    237:   PetscRandomSetFromOptions(rand);

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

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

    246:   PetscMalloc((2*dof*dof+1)*sizeof(PetscScalar),&v);
    247:   v_neighbor = v + dof*dof;
    248:   PetscMemzero(v,(2*dof*dof+1)*sizeof(PetscScalar));
    249:   k3         = 0;
    250:   for (k1=0; k1<dof; k1++) {
    251:     for (k2=0; k2<dof; k2++) {
    252:       if (k1 == k2) {
    253:         v[k3]          = 2.0*(HxHydHz + HxHzdHy + HyHzdHx);
    254:         v_neighbor[k3] = -HxHydHz;
    255:       } else {
    256:         PetscRandomGetValue(rand,&r1);
    257:         PetscRandomGetValue(rand,&r2);

    259:         v[k3]          = r1;
    260:         v_neighbor[k3] = r2;
    261:       }
    262:       k3++;
    263:     }
    264:   }
    265:   DMDAGetCorners(da,&xs,&ys,&zs,&xm,&ym,&zm);

    267:   for (k=zs; k<zs+zm; k++) {
    268:     for (j=ys; j<ys+ym; j++) {
    269:       for (i=xs; i<xs+xm; i++) {
    270:         row.i = i; row.j = j; row.k = k;
    271:         if (i==0 || j==0 || k==0 || i==mx-1 || j==my-1 || k==mz-1) { /* boudary points */
    272:           MatSetValuesBlockedStencil(B,1,&row,1,&row,v,INSERT_VALUES);
    273:         } else { /* interior points */
    274:           /* center */
    275:           col.i = i; col.j = j; col.k = k;
    276:           MatSetValuesBlockedStencil(B,1,&row,1,&col,v,INSERT_VALUES);

    278:           /* x neighbors */
    279:           col.i = i-1; col.j = j; col.k = k;
    280:           MatSetValuesBlockedStencil(B,1,&row,1,&col,v_neighbor,INSERT_VALUES);
    281:           col.i = i+1; col.j = j; col.k = k;
    282:           MatSetValuesBlockedStencil(B,1,&row,1,&col,v_neighbor,INSERT_VALUES);

    284:           /* y neighbors */
    285:           col.i = i; col.j = j-1; col.k = k;
    286:           MatSetValuesBlockedStencil(B,1,&row,1,&col,v_neighbor,INSERT_VALUES);
    287:           col.i = i; col.j = j+1; col.k = k;
    288:           MatSetValuesBlockedStencil(B,1,&row,1,&col,v_neighbor,INSERT_VALUES);

    290:           /* z neighbors */
    291:           col.i = i; col.j = j; col.k = k-1;
    292:           MatSetValuesBlockedStencil(B,1,&row,1,&col,v_neighbor,INSERT_VALUES);
    293:           col.i = i; col.j = j; col.k = k+1;
    294:           MatSetValuesBlockedStencil(B,1,&row,1,&col,v_neighbor,INSERT_VALUES);
    295:         }
    296:       }
    297:     }
    298:   }
    299:   MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);
    300:   MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);
    301:   PetscFree(v);
    302:   PetscRandomDestroy(&rand);
    303:   return(0);
    304: }