Actual source code: test3.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 MFN interface functions.\n\n"
     12:   "The command line options are:\n"
     13:   "  -n <n>, where <n> = number of grid subdivisions in x dimension.\n\n";

     15: #include <slepcmfn.h>

     17: int main(int argc,char **argv)
     18: {
     19:   Mat                  A,B;
     20:   MFN                  mfn;
     21:   FN                   f;
     22:   MFNConvergedReason   reason;
     23:   MFNType              type;
     24:   PetscReal            norm,tol;
     25:   Vec                  v,y;
     26:   PetscInt             N,n=4,Istart,Iend,i,j,II,ncv,its,maxit;
     27:   PetscBool            flg,testprefix=PETSC_FALSE;
     28:   const char           *prefix;
     29:   PetscViewerAndFormat *vf;

     31:   PetscFunctionBeginUser;
     32:   PetscCall(SlepcInitialize(&argc,&argv,NULL,help));
     33:   PetscCall(PetscOptionsGetInt(NULL,NULL,"-n",&n,NULL));
     34:   N = n*n;
     35:   PetscCall(PetscPrintf(PETSC_COMM_WORLD,"\nSquare root of Laplacian y=sqrt(A)*e_1, N=%" PetscInt_FMT " (%" PetscInt_FMT "x%" PetscInt_FMT " grid)\n\n",N,n,n));
     36:   PetscCall(PetscOptionsGetBool(NULL,NULL,"-test_prefix",&testprefix,NULL));

     38:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
     39:                  Compute the discrete 2-D Laplacian, A
     40:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

     42:   PetscCall(MatCreate(PETSC_COMM_WORLD,&A));
     43:   PetscCall(MatSetSizes(A,PETSC_DECIDE,PETSC_DECIDE,N,N));
     44:   PetscCall(MatSetFromOptions(A));

     46:   PetscCall(MatGetOwnershipRange(A,&Istart,&Iend));
     47:   for (II=Istart;II<Iend;II++) {
     48:     i = II/n; j = II-i*n;
     49:     if (i>0) PetscCall(MatSetValue(A,II,II-n,-1.0,INSERT_VALUES));
     50:     if (i<n-1) PetscCall(MatSetValue(A,II,II+n,-1.0,INSERT_VALUES));
     51:     if (j>0) PetscCall(MatSetValue(A,II,II-1,-1.0,INSERT_VALUES));
     52:     if (j<n-1) PetscCall(MatSetValue(A,II,II+1,-1.0,INSERT_VALUES));
     53:     PetscCall(MatSetValue(A,II,II,4.0,INSERT_VALUES));
     54:   }

     56:   PetscCall(MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY));
     57:   PetscCall(MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY));
     58:   PetscCall(MatSetOption(A,MAT_HERMITIAN,PETSC_TRUE));

     60:   PetscCall(MatCreateVecs(A,&v,&y));
     61:   PetscCall(VecSetValue(v,0,1.0,INSERT_VALUES));
     62:   PetscCall(VecAssemblyBegin(v));
     63:   PetscCall(VecAssemblyEnd(v));

     65:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
     66:              Create the solver, set the matrix and the function
     67:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
     68:   PetscCall(MFNCreate(PETSC_COMM_WORLD,&mfn));
     69:   PetscCall(MFNSetOperator(mfn,A));
     70:   PetscCall(MFNGetFN(mfn,&f));
     71:   PetscCall(FNSetType(f,FNSQRT));

     73:   PetscCall(MFNSetType(mfn,MFNKRYLOV));
     74:   PetscCall(MFNGetType(mfn,&type));
     75:   PetscCall(PetscPrintf(PETSC_COMM_WORLD," Type set to %s\n",type));

     77:   /* test prefix usage */
     78:   if (testprefix) {
     79:     PetscCall(MFNSetOptionsPrefix(mfn,"check_"));
     80:     PetscCall(MFNAppendOptionsPrefix(mfn,"myprefix_"));
     81:     PetscCall(MFNGetOptionsPrefix(mfn,&prefix));
     82:     PetscCall(PetscPrintf(PETSC_COMM_WORLD," MFN prefix is currently: %s\n",prefix));
     83:   }

     85:   /* test some interface functions */
     86:   PetscCall(MFNGetOperator(mfn,&B));
     87:   PetscCall(MatView(B,PETSC_VIEWER_STDOUT_WORLD));
     88:   PetscCall(MFNSetTolerances(mfn,1e-4,500));
     89:   PetscCall(MFNSetDimensions(mfn,6));
     90:   PetscCall(MFNSetErrorIfNotConverged(mfn,PETSC_TRUE));
     91:   /* test monitors */
     92:   PetscCall(PetscViewerAndFormatCreate(PETSC_VIEWER_STDOUT_WORLD,PETSC_VIEWER_DEFAULT,&vf));
     93:   PetscCall(MFNMonitorSet(mfn,(PetscErrorCode (*)(MFN,PetscInt,PetscReal,void*))MFNMonitorDefault,vf,(PetscCtxDestroyFn*)PetscViewerAndFormatDestroy));
     94:   /* PetscCall(MFNMonitorCancel(mfn)); */
     95:   PetscCall(MFNSetFromOptions(mfn));

     97:   /* query properties and print them */
     98:   PetscCall(MFNGetTolerances(mfn,&tol,&maxit));
     99:   PetscCall(PetscPrintf(PETSC_COMM_WORLD," Tolerance: %g, max iterations: %" PetscInt_FMT "\n",(double)tol,maxit));
    100:   PetscCall(MFNGetDimensions(mfn,&ncv));
    101:   PetscCall(PetscPrintf(PETSC_COMM_WORLD," Subspace dimension: %" PetscInt_FMT "\n",ncv));
    102:   PetscCall(MFNGetErrorIfNotConverged(mfn,&flg));
    103:   if (flg) PetscCall(PetscPrintf(PETSC_COMM_WORLD," Erroring out if convergence fails\n"));

    105:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    106:                            Solve  y=sqrt(A)*v
    107:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

    109:   PetscCall(MFNSolve(mfn,v,y));
    110:   PetscCall(MFNGetConvergedReason(mfn,&reason));
    111:   PetscCall(MFNGetIterationNumber(mfn,&its));
    112:   PetscCall(PetscPrintf(PETSC_COMM_WORLD," Finished - converged reason = %d\n",(int)reason));
    113:   /* PetscCall(PetscPrintf(PETSC_COMM_WORLD," its = %" PetscInt_FMT "\n",its)); */
    114:   PetscCall(VecNorm(y,NORM_2,&norm));
    115:   PetscCall(PetscPrintf(PETSC_COMM_WORLD," sqrt(A)*v has norm %g\n",(double)norm));

    117:   /*
    118:      Free work space
    119:   */
    120:   PetscCall(MFNDestroy(&mfn));
    121:   PetscCall(MatDestroy(&A));
    122:   PetscCall(VecDestroy(&v));
    123:   PetscCall(VecDestroy(&y));
    124:   PetscCall(SlepcFinalize());
    125:   return 0;
    126: }

    128: /*TEST

    130:    testset:
    131:       args: -mfn_monitor_cancel -mfn_converged_reason -mfn_view -log_exclude mfn,bv,fn
    132:       output_file: output/test3_1.out
    133:       test:
    134:          suffix: 1
    135:       test:
    136:          suffix: 1_x
    137:          args: -mfn_monitor draw::draw_lg -draw_virtual
    138:          requires: x

    140:    test:
    141:       suffix: 2
    142:       args: -test_prefix -check_myprefix_mfn_monitor
    143:       filter: sed -e "s/estimate [0-9]\.[0-9]*e[+-]\([0-9]*\)/estimate (removed)/g" | sed -e "s/4.0[0-9]*e-10/4.03e-10/"

    145: TEST*/