LCOV - code coverage report
Current view: top level - pep/tutorials/nlevp - pdde_stability.c (source / functions) Hit Total Coverage
Test: SLEPc Lines: 60 64 93.8 %
Date: 2024-12-18 00:51:33 Functions: 2 2 100.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             :    This example implements one of the problems found at
      12             :        NLEVP: A Collection of Nonlinear Eigenvalue Problems,
      13             :        The University of Manchester.
      14             :    The details of the collection can be found at:
      15             :        [1] T. Betcke et al., "NLEVP: A Collection of Nonlinear Eigenvalue
      16             :            Problems", ACM Trans. Math. Software 39(2), Article 7, 2013.
      17             : 
      18             :    The pdde_stability problem is a complex-symmetric QEP from the stability
      19             :    analysis of a discretized partial delay-differential equation. It requires
      20             :    complex scalars.
      21             : */
      22             : 
      23             : static char help[] = "Stability analysis of a discretized partial delay-differential equation.\n\n"
      24             :   "The command line options are:\n"
      25             :   "  -m <m>, grid size, the matrices have dimension n=m*m.\n"
      26             :   "  -c <a0,b0,a1,b1,a2,b2,phi1>, comma-separated list of 7 real parameters.\n\n";
      27             : 
      28             : #include <slepcpep.h>
      29             : 
      30             : #define NMAT 3
      31             : 
      32             : /*
      33             :     Function for user-defined eigenvalue ordering criterion.
      34             : 
      35             :     Given two eigenvalues ar+i*ai and br+i*bi, the subroutine must choose
      36             :     one of them as the preferred one according to the criterion.
      37             :     In this example, the preferred value is the one with absolute value closest to 1.
      38             : */
      39       16630 : PetscErrorCode MyEigenSort(PetscScalar ar,PetscScalar ai,PetscScalar br,PetscScalar bi,PetscInt *r,void *ctx)
      40             : {
      41       16630 :   PetscReal aa,ab;
      42             : 
      43       16630 :   PetscFunctionBeginUser;
      44       16630 :   aa = PetscAbsReal(SlepcAbsEigenvalue(ar,ai)-PetscRealConstant(1.0));
      45       16630 :   ab = PetscAbsReal(SlepcAbsEigenvalue(br,bi)-PetscRealConstant(1.0));
      46       16630 :   *r = aa > ab ? 1 : (aa < ab ? -1 : 0);
      47       16630 :   PetscFunctionReturn(PETSC_SUCCESS);
      48             : }
      49             : 
      50           3 : int main(int argc,char **argv)
      51             : {
      52           3 :   Mat            A[NMAT];         /* problem matrices */
      53           3 :   PEP            pep;             /* polynomial eigenproblem solver context */
      54           3 :   PetscInt       m=15,n,II,Istart,Iend,i,j,k;
      55           3 :   PetscReal      h,xi,xj,c[7] = { 2, .3, -2, .2, -2, -.3, -PETSC_PI/2 };
      56           3 :   PetscScalar    alpha,beta,gamma;
      57           3 :   PetscBool      flg,terse;
      58             : 
      59           3 :   PetscFunctionBeginUser;
      60           3 :   PetscCall(SlepcInitialize(&argc,&argv,NULL,help));
      61             : #if !defined(PETSC_USE_COMPLEX)
      62             :   SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_SUP,"This example requires complex scalars");
      63             : #endif
      64             : 
      65           3 :   PetscCall(PetscOptionsGetInt(NULL,NULL,"-m",&m,NULL));
      66           3 :   n = m*m;
      67           3 :   h = PETSC_PI/(m+1);
      68           3 :   gamma = PetscExpScalar(PETSC_i*c[6]);
      69           3 :   gamma = gamma/PetscAbsScalar(gamma);
      70           3 :   k = 7;
      71           3 :   PetscCall(PetscOptionsGetRealArray(NULL,NULL,"-c",c,&k,&flg));
      72           3 :   PetscCheck(!flg || k==7,PETSC_COMM_WORLD,PETSC_ERR_USER,"The number of parameters -c should be 7, you provided %" PetscInt_FMT,k);
      73           3 :   PetscCall(PetscPrintf(PETSC_COMM_WORLD,"\nPDDE stability, n=%" PetscInt_FMT " (m=%" PetscInt_FMT ")\n\n",n,m));
      74             : 
      75             :   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
      76             :                      Compute the polynomial matrices
      77             :      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
      78             : 
      79             :   /* initialize matrices */
      80          12 :   for (i=0;i<NMAT;i++) {
      81           9 :     PetscCall(MatCreate(PETSC_COMM_WORLD,&A[i]));
      82           9 :     PetscCall(MatSetSizes(A[i],PETSC_DECIDE,PETSC_DECIDE,n,n));
      83           9 :     PetscCall(MatSetFromOptions(A[i]));
      84             :   }
      85           3 :   PetscCall(MatGetOwnershipRange(A[0],&Istart,&Iend));
      86             : 
      87             :   /* A[1] has a pattern similar to the 2D Laplacian */
      88         678 :   for (II=Istart;II<Iend;II++) {
      89         675 :     i = II/m; j = II-i*m;
      90         675 :     xi = (i+1)*h; xj = (j+1)*h;
      91         675 :     alpha = c[0]+c[1]*PetscSinReal(xi)+gamma*(c[2]+c[3]*xi*(1.0-PetscExpReal(xi-PETSC_PI)));
      92         675 :     beta = c[0]+c[1]*PetscSinReal(xj)-gamma*(c[2]+c[3]*xj*(1.0-PetscExpReal(xj-PETSC_PI)));
      93         675 :     PetscCall(MatSetValue(A[1],II,II,alpha+beta-4.0/(h*h),INSERT_VALUES));
      94         675 :     if (j>0) PetscCall(MatSetValue(A[1],II,II-1,1.0/(h*h),INSERT_VALUES));
      95         675 :     if (j<m-1) PetscCall(MatSetValue(A[1],II,II+1,1.0/(h*h),INSERT_VALUES));
      96         675 :     if (i>0) PetscCall(MatSetValue(A[1],II,II-m,1.0/(h*h),INSERT_VALUES));
      97         675 :     if (i<m-1) PetscCall(MatSetValue(A[1],II,II+m,1.0/(h*h),INSERT_VALUES));
      98             :   }
      99             : 
     100             :   /* A[0] and A[2] are diagonal */
     101         678 :   for (II=Istart;II<Iend;II++) {
     102         675 :     i = II/m; j = II-i*m;
     103         675 :     xi = (i+1)*h; xj = (j+1)*h;
     104         675 :     alpha = c[4]+c[5]*xi*(PETSC_PI-xi);
     105         675 :     beta = c[4]+c[5]*xj*(PETSC_PI-xj);
     106         675 :     PetscCall(MatSetValue(A[0],II,II,alpha,INSERT_VALUES));
     107         675 :     PetscCall(MatSetValue(A[2],II,II,beta,INSERT_VALUES));
     108             :   }
     109             : 
     110             :   /* assemble matrices */
     111          12 :   for (i=0;i<NMAT;i++) PetscCall(MatAssemblyBegin(A[i],MAT_FINAL_ASSEMBLY));
     112          12 :   for (i=0;i<NMAT;i++) PetscCall(MatAssemblyEnd(A[i],MAT_FINAL_ASSEMBLY));
     113             : 
     114             :   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
     115             :                 Create the eigensolver and solve the problem
     116             :      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
     117             : 
     118           3 :   PetscCall(PEPCreate(PETSC_COMM_WORLD,&pep));
     119           3 :   PetscCall(PEPSetOperators(pep,NMAT,A));
     120           3 :   PetscCall(PEPSetEigenvalueComparison(pep,MyEigenSort,NULL));
     121           3 :   PetscCall(PEPSetDimensions(pep,4,PETSC_DETERMINE,PETSC_DETERMINE));
     122           3 :   PetscCall(PEPSetFromOptions(pep));
     123           3 :   PetscCall(PEPSolve(pep));
     124             : 
     125             :   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
     126             :                     Display solution and clean up
     127             :      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
     128             : 
     129             :   /* show detailed info unless -terse option is given by user */
     130           3 :   PetscCall(PetscOptionsHasName(NULL,NULL,"-terse",&terse));
     131           3 :   if (terse) PetscCall(PEPErrorView(pep,PEP_ERROR_BACKWARD,NULL));
     132             :   else {
     133           0 :     PetscCall(PetscViewerPushFormat(PETSC_VIEWER_STDOUT_WORLD,PETSC_VIEWER_ASCII_INFO_DETAIL));
     134           0 :     PetscCall(PEPConvergedReasonView(pep,PETSC_VIEWER_STDOUT_WORLD));
     135           0 :     PetscCall(PEPErrorView(pep,PEP_ERROR_BACKWARD,PETSC_VIEWER_STDOUT_WORLD));
     136           0 :     PetscCall(PetscViewerPopFormat(PETSC_VIEWER_STDOUT_WORLD));
     137             :   }
     138           3 :   PetscCall(PEPDestroy(&pep));
     139          12 :   for (i=0;i<NMAT;i++) PetscCall(MatDestroy(&A[i]));
     140           3 :   PetscCall(SlepcFinalize());
     141             :   return 0;
     142             : }
     143             : 
     144             : /*TEST
     145             : 
     146             :    build:
     147             :       requires: complex
     148             : 
     149             :    test:
     150             :       suffix: 1
     151             :       args: -pep_type {{toar qarnoldi linear}} -pep_ncv 25 -terse
     152             :       requires: complex double
     153             : 
     154             : TEST*/

Generated by: LCOV version 1.14