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[] = "Test EPSSetArbitrarySelection.\n\n";
      12             : 
      13             : #include <slepceps.h>
      14             : 
      15       16260 : PetscErrorCode MyArbitrarySelection(PetscScalar eigr,PetscScalar eigi,Vec xr,Vec xi,PetscScalar *rr,PetscScalar *ri,void *ctx)
      16             : {
      17       16260 :   Vec             xref = *(Vec*)ctx;
      18             : 
      19       16260 :   PetscFunctionBeginUser;
      20       16260 :   PetscCall(VecDot(xr,xref,rr));
      21       16260 :   *rr = PetscAbsScalar(*rr);
      22       16260 :   if (ri) *ri = 0.0;
      23       16260 :   PetscFunctionReturn(PETSC_SUCCESS);
      24             : }
      25             : 
      26           8 : int main(int argc,char **argv)
      27             : {
      28           8 :   Mat            A;           /* problem matrices */
      29           8 :   EPS            eps;         /* eigenproblem solver context */
      30           8 :   PetscScalar    seigr,seigi;
      31           8 :   PetscReal      tol=1000*PETSC_MACHINE_EPSILON;
      32           8 :   Vec            sxr,sxi;
      33           8 :   PetscInt       n=30,i,Istart,Iend,nconv;
      34             : 
      35           8 :   PetscFunctionBeginUser;
      36           8 :   PetscCall(SlepcInitialize(&argc,&argv,NULL,help));
      37             : 
      38           8 :   PetscCall(PetscOptionsGetInt(NULL,NULL,"-n",&n,NULL));
      39           8 :   PetscCall(PetscPrintf(PETSC_COMM_WORLD,"\nTridiagonal with zero diagonal, n=%" PetscInt_FMT "\n\n",n));
      40             : 
      41             :   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
      42             :            Create matrix tridiag([-1 0 -1])
      43             :      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
      44           8 :   PetscCall(MatCreate(PETSC_COMM_WORLD,&A));
      45           8 :   PetscCall(MatSetSizes(A,PETSC_DECIDE,PETSC_DECIDE,n,n));
      46           8 :   PetscCall(MatSetFromOptions(A));
      47             : 
      48           8 :   PetscCall(MatGetOwnershipRange(A,&Istart,&Iend));
      49         248 :   for (i=Istart;i<Iend;i++) {
      50         240 :     if (i>0) PetscCall(MatSetValue(A,i,i-1,-1.0,INSERT_VALUES));
      51         240 :     if (i<n-1) PetscCall(MatSetValue(A,i,i+1,-1.0,INSERT_VALUES));
      52             :   }
      53           8 :   PetscCall(MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY));
      54           8 :   PetscCall(MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY));
      55             : 
      56             :   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
      57             :                         Create the eigensolver
      58             :      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
      59           8 :   PetscCall(EPSCreate(PETSC_COMM_WORLD,&eps));
      60           8 :   PetscCall(EPSSetProblemType(eps,EPS_HEP));
      61           8 :   PetscCall(EPSSetTolerances(eps,tol,PETSC_CURRENT));
      62           8 :   PetscCall(EPSSetOperators(eps,A,NULL));
      63           8 :   PetscCall(EPSSetWhichEigenpairs(eps,EPS_SMALLEST_REAL));
      64           8 :   PetscCall(EPSSetFromOptions(eps));
      65             : 
      66             :   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
      67             :                 Solve eigenproblem and store some solution
      68             :      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
      69           8 :   PetscCall(EPSSolve(eps));
      70           8 :   PetscCall(MatCreateVecs(A,&sxr,NULL));
      71           8 :   PetscCall(MatCreateVecs(A,&sxi,NULL));
      72           8 :   PetscCall(EPSGetConverged(eps,&nconv));
      73           8 :   if (nconv>0) {
      74           8 :     PetscCall(EPSGetEigenpair(eps,0,&seigr,&seigi,sxr,sxi));
      75           8 :     PetscCall(EPSErrorView(eps,EPS_ERROR_RELATIVE,NULL));
      76             : 
      77             :     /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
      78             :                  Solve eigenproblem using an arbitrary selection
      79             :        - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
      80           8 :     PetscCall(EPSSetArbitrarySelection(eps,MyArbitrarySelection,&sxr));
      81           8 :     PetscCall(EPSSetWhichEigenpairs(eps,EPS_LARGEST_MAGNITUDE));
      82           8 :     PetscCall(EPSSolve(eps));
      83           8 :     PetscCall(EPSErrorView(eps,EPS_ERROR_RELATIVE,NULL));
      84           0 :   } else PetscCall(PetscPrintf(PETSC_COMM_WORLD,"Problem: no eigenpairs converged.\n"));
      85             : 
      86           8 :   PetscCall(EPSDestroy(&eps));
      87           8 :   PetscCall(VecDestroy(&sxr));
      88           8 :   PetscCall(VecDestroy(&sxi));
      89           8 :   PetscCall(MatDestroy(&A));
      90           8 :   PetscCall(SlepcFinalize());
      91             :   return 0;
      92             : }
      93             : 
      94             : /*TEST
      95             : 
      96             :    testset:
      97             :       args: -eps_max_it 5000 -st_pc_type jacobi
      98             :       output_file: output/test13_1.out
      99             :       filter: sed -e "s/-1.98975/-1.98974/"
     100             :       test:
     101             :          suffix: 1
     102             :          args: -eps_type {{krylovschur gd}}
     103             :       test:
     104             :          suffix: 1_jd
     105             :          args: -eps_type jd -st_ksp_type gmres
     106             :       test:
     107             :          suffix: 1_gd2
     108             :          args: -eps_type gd -eps_gd_double_expansion
     109             :       test:
     110             :          suffix: 2
     111             :          args: -eps_non_hermitian -eps_type {{krylovschur gd jd}}
     112             :       test:
     113             :          suffix: 2_gd2
     114             :          args: -eps_non_hermitian -eps_type gd -eps_gd_double_expansion
     115             :          timeoutfactor: 2
     116             : 
     117             : TEST*/