LCOV - code coverage report
Current view: top level - pep/tutorials - ex38.c (source / functions) Hit Total Coverage
Test: SLEPc Lines: 84 90 93.3 %
Date: 2024-11-21 00:34:55 Functions: 1 1 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             : static char help[] = "Spectrum slicing on quadratic symmetric eigenproblem.\n\n"
      12             :   "The command line options are:\n"
      13             :   "  -n <n> ... dimension of the matrices.\n\n";
      14             : 
      15             : #include <slepcpep.h>
      16             : 
      17           2 : int main(int argc,char **argv)
      18             : {
      19           2 :   Mat            M,C,K,A[3]; /* problem matrices */
      20           2 :   PEP            pep;        /* polynomial eigenproblem solver context */
      21           2 :   ST             st;         /* spectral transformation context */
      22           2 :   KSP            ksp;
      23           2 :   PC             pc;
      24           2 :   PEPType        type;
      25           2 :   PetscBool      show=PETSC_FALSE,terse;
      26           2 :   PetscInt       n=100,Istart,Iend,nev,i,*inertias,ns;
      27           2 :   PetscReal      mu=1,tau=10,kappa=5,inta,intb,*shifts;
      28             : 
      29           2 :   PetscFunctionBeginUser;
      30           2 :   PetscCall(SlepcInitialize(&argc,&argv,NULL,help));
      31             : 
      32           2 :   PetscCall(PetscOptionsGetInt(NULL,NULL,"-n",&n,NULL));
      33           2 :   PetscCall(PetscOptionsGetBool(NULL,NULL,"-show_inertias",&show,NULL));
      34           2 :   PetscCall(PetscOptionsHasName(NULL,NULL,"-terse",&terse));
      35           2 :   PetscCall(PetscPrintf(PETSC_COMM_WORLD,"\nSpectrum slicing on PEP, n=%" PetscInt_FMT "\n\n",n));
      36             : 
      37             :   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
      38             :      Compute the matrices that define the eigensystem, (k^2*M+k*C+K)x=0
      39             :      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
      40             : 
      41             :   /* K is a tridiagonal */
      42           2 :   PetscCall(MatCreate(PETSC_COMM_WORLD,&K));
      43           2 :   PetscCall(MatSetSizes(K,PETSC_DECIDE,PETSC_DECIDE,n,n));
      44           2 :   PetscCall(MatSetFromOptions(K));
      45             : 
      46           2 :   PetscCall(MatGetOwnershipRange(K,&Istart,&Iend));
      47         202 :   for (i=Istart;i<Iend;i++) {
      48         200 :     if (i>0) PetscCall(MatSetValue(K,i,i-1,-kappa,INSERT_VALUES));
      49         200 :     PetscCall(MatSetValue(K,i,i,kappa*3.0,INSERT_VALUES));
      50         200 :     if (i<n-1) PetscCall(MatSetValue(K,i,i+1,-kappa,INSERT_VALUES));
      51             :   }
      52             : 
      53           2 :   PetscCall(MatAssemblyBegin(K,MAT_FINAL_ASSEMBLY));
      54           2 :   PetscCall(MatAssemblyEnd(K,MAT_FINAL_ASSEMBLY));
      55             : 
      56             :   /* C is a tridiagonal */
      57           2 :   PetscCall(MatCreate(PETSC_COMM_WORLD,&C));
      58           2 :   PetscCall(MatSetSizes(C,PETSC_DECIDE,PETSC_DECIDE,n,n));
      59           2 :   PetscCall(MatSetFromOptions(C));
      60             : 
      61           2 :   PetscCall(MatGetOwnershipRange(C,&Istart,&Iend));
      62         202 :   for (i=Istart;i<Iend;i++) {
      63         200 :     if (i>0) PetscCall(MatSetValue(C,i,i-1,-tau,INSERT_VALUES));
      64         200 :     PetscCall(MatSetValue(C,i,i,tau*3.0,INSERT_VALUES));
      65         200 :     if (i<n-1) PetscCall(MatSetValue(C,i,i+1,-tau,INSERT_VALUES));
      66             :   }
      67             : 
      68           2 :   PetscCall(MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY));
      69           2 :   PetscCall(MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY));
      70             : 
      71             :   /* M is a diagonal matrix */
      72           2 :   PetscCall(MatCreate(PETSC_COMM_WORLD,&M));
      73           2 :   PetscCall(MatSetSizes(M,PETSC_DECIDE,PETSC_DECIDE,n,n));
      74           2 :   PetscCall(MatSetFromOptions(M));
      75           2 :   PetscCall(MatGetOwnershipRange(M,&Istart,&Iend));
      76         202 :   for (i=Istart;i<Iend;i++) PetscCall(MatSetValue(M,i,i,mu,INSERT_VALUES));
      77           2 :   PetscCall(MatAssemblyBegin(M,MAT_FINAL_ASSEMBLY));
      78           2 :   PetscCall(MatAssemblyEnd(M,MAT_FINAL_ASSEMBLY));
      79             : 
      80             :   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
      81             :                 Create the eigensolver and solve the problem
      82             :      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
      83             : 
      84             :   /*
      85             :      Create eigensolver context
      86             :   */
      87           2 :   PetscCall(PEPCreate(PETSC_COMM_WORLD,&pep));
      88             : 
      89             :   /*
      90             :      Set operators and set problem type
      91             :   */
      92           2 :   A[0] = K; A[1] = C; A[2] = M;
      93           2 :   PetscCall(PEPSetOperators(pep,3,A));
      94           2 :   PetscCall(PEPSetProblemType(pep,PEP_HYPERBOLIC));
      95             : 
      96             :   /*
      97             :      Set interval for spectrum slicing
      98             :   */
      99           2 :   inta = -11.3;
     100           2 :   intb = -9.5;
     101           2 :   PetscCall(PEPSetInterval(pep,inta,intb));
     102           2 :   PetscCall(PEPSetWhichEigenpairs(pep,PEP_ALL));
     103             : 
     104             :   /*
     105             :      Spectrum slicing requires STOAR
     106             :   */
     107           2 :   PetscCall(PEPSetType(pep,PEPSTOAR));
     108             : 
     109             :   /*
     110             :      Set shift-and-invert with Cholesky; select MUMPS if available
     111             :   */
     112           2 :   PetscCall(PEPGetST(pep,&st));
     113           2 :   PetscCall(STSetType(st,STSINVERT));
     114             : 
     115           2 :   PetscCall(STGetKSP(st,&ksp));
     116           2 :   PetscCall(KSPSetType(ksp,KSPPREONLY));
     117           2 :   PetscCall(KSPGetPC(ksp,&pc));
     118           2 :   PetscCall(PCSetType(pc,PCCHOLESKY));
     119             : 
     120             :   /*
     121             :      Use MUMPS if available.
     122             :      Note that in complex scalars we cannot use MUMPS for spectrum slicing,
     123             :      because MatGetInertia() is not available in that case.
     124             :   */
     125             : #if defined(PETSC_HAVE_MUMPS) && !defined(PETSC_USE_COMPLEX)
     126             :   PetscCall(PEPSTOARSetDetectZeros(pep,PETSC_TRUE));  /* enforce zero detection */
     127             :   PetscCall(PCFactorSetMatSolverType(pc,MATSOLVERMUMPS));
     128             :   /*
     129             :      Add several MUMPS options (see ex43.c for a better way of setting them in program):
     130             :      '-st_mat_mumps_icntl_13 1': turn off ScaLAPACK for matrix inertia
     131             :      '-st_mat_mumps_icntl_24 1': detect null pivots in factorization (for the case that a shift is equal to an eigenvalue)
     132             :      '-st_mat_mumps_cntl_3 <tol>': a tolerance used for null pivot detection (must be larger than machine epsilon)
     133             : 
     134             :      Note: depending on the interval, it may be necessary also to increase the workspace:
     135             :      '-st_mat_mumps_icntl_14 <percentage>': increase workspace with a percentage (50, 100 or more)
     136             :   */
     137             :   PetscCall(PetscOptionsInsertString(NULL,"-st_mat_mumps_icntl_13 1 -st_mat_mumps_icntl_24 1 -st_mat_mumps_cntl_3 1e-12"));
     138             : #endif
     139             : 
     140             :   /*
     141             :      Set solver parameters at runtime
     142             :   */
     143           2 :   PetscCall(PEPSetFromOptions(pep));
     144             : 
     145             :   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
     146             :                       Solve the eigensystem
     147             :      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
     148           2 :   PetscCall(PEPSetUp(pep));
     149           2 :   if (show) {
     150           1 :     PetscCall(PEPSTOARGetInertias(pep,&ns,&shifts,&inertias));
     151           1 :     PetscCall(PetscPrintf(PETSC_COMM_WORLD,"Subintervals (after setup):\n"));
     152           3 :     for (i=0;i<ns;i++) PetscCall(PetscPrintf(PETSC_COMM_WORLD,"Shift %g  Inertia %" PetscInt_FMT " \n",(double)shifts[i],inertias[i]));
     153           1 :     PetscCall(PetscPrintf(PETSC_COMM_WORLD,"\n"));
     154           1 :     PetscCall(PetscFree(shifts));
     155           1 :     PetscCall(PetscFree(inertias));
     156             :   }
     157           2 :   PetscCall(PEPSolve(pep));
     158           2 :   if (show && !terse) {
     159           0 :     PetscCall(PEPSTOARGetInertias(pep,&ns,&shifts,&inertias));
     160           0 :     PetscCall(PetscPrintf(PETSC_COMM_WORLD,"All shifts (after solve):\n"));
     161           0 :     for (i=0;i<ns;i++) PetscCall(PetscPrintf(PETSC_COMM_WORLD,"Shift %g  Inertia %" PetscInt_FMT " \n",(double)shifts[i],inertias[i]));
     162           0 :     PetscCall(PetscPrintf(PETSC_COMM_WORLD,"\n"));
     163           0 :     PetscCall(PetscFree(shifts));
     164           0 :     PetscCall(PetscFree(inertias));
     165             :   }
     166             : 
     167             :   /*
     168             :      Show eigenvalues in interval and print solution
     169             :   */
     170           2 :   PetscCall(PEPGetType(pep,&type));
     171           2 :   PetscCall(PetscPrintf(PETSC_COMM_WORLD," Solution method: %s\n\n",type));
     172           2 :   PetscCall(PEPGetDimensions(pep,&nev,NULL,NULL));
     173           2 :   PetscCall(PEPGetInterval(pep,&inta,&intb));
     174           2 :   PetscCall(PetscPrintf(PETSC_COMM_WORLD," %" PetscInt_FMT " eigenvalues found in [%g, %g]\n",nev,(double)inta,(double)intb));
     175             : 
     176             :   /*
     177             :      Show detailed info unless -terse option is given by user
     178             :    */
     179           2 :   if (terse) PetscCall(PEPErrorView(pep,PEP_ERROR_BACKWARD,NULL));
     180             :   else {
     181           1 :     PetscCall(PetscViewerPushFormat(PETSC_VIEWER_STDOUT_WORLD,PETSC_VIEWER_ASCII_INFO_DETAIL));
     182           1 :     PetscCall(PEPConvergedReasonView(pep,PETSC_VIEWER_STDOUT_WORLD));
     183           1 :     PetscCall(PEPErrorView(pep,PEP_ERROR_BACKWARD,PETSC_VIEWER_STDOUT_WORLD));
     184           1 :     PetscCall(PetscViewerPopFormat(PETSC_VIEWER_STDOUT_WORLD));
     185             :   }
     186             : 
     187             :   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
     188             :                     Clean up
     189             :      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
     190           2 :   PetscCall(PEPDestroy(&pep));
     191           2 :   PetscCall(MatDestroy(&M));
     192           2 :   PetscCall(MatDestroy(&C));
     193           2 :   PetscCall(MatDestroy(&K));
     194           2 :   PetscCall(SlepcFinalize());
     195             :   return 0;
     196             : }
     197             : 
     198             : /*TEST
     199             : 
     200             :    testset:
     201             :       requires: !single
     202             :       args: -show_inertias -terse
     203             :       output_file: output/ex38_1.out
     204             :       test:
     205             :          suffix: 1
     206             :          requires: !complex
     207             :       test:
     208             :          suffix: 1_complex
     209             :          requires: complex !mumps
     210             : 
     211             :    test:
     212             :       suffix: 2
     213             :       args: -pep_interval 0,2
     214             : 
     215             : TEST*/

Generated by: LCOV version 1.14