Actual source code: ex67.c
petsc-3.12.4 2020-02-04
Report Typos and Errors

      2: static char help[] = "Krylov methods to solve u''  = f in parallel with periodic boundary conditions,\n\
      3:                       with a singular, inconsistent system.\n\n";

      5: /*T
      6:    Concepts: KSP^basic parallel example
      7:    Concepts: periodic boundary conditions
      8:    Processors: n
      9: T*/

     13: /*

     15:    This tests solving singular inconsistent systems with GMRES

     17:    Default: Solves a symmetric system
     18:    -symmetric false: Solves a non-symmetric system where nullspace(A) != nullspace(A')

     20:   -ksp_pc_side left or right

     22:    See the KSPSolve() for a discussion of when right preconditioning with nullspace(A) != nullspace(A') can fail to produce the
     23:    norm minimizing solution.

     25:    Note that though this example does solve the system with right preconditioning and nullspace(A) != nullspace(A') it does not produce the
     26:    norm minimizing solution, that is the computed solution is not orthogonal to the nullspace(A).

     29:    Include "petscdmda.h" so that we can use distributed arrays (DMDAs).
     30:    Include "petscksp.h" so that we can use KSP solvers.  Note that this
     31:    file automatically includes:
     32:      petscsys.h       - base PETSc routines   petscvec.h - vectors
     33:      petscmat.h - matrices
     34:      petscis.h     - index sets            petscksp.h - Krylov subspace methods
     35:      petscviewer.h - viewers               petscpc.h  - preconditioners
     36:      petscksp.h   - linear solvers
     37: */

     39:  #include <petscdm.h>
     40:  #include <petscdmda.h>
     41:  #include <petscsnes.h>
     42:  #include <petsc/private/kspimpl.h>

     44: PetscBool symmetric = PETSC_TRUE;

     46: PetscErrorCode FormMatrix(Mat,void*);
     47: PetscErrorCode FormRightHandSide(Vec,void*);

     49: int main(int argc,char **argv)
     50: {
     51:   KSP            ksp;
     52:   Mat            J;
     53:   DM             da;
     54:   Vec            x,r;              /* vectors */
     56:   PetscInt       M = 10;
     57:   MatNullSpace   constants,nconstants;

     59:   PetscInitialize(&argc,&argv,(char*)0,help);if (ierr) return ierr;
     60:   PetscOptionsGetInt(NULL,NULL,"-M",&M,NULL);
     61:   PetscOptionsGetBool(NULL,NULL,"-symmetric",&symmetric,NULL);

     63:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
     64:      Create linear solver context
     65:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

     67:   KSPCreate(PETSC_COMM_WORLD,&ksp);

     69:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
     70:      Create vector data structures; set function evaluation routine
     71:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

     73:   /*
     74:      Create distributed array (DMDA) to manage parallel grid and vectors
     75:   */
     76:   DMDACreate1d(PETSC_COMM_WORLD,DM_BOUNDARY_PERIODIC,M,1,2,NULL,&da);
     77:   DMSetFromOptions(da);
     78:   DMSetUp(da);

     80:   /*
     81:      Extract global and local vectors from DMDA; then duplicate for remaining
     82:      vectors that are the same types
     83:   */
     84:   DMCreateGlobalVector(da,&x);
     85:   VecDuplicate(x,&r);

     87:   /*
     88:      Set function evaluation routine and vector.  Whenever the nonlinear
     89:      solver needs to compute the nonlinear function, it will call this
     90:      routine.
     91:       - Note that the final routine argument is the user-defined
     92:         context that provides application-specific data for the
     93:         function evaluation routine.
     94:   */
     95:   FormRightHandSide(r,da);

     97:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
     98:      Create matrix data structure;
     99:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
    100:   DMCreateMatrix(da,&J);
    101:   MatNullSpaceCreate(PETSC_COMM_WORLD,PETSC_TRUE,0,NULL,&constants);
    102:   if (symmetric) {
    103:     MatSetOption(J,MAT_SYMMETRY_ETERNAL,PETSC_TRUE);
    104:     MatSetOption(J,MAT_SYMMETRIC,PETSC_TRUE);
    105:   } else {
    106:     Vec         n;
    107:     PetscInt    zero = 0;
    108:     PetscScalar zeros = 0.0;
    109:     VecDuplicate(x,&n);
    110:     /* the nullspace(A') is the constant vector but with a zero in the very first entry; hence nullspace(A') != nullspace(A) */
    111:     VecSet(n,1.0);
    112:     VecSetValues(n,1,&zero,&zeros,INSERT_VALUES);
    113:     VecAssemblyBegin(n);
    114:     VecAssemblyEnd(n);
    115:     VecNormalize(n,NULL);
    116:     MatNullSpaceCreate(PETSC_COMM_WORLD,PETSC_FALSE,1,&n,&nconstants);
    117:     MatSetTransposeNullSpace(J,nconstants);
    118:     MatNullSpaceDestroy(&nconstants);
    119:     VecDestroy(&n);
    120:   }
    121:   MatSetNullSpace(J,constants);
    122:   FormMatrix(J,da);

    124:   KSPSetOperators(ksp,J,J);

    126:   KSPSetFromOptions(ksp);
    127:   KSPSolve(ksp,r,x);
    128:   KSPSolveTranspose(ksp,r,x);

    130:   VecDestroy(&x);
    131:   VecDestroy(&r);
    132:   MatDestroy(&J);
    133:   MatNullSpaceDestroy(&constants);
    134:   KSPDestroy(&ksp);
    135:   DMDestroy(&da);
    136:   PetscFinalize();
    137:   return ierr;
    138: }

    140: /*

    142:      This intentionally includes something in the right hand side that is not in the range of the matrix A.
    143:      Since MatSetNullSpace() is called and the matrix is symmetric; the Krylov method automatically removes this
    144:      portion of the right hand side before solving the linear system.
    145: */
    146: PetscErrorCode FormRightHandSide(Vec f,void *ctx)
    147: {
    148:   DM             da    = (DM) ctx;
    149:   PetscScalar    *ff;
    151:   PetscInt       i,M,xs,xm;
    152:   PetscReal      h;

    155:   DMDAVecGetArray(da,f,&ff);

    157:   /*
    158:      Get local grid boundaries (for 1-dimensional DMDA):
    159:        xs, xm  - starting grid index, width of local grid (no ghost points)
    160:   */
    161:   DMDAGetCorners(da,&xs,NULL,NULL,&xm,NULL,NULL);
    162:   DMDAGetInfo(da,NULL,&M,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL);

    164:   /*
    165:      Compute function over locally owned part of the grid
    166:      Note the [i-1] and [i+1] will automatically access the ghost points from other processes or the periodic points.
    167:   */
    168:   h = 1.0/M;
    169:   for (i=xs; i<xs+xm; i++) ff[i] = - PetscSinReal(2.0*PETSC_PI*i*h) + 1.0;

    171:   /*
    172:      Restore vectors
    173:   */
    174:   DMDAVecRestoreArray(da,f,&ff);
    175:   return(0);
    176: }
    177: /* ------------------------------------------------------------------- */
    178: PetscErrorCode FormMatrix(Mat jac,void *ctx)
    179: {
    180:   PetscScalar    A[3];
    182:   PetscInt       i,M,xs,xm;
    183:   DM             da = (DM) ctx;
    184:   MatStencil     row,cols[3];
    185:   PetscReal      h;

    188:   DMDAGetCorners(da,&xs,NULL,NULL,&xm,NULL,NULL);

    190:   /*
    191:     Get range of locally owned matrix
    192:   */
    193:   DMDAGetInfo(da,NULL,&M,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL);

    195:   MatZeroEntries(jac);
    196:   h = 1.0/M;
    197:   /* because of periodic boundary conditions we can simply loop over all local nodes and access to the left and right */
    198:   if (symmetric) {
    199:     for (i=xs; i<xs+xm; i++) {
    200:       row.i = i;
    201:       cols[0].i = i - 1;
    202:       cols[1].i = i;
    203:       cols[2].i = i + 1;
    204:       A[0] = A[2] = 1.0/(h*h); A[1] = -2.0/(h*h);
    205:       MatSetValuesStencil(jac,1,&row,3,cols,A,ADD_VALUES);
    206:     }
    207:   } else {
    208:     MatStencil  *acols;
    209:     PetscScalar *avals;

    211:     /* only works for one process */
    212:     MatSetOption(jac,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_FALSE);
    213:     row.i = 0;
    214:     PetscMalloc1(M,&acols);
    215:     PetscMalloc1(M,&avals);
    216:     for (i=0; i<M; i++) {
    217:       acols[i].i = i;
    218:       avals[i]  = (i % 2) ? 1 : -1;
    219:     }
    220:     MatSetValuesStencil(jac,1,&row,M,acols,avals,ADD_VALUES);
    221:     PetscFree(acols);
    222:     PetscFree(avals);
    223:     row.i = 1;
    224:     cols[0].i = - 1;
    225:     cols[1].i = 1;
    226:     cols[2].i = 1 + 1;
    227:     A[0] = A[2] = 1.0/(h*h); A[1] = -2.0/(h*h);
    228:     MatSetValuesStencil(jac,1,&row,3,cols,A,ADD_VALUES);
    229:     for (i=2; i<xs+xm-1; i++) {
    230:       row.i = i;
    231:       cols[0].i = i - 1;
    232:       cols[1].i = i;
    233:       cols[2].i = i + 1;
    234:       A[0] = A[2] = 1.0/(h*h); A[1] = -2.0/(h*h);
    235:       MatSetValuesStencil(jac,1,&row,3,cols,A,ADD_VALUES);
    236:     }
    237:     row.i = M - 1 ;
    238:     cols[0].i = M-2;
    239:     cols[1].i = M-1;
    240:     cols[2].i = M+1;
    241:     A[0] = A[2] = 1.0/(h*h); A[1] = -2.0/(h*h);
    242:     MatSetValuesStencil(jac,1,&row,3,cols,A,ADD_VALUES);
    243:   }
    244:   MatAssemblyBegin(jac,MAT_FINAL_ASSEMBLY);
    245:   MatAssemblyEnd(jac,MAT_FINAL_ASSEMBLY);
    246:   return(0);
    247: }

    250: /*TEST

    252:    test:
    253:       suffix: nonsymmetric_left
    254:       args: -symmetric false -ksp_view -ksp_converged_reason -pc_type jacobi -mat_no_inode -ksp_monitor_true_residual -ksp_rtol 1.e-14 -ksp_max_it 12 -ksp_pc_side left
    255:       filter: sed 's/ATOL/RTOL/g'
    256:       requires: !single

    258:    test:
    259:       suffix: nonsymmetric_right
    260:       args: -symmetric false -ksp_view -ksp_converged_reason -pc_type jacobi -mat_no_inode -ksp_monitor_true_residual -ksp_rtol 1.e-14 -ksp_max_it 12 -ksp_pc_side right
    261:       filter: sed 's/ATOL/RTOL/g'
    262:       requires: !single

    264:    test:
    265:       suffix: symmetric_left
    266:       args: -ksp_view -ksp_converged_reason -pc_type sor -mat_no_inode -ksp_monitor_true_residual -ksp_rtol 1.e-14 -ksp_max_it 12 -ksp_pc_side left
    267:       requires: !single

    269:    test:
    270:       suffix: symmetric_right
    271:       args: -ksp_view -ksp_converged_reason -pc_type sor -mat_no_inode -ksp_monitor_true_residual -ksp_rtol 1.e-14 -ksp_max_it 12 -ksp_pc_side right
    272:       filter: sed 's/ATOL/RTOL/g'
    273:       requires: !single

    275:    test:
    276:       suffix: transpose_asm
    277:       args: -symmetric false -ksp_monitor -ksp_view -pc_type asm -sub_pc_type lu -sub_pc_factor_zeropivot 1.e-33 -ksp_converged_reason
    278:       filter: sed 's/ATOL/RTOL/g'

    280: TEST*/