Line data    Source code

       1             : /*
       2             :    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
       3             :    SLEPc - Scalable Library for Eigenvalue Problem Computations
       4             :    Copyright (c) 2002-, Universitat Politecnica de Valencia, Spain
       5             : 
       6             :    This file is part of SLEPc.
       7             :    SLEPc is distributed under a 2-clause BSD license (see LICENSE).
       8             :    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
       9             : */
      10             : 
      11             : static char help[] = "Tests a user-defined convergence test in NEP.\n\n"
      12             :   "The command line options are:\n"
      13             :   "  -n <n>, where <n> = matrix dimension.\n";
      14             : 
      15             : /*
      16             :    Solve T(lambda)x=0 with T(lambda) = -D+sqrt(lambda)*I
      17             :       where D is the Laplacian operator in 1 dimension
      18             : */
      19             : 
      20             : #include <slepcnep.h>
      21             : 
      22             : /*
      23             :   MyConvergedRel - Convergence test relative to the norm of D (given in ctx).
      24             : */
      25           7 : PetscErrorCode MyConvergedRel(NEP nep,PetscScalar eigr,PetscScalar eigi,PetscReal res,PetscReal *errest,void *ctx)
      26             : {
      27           7 :   PetscReal norm = *(PetscReal*)ctx;
      28             : 
      29           7 :   PetscFunctionBegin;
      30           7 :   *errest = res/norm;
      31           7 :   PetscFunctionReturn(PETSC_SUCCESS);
      32             : }
      33             : 
      34           1 : int main(int argc,char **argv)
      35             : {
      36           1 :   NEP            nep;             /* nonlinear eigensolver context */
      37           1 :   Mat            A[2];
      38           1 :   PetscInt       n=100,Istart,Iend,i;
      39           1 :   PetscBool      terse;
      40           1 :   PetscReal      norm;
      41           1 :   FN             f[2];
      42           1 :   PetscScalar    coeffs;
      43             : 
      44           1 :   PetscFunctionBeginUser;
      45           1 :   PetscCall(SlepcInitialize(&argc,&argv,NULL,help));
      46           1 :   PetscCall(PetscOptionsGetInt(NULL,NULL,"-n",&n,NULL));
      47           1 :   PetscCall(PetscPrintf(PETSC_COMM_WORLD,"\nSquare root eigenproblem, n=%" PetscInt_FMT "\n\n",n));
      48             : 
      49             :   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
      50             :      Create nonlinear eigensolver, define problem in split form
      51             :      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
      52             : 
      53           1 :   PetscCall(NEPCreate(PETSC_COMM_WORLD,&nep));
      54             : 
      55             :   /* Create matrices */
      56           1 :   PetscCall(MatCreate(PETSC_COMM_WORLD,&A[0]));
      57           1 :   PetscCall(MatSetSizes(A[0],PETSC_DECIDE,PETSC_DECIDE,n,n));
      58           1 :   PetscCall(MatSetFromOptions(A[0]));
      59           1 :   PetscCall(MatGetOwnershipRange(A[0],&Istart,&Iend));
      60         101 :   for (i=Istart;i<Iend;i++) {
      61         100 :     if (i>0) PetscCall(MatSetValue(A[0],i,i-1,1.0,INSERT_VALUES));
      62         100 :     if (i<n-1) PetscCall(MatSetValue(A[0],i,i+1,1.0,INSERT_VALUES));
      63         100 :     PetscCall(MatSetValue(A[0],i,i,-2.0,INSERT_VALUES));
      64             :   }
      65           1 :   PetscCall(MatAssemblyBegin(A[0],MAT_FINAL_ASSEMBLY));
      66           1 :   PetscCall(MatAssemblyEnd(A[0],MAT_FINAL_ASSEMBLY));
      67             : 
      68           1 :   PetscCall(MatCreateConstantDiagonal(PETSC_COMM_WORLD,PETSC_DECIDE,PETSC_DECIDE,n,n,1.0,&A[1]));
      69             : 
      70             :   /* Define functions */
      71           1 :   PetscCall(FNCreate(PETSC_COMM_WORLD,&f[0]));
      72           1 :   PetscCall(FNSetType(f[0],FNRATIONAL));
      73           1 :   coeffs = 1.0;
      74           1 :   PetscCall(FNRationalSetNumerator(f[0],1,&coeffs));
      75           1 :   PetscCall(FNCreate(PETSC_COMM_WORLD,&f[1]));
      76           1 :   PetscCall(FNSetType(f[1],FNSQRT));
      77           1 :   PetscCall(NEPSetSplitOperator(nep,2,A,f,SUBSET_NONZERO_PATTERN));
      78             : 
      79             :   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
      80             :                    Set some options and solve
      81             :      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
      82             : 
      83           1 :   PetscCall(NEPSetTarget(nep,1.1));
      84             : 
      85             :   /* setup convergence test relative to the norm of D */
      86           1 :   PetscCall(MatNorm(A[0],NORM_1,&norm));
      87           1 :   PetscCall(NEPSetConvergenceTestFunction(nep,MyConvergedRel,&norm,NULL));
      88             : 
      89           1 :   PetscCall(NEPSetFromOptions(nep));
      90           1 :   PetscCall(NEPSolve(nep));
      91             : 
      92             :   /* show detailed info unless -terse option is given by user */
      93           1 :   PetscCall(PetscOptionsHasName(NULL,NULL,"-terse",&terse));
      94           1 :   if (terse) PetscCall(NEPErrorView(nep,NEP_ERROR_BACKWARD,NULL));
      95             :   else {
      96           0 :     PetscCall(PetscViewerPushFormat(PETSC_VIEWER_STDOUT_WORLD,PETSC_VIEWER_ASCII_INFO_DETAIL));
      97           0 :     PetscCall(NEPConvergedReasonView(nep,PETSC_VIEWER_STDOUT_WORLD));
      98           0 :     PetscCall(NEPErrorView(nep,NEP_ERROR_BACKWARD,PETSC_VIEWER_STDOUT_WORLD));
      99           0 :     PetscCall(PetscViewerPopFormat(PETSC_VIEWER_STDOUT_WORLD));
     100             :   }
     101           1 :   PetscCall(NEPDestroy(&nep));
     102           1 :   PetscCall(MatDestroy(&A[0]));
     103           1 :   PetscCall(MatDestroy(&A[1]));
     104           1 :   PetscCall(FNDestroy(&f[0]));
     105           1 :   PetscCall(FNDestroy(&f[1]));
     106           1 :   PetscCall(SlepcFinalize());
     107             :   return 0;
     108             : }
     109             : 
     110             : /*TEST
     111             : 
     112             :    test:
     113             :       suffix: 1
     114             :       args: -nep_type slp -nep_nev 2 -terse
     115             : 
     116             : TEST*/