LCOV - code coverage report
Current view: top level - eps/tutorials - ex30.c (source / functions) Hit Total Coverage
Test: SLEPc Lines: 117 136 86.0 %
Date: 2024-11-21 00:34:55 Functions: 3 4 75.0 %
Legend: Lines: hit not hit

          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[] = "Illustrates the use of a region for filtering; the number of wanted eigenvalues is not known a priori.\n\n"
      12             :   "The problem is the Brusselator wave model as in ex9.c.\n"
      13             :   "The command line options are:\n"
      14             :   "  -n <n>, where <n> = block dimension of the 2x2 block matrix.\n"
      15             :   "  -L <L>, where <L> = bifurcation parameter.\n"
      16             :   "  -alpha <alpha>, -beta <beta>, -delta1 <delta1>,  -delta2 <delta2>,\n"
      17             :   "       where <alpha> <beta> <delta1> <delta2> = model parameters.\n\n";
      18             : 
      19             : #include <slepceps.h>
      20             : 
      21             : /*
      22             :    This example tries to compute all eigenvalues lying outside the real axis.
      23             :    This could be achieved by computing LARGEST_IMAGINARY eigenvalues, but
      24             :    here we take a different route: define a region of the complex plane where
      25             :    eigenvalues must be emphasized (eigenvalues outside the region are filtered
      26             :    out). In this case, we select the region as the complement of a thin stripe
      27             :    around the real axis.
      28             :  */
      29             : 
      30             : PetscErrorCode MatMult_Brussel(Mat,Vec,Vec);
      31             : PetscErrorCode MatGetDiagonal_Brussel(Mat,Vec);
      32             : PetscErrorCode MyStoppingTest(EPS,PetscInt,PetscInt,PetscInt,PetscInt,EPSConvergedReason*,void*);
      33             : 
      34             : typedef struct {
      35             :   Mat         T;
      36             :   Vec         x1,x2,y1,y2;
      37             :   PetscScalar alpha,beta,tau1,tau2,sigma;
      38             :   PetscInt    lastnconv;      /* last value of nconv; used in stopping test */
      39             :   PetscInt    nreps;          /* number of repetitions of nconv; used in stopping test */
      40             : } CTX_BRUSSEL;
      41             : 
      42           1 : int main(int argc,char **argv)
      43             : {
      44           1 :   Mat            A;               /* eigenvalue problem matrix */
      45           1 :   EPS            eps;             /* eigenproblem solver context */
      46           1 :   RG             rg;              /* region object */
      47           1 :   PetscScalar    delta1,delta2,L,h;
      48           1 :   PetscInt       N=30,n,i,Istart,Iend,mpd;
      49           1 :   CTX_BRUSSEL    *ctx;
      50           1 :   PetscBool      terse;
      51           1 :   PetscViewer    viewer;
      52             : 
      53           1 :   PetscFunctionBeginUser;
      54           1 :   PetscCall(SlepcInitialize(&argc,&argv,NULL,help));
      55             : 
      56           1 :   PetscCall(PetscOptionsGetInt(NULL,NULL,"-n",&N,NULL));
      57           1 :   PetscCall(PetscPrintf(PETSC_COMM_WORLD,"\nBrusselator wave model, n=%" PetscInt_FMT "\n\n",N));
      58             : 
      59             :   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
      60             :         Generate the matrix
      61             :      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
      62             : 
      63             :   /*
      64             :      Create shell matrix context and set default parameters
      65             :   */
      66           1 :   PetscCall(PetscNew(&ctx));
      67           1 :   ctx->alpha = 2.0;
      68           1 :   ctx->beta  = 5.45;
      69           1 :   delta1     = 0.008;
      70           1 :   delta2     = 0.004;
      71           1 :   L          = 0.51302;
      72             : 
      73             :   /*
      74             :      Look the command line for user-provided parameters
      75             :   */
      76           1 :   PetscCall(PetscOptionsGetScalar(NULL,NULL,"-L",&L,NULL));
      77           1 :   PetscCall(PetscOptionsGetScalar(NULL,NULL,"-alpha",&ctx->alpha,NULL));
      78           1 :   PetscCall(PetscOptionsGetScalar(NULL,NULL,"-beta",&ctx->beta,NULL));
      79           1 :   PetscCall(PetscOptionsGetScalar(NULL,NULL,"-delta1",&delta1,NULL));
      80           1 :   PetscCall(PetscOptionsGetScalar(NULL,NULL,"-delta2",&delta2,NULL));
      81             : 
      82             :   /*
      83             :      Create matrix T
      84             :   */
      85           1 :   PetscCall(MatCreate(PETSC_COMM_WORLD,&ctx->T));
      86           1 :   PetscCall(MatSetSizes(ctx->T,PETSC_DECIDE,PETSC_DECIDE,N,N));
      87           1 :   PetscCall(MatSetFromOptions(ctx->T));
      88             : 
      89           1 :   PetscCall(MatGetOwnershipRange(ctx->T,&Istart,&Iend));
      90         101 :   for (i=Istart;i<Iend;i++) {
      91         100 :     if (i>0) PetscCall(MatSetValue(ctx->T,i,i-1,1.0,INSERT_VALUES));
      92         100 :     if (i<N-1) PetscCall(MatSetValue(ctx->T,i,i+1,1.0,INSERT_VALUES));
      93         100 :     PetscCall(MatSetValue(ctx->T,i,i,-2.0,INSERT_VALUES));
      94             :   }
      95           1 :   PetscCall(MatAssemblyBegin(ctx->T,MAT_FINAL_ASSEMBLY));
      96           1 :   PetscCall(MatAssemblyEnd(ctx->T,MAT_FINAL_ASSEMBLY));
      97           1 :   PetscCall(MatGetLocalSize(ctx->T,&n,NULL));
      98             : 
      99             :   /*
     100             :      Fill the remaining information in the shell matrix context
     101             :      and create auxiliary vectors
     102             :   */
     103           1 :   h = 1.0 / (PetscReal)(N+1);
     104           1 :   ctx->tau1 = delta1 / ((h*L)*(h*L));
     105           1 :   ctx->tau2 = delta2 / ((h*L)*(h*L));
     106           1 :   ctx->sigma = 0.0;
     107           1 :   PetscCall(VecCreateMPIWithArray(PETSC_COMM_WORLD,1,n,PETSC_DECIDE,NULL,&ctx->x1));
     108           1 :   PetscCall(VecCreateMPIWithArray(PETSC_COMM_WORLD,1,n,PETSC_DECIDE,NULL,&ctx->x2));
     109           1 :   PetscCall(VecCreateMPIWithArray(PETSC_COMM_WORLD,1,n,PETSC_DECIDE,NULL,&ctx->y1));
     110           1 :   PetscCall(VecCreateMPIWithArray(PETSC_COMM_WORLD,1,n,PETSC_DECIDE,NULL,&ctx->y2));
     111             : 
     112             :   /*
     113             :      Create the shell matrix
     114             :   */
     115           1 :   PetscCall(MatCreateShell(PETSC_COMM_WORLD,2*n,2*n,2*N,2*N,(void*)ctx,&A));
     116           1 :   PetscCall(MatShellSetOperation(A,MATOP_MULT,(void(*)(void))MatMult_Brussel));
     117           1 :   PetscCall(MatShellSetOperation(A,MATOP_GET_DIAGONAL,(void(*)(void))MatGetDiagonal_Brussel));
     118             : 
     119             :   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
     120             :                 Create the eigensolver and configure the region
     121             :      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
     122             : 
     123           1 :   PetscCall(EPSCreate(PETSC_COMM_WORLD,&eps));
     124           1 :   PetscCall(EPSSetOperators(eps,A,NULL));
     125           1 :   PetscCall(EPSSetProblemType(eps,EPS_NHEP));
     126             : 
     127             :   /*
     128             :      Define the region containing the eigenvalues of interest
     129             :   */
     130           1 :   PetscCall(EPSGetRG(eps,&rg));
     131           1 :   PetscCall(RGSetType(rg,RGINTERVAL));
     132           1 :   PetscCall(RGIntervalSetEndpoints(rg,-PETSC_INFINITY,PETSC_INFINITY,-0.01,0.01));
     133           1 :   PetscCall(RGSetComplement(rg,PETSC_TRUE));
     134             :   /* sort eigenvalue approximations wrt a target, otherwise convergence will be erratic */
     135           1 :   PetscCall(EPSSetTarget(eps,0.0));
     136           1 :   PetscCall(EPSSetWhichEigenpairs(eps,EPS_TARGET_MAGNITUDE));
     137             : 
     138             :   /*
     139             :      Set solver options. In particular, we must allocate sufficient
     140             :      storage for all eigenpairs that may converge (ncv). This is
     141             :      application-dependent.
     142             :   */
     143           1 :   mpd = 40;
     144           1 :   PetscCall(EPSSetDimensions(eps,2*mpd,3*mpd,mpd));
     145           1 :   PetscCall(EPSSetTolerances(eps,1e-7,2000));
     146           1 :   ctx->lastnconv = 0;
     147           1 :   ctx->nreps     = 0;
     148           1 :   PetscCall(EPSSetStoppingTestFunction(eps,MyStoppingTest,(void*)ctx,NULL));
     149           1 :   PetscCall(EPSSetFromOptions(eps));
     150             : 
     151             :   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
     152             :                 Solve the eigensystem and display solution
     153             :      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
     154             : 
     155           1 :   PetscCall(EPSSolve(eps));
     156             : 
     157             :   /* show detailed info unless -terse option is given by user */
     158           1 :   PetscCall(PetscViewerASCIIGetStdout(PETSC_COMM_WORLD,&viewer));
     159           1 :   PetscCall(PetscViewerPushFormat(viewer,PETSC_VIEWER_ASCII_INFO_DETAIL));
     160           1 :   PetscCall(EPSConvergedReasonView(eps,viewer));
     161           1 :   PetscCall(PetscOptionsHasName(NULL,NULL,"-terse",&terse));
     162           1 :   if (!terse) PetscCall(EPSErrorView(eps,EPS_ERROR_RELATIVE,viewer));
     163           1 :   PetscCall(PetscViewerPopFormat(viewer));
     164             : 
     165           1 :   PetscCall(EPSDestroy(&eps));
     166           1 :   PetscCall(MatDestroy(&A));
     167           1 :   PetscCall(MatDestroy(&ctx->T));
     168           1 :   PetscCall(VecDestroy(&ctx->x1));
     169           1 :   PetscCall(VecDestroy(&ctx->x2));
     170           1 :   PetscCall(VecDestroy(&ctx->y1));
     171           1 :   PetscCall(VecDestroy(&ctx->y2));
     172           1 :   PetscCall(PetscFree(ctx));
     173           1 :   PetscCall(SlepcFinalize());
     174             :   return 0;
     175             : }
     176             : 
     177         550 : PetscErrorCode MatMult_Brussel(Mat A,Vec x,Vec y)
     178             : {
     179         550 :   PetscInt          n;
     180         550 :   const PetscScalar *px;
     181         550 :   PetscScalar       *py;
     182         550 :   CTX_BRUSSEL       *ctx;
     183             : 
     184         550 :   PetscFunctionBeginUser;
     185         550 :   PetscCall(MatShellGetContext(A,&ctx));
     186         550 :   PetscCall(MatGetLocalSize(ctx->T,&n,NULL));
     187         550 :   PetscCall(VecGetArrayRead(x,&px));
     188         550 :   PetscCall(VecGetArray(y,&py));
     189         550 :   PetscCall(VecPlaceArray(ctx->x1,px));
     190         550 :   PetscCall(VecPlaceArray(ctx->x2,px+n));
     191         550 :   PetscCall(VecPlaceArray(ctx->y1,py));
     192         550 :   PetscCall(VecPlaceArray(ctx->y2,py+n));
     193             : 
     194         550 :   PetscCall(MatMult(ctx->T,ctx->x1,ctx->y1));
     195         550 :   PetscCall(VecScale(ctx->y1,ctx->tau1));
     196         550 :   PetscCall(VecAXPY(ctx->y1,ctx->beta - 1.0 + ctx->sigma,ctx->x1));
     197         550 :   PetscCall(VecAXPY(ctx->y1,ctx->alpha * ctx->alpha,ctx->x2));
     198             : 
     199         550 :   PetscCall(MatMult(ctx->T,ctx->x2,ctx->y2));
     200         550 :   PetscCall(VecScale(ctx->y2,ctx->tau2));
     201         550 :   PetscCall(VecAXPY(ctx->y2,-ctx->beta,ctx->x1));
     202         550 :   PetscCall(VecAXPY(ctx->y2,-ctx->alpha * ctx->alpha + ctx->sigma,ctx->x2));
     203             : 
     204         550 :   PetscCall(VecRestoreArrayRead(x,&px));
     205         550 :   PetscCall(VecRestoreArray(y,&py));
     206         550 :   PetscCall(VecResetArray(ctx->x1));
     207         550 :   PetscCall(VecResetArray(ctx->x2));
     208         550 :   PetscCall(VecResetArray(ctx->y1));
     209         550 :   PetscCall(VecResetArray(ctx->y2));
     210         550 :   PetscFunctionReturn(PETSC_SUCCESS);
     211             : }
     212             : 
     213           0 : PetscErrorCode MatGetDiagonal_Brussel(Mat A,Vec diag)
     214             : {
     215           0 :   Vec            d1,d2;
     216           0 :   PetscInt       n;
     217           0 :   PetscScalar    *pd;
     218           0 :   MPI_Comm       comm;
     219           0 :   CTX_BRUSSEL    *ctx;
     220             : 
     221           0 :   PetscFunctionBeginUser;
     222           0 :   PetscCall(MatShellGetContext(A,&ctx));
     223           0 :   PetscCall(PetscObjectGetComm((PetscObject)A,&comm));
     224           0 :   PetscCall(MatGetLocalSize(ctx->T,&n,NULL));
     225           0 :   PetscCall(VecGetArray(diag,&pd));
     226           0 :   PetscCall(VecCreateMPIWithArray(comm,1,n,PETSC_DECIDE,pd,&d1));
     227           0 :   PetscCall(VecCreateMPIWithArray(comm,1,n,PETSC_DECIDE,pd+n,&d2));
     228             : 
     229           0 :   PetscCall(VecSet(d1,-2.0*ctx->tau1 + ctx->beta - 1.0 + ctx->sigma));
     230           0 :   PetscCall(VecSet(d2,-2.0*ctx->tau2 - ctx->alpha*ctx->alpha + ctx->sigma));
     231             : 
     232           0 :   PetscCall(VecDestroy(&d1));
     233           0 :   PetscCall(VecDestroy(&d2));
     234           0 :   PetscCall(VecRestoreArray(diag,&pd));
     235           0 :   PetscFunctionReturn(PETSC_SUCCESS);
     236             : }
     237             : 
     238             : /*
     239             :     Function for user-defined stopping test.
     240             : 
     241             :     Ignores the value of nev. It only takes into account the number of
     242             :     eigenpairs that have converged in recent outer iterations (restarts);
     243             :     if no new eigenvalues have converged in the last few restarts,
     244             :     we stop the iteration, assuming that no more eigenvalues are present
     245             :     inside the region.
     246             : */
     247          26 : PetscErrorCode MyStoppingTest(EPS eps,PetscInt its,PetscInt max_it,PetscInt nconv,PetscInt nev,EPSConvergedReason *reason,void *ptr)
     248             : {
     249          26 :   CTX_BRUSSEL    *ctx = (CTX_BRUSSEL*)ptr;
     250             : 
     251          26 :   PetscFunctionBeginUser;
     252             :   /* check usual termination conditions, but ignoring the case nconv>=nev */
     253          26 :   PetscCall(EPSStoppingBasic(eps,its,max_it,nconv,PETSC_INT_MAX,reason,NULL));
     254          26 :   if (*reason==EPS_CONVERGED_ITERATING) {
     255             :     /* check if nconv is the same as before */
     256          26 :     if (nconv==ctx->lastnconv) ctx->nreps++;
     257             :     else {
     258           4 :       ctx->lastnconv = nconv;
     259           4 :       ctx->nreps     = 0;
     260             :     }
     261             :     /* check if no eigenvalues converged in last 10 restarts */
     262          26 :     if (nconv && ctx->nreps>10) *reason = EPS_CONVERGED_USER;
     263             :   }
     264          26 :   PetscFunctionReturn(PETSC_SUCCESS);
     265             : }
     266             : 
     267             : /*TEST
     268             : 
     269             :    test:
     270             :       suffix: 1
     271             :       args: -n 100 -terse
     272             :       requires: !single
     273             : 
     274             : TEST*/

Generated by: LCOV version 1.14