Actual source code: test6.c

      1: /*
      2:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
      3:    SLEPc - Scalable Library for Eigenvalue Problem Computations
      4:    Copyright (c) 2002-, Universitat Politecnica de Valencia, Spain

      6:    This file is part of SLEPc.
      7:    SLEPc is distributed under a 2-clause BSD license (see LICENSE).
      8:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
      9: */

     11: static char help[] = "Test STPRECOND operations.\n\n";

     13: #include <slepcst.h>

     15: int main(int argc,char **argv)
     16: {
     17:   Mat            A,P,mat[1];
     18:   ST             st;
     19:   KSP            ksp;
     20:   Vec            v,w;
     21:   PetscScalar    sigma;
     22:   PetscInt       n=10,i,Istart,Iend;
     23:   STMatMode      matmode;

     25:   PetscFunctionBeginUser;
     26:   PetscCall(SlepcInitialize(&argc,&argv,NULL,help));
     27:   PetscCall(PetscOptionsGetInt(NULL,NULL,"-n",&n,NULL));
     28:   PetscCall(PetscPrintf(PETSC_COMM_WORLD,"\nPreconditioner for 1-D Laplacian, n=%" PetscInt_FMT "\n\n",n));

     30:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
     31:              Compute the operator matrix for the 1-D Laplacian
     32:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

     34:   PetscCall(MatCreate(PETSC_COMM_WORLD,&A));
     35:   PetscCall(MatSetSizes(A,PETSC_DECIDE,PETSC_DECIDE,n,n));
     36:   PetscCall(MatSetFromOptions(A));

     38:   PetscCall(MatGetOwnershipRange(A,&Istart,&Iend));
     39:   for (i=Istart;i<Iend;i++) {
     40:     if (i>0) PetscCall(MatSetValue(A,i,i-1,-1.0,INSERT_VALUES));
     41:     if (i<n-1) PetscCall(MatSetValue(A,i,i+1,-1.0,INSERT_VALUES));
     42:     PetscCall(MatSetValue(A,i,i,2.0,INSERT_VALUES));
     43:   }
     44:   PetscCall(MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY));
     45:   PetscCall(MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY));
     46:   PetscCall(MatCreateVecs(A,&v,&w));
     47:   PetscCall(VecSet(v,1.0));

     49:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
     50:                 Create the spectral transformation object
     51:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

     53:   PetscCall(STCreate(PETSC_COMM_WORLD,&st));
     54:   mat[0] = A;
     55:   PetscCall(STSetMatrices(st,1,mat));
     56:   PetscCall(STSetType(st,STPRECOND));
     57:   PetscCall(STGetKSP(st,&ksp));
     58:   PetscCall(KSPSetType(ksp,KSPPREONLY));
     59:   PetscCall(STSetFromOptions(st));

     61:   /* set up */
     62:   /* - the transform flag is necessary so that A-sigma*I is built explicitly */
     63:   PetscCall(STSetTransform(st,PETSC_TRUE));
     64:   /* - the ksphasmat flag is necessary when using STApply(), otherwise can only use PCApply() */
     65:   PetscCall(STPrecondSetKSPHasMat(st,PETSC_TRUE));
     66:   /* no need to call STSetUp() explicitly */
     67:   PetscCall(STSetUp(st));

     69:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
     70:                        Apply the preconditioner
     71:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

     73:   /* default shift */
     74:   PetscCall(PetscPrintf(PETSC_COMM_WORLD,"With default shift\n"));
     75:   PetscCall(STApply(st,v,w));
     76:   PetscCall(VecView(w,NULL));

     78:   /* change shift */
     79:   sigma = 0.1;
     80:   PetscCall(STSetShift(st,sigma));
     81:   PetscCall(STGetShift(st,&sigma));
     82:   PetscCall(PetscPrintf(PETSC_COMM_WORLD,"With shift=%g\n",(double)PetscRealPart(sigma)));
     83:   PetscCall(STApply(st,v,w));
     84:   PetscCall(VecView(w,NULL));
     85:   PetscCall(STPostSolve(st));   /* undo changes if inplace */

     87:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
     88:                  Test a user-provided preconditioner matrix
     89:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

     91:   PetscCall(MatCreate(PETSC_COMM_WORLD,&P));
     92:   PetscCall(MatSetSizes(P,PETSC_DECIDE,PETSC_DECIDE,n,n));
     93:   PetscCall(MatSetFromOptions(P));

     95:   PetscCall(MatGetOwnershipRange(P,&Istart,&Iend));
     96:   for (i=Istart;i<Iend;i++) PetscCall(MatSetValue(P,i,i,2.0,INSERT_VALUES));
     97:   if (Istart==0) {
     98:     PetscCall(MatSetValue(P,1,0,-1.0,INSERT_VALUES));
     99:     PetscCall(MatSetValue(P,0,1,-1.0,INSERT_VALUES));
    100:   }
    101:   PetscCall(MatAssemblyBegin(P,MAT_FINAL_ASSEMBLY));
    102:   PetscCall(MatAssemblyEnd(P,MAT_FINAL_ASSEMBLY));

    104:   /* apply new preconditioner */
    105:   PetscCall(STSetPreconditionerMat(st,P));
    106:   PetscCall(PetscPrintf(PETSC_COMM_WORLD,"With user-provided matrix\n"));
    107:   PetscCall(STApply(st,v,w));
    108:   PetscCall(VecView(w,NULL));

    110:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    111:             Test a user-provided preconditioner in split form
    112:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

    114:   PetscCall(STGetMatMode(st,&matmode));
    115:   if (matmode==ST_MATMODE_COPY) {
    116:     PetscCall(STSetPreconditionerMat(st,NULL));
    117:     mat[0] = P;
    118:     PetscCall(STSetSplitPreconditioner(st,1,mat,SAME_NONZERO_PATTERN));

    120:     /* apply new preconditioner */
    121:     PetscCall(PetscPrintf(PETSC_COMM_WORLD,"With split preconditioner\n"));
    122:     PetscCall(STApply(st,v,w));
    123:     PetscCall(VecView(w,NULL));
    124:   }

    126:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    127:                              Clean up
    128:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
    129:   PetscCall(STDestroy(&st));
    130:   PetscCall(MatDestroy(&A));
    131:   PetscCall(MatDestroy(&P));
    132:   PetscCall(VecDestroy(&v));
    133:   PetscCall(VecDestroy(&w));
    134:   PetscCall(SlepcFinalize());
    135:   return 0;
    136: }

    138: /*TEST

    140:    test:
    141:       suffix: 1
    142:       args: -st_matmode {{copy inplace shell}separate output}
    143:       requires: !single

    145: TEST*/