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[] = "Block-diagonal orthogonal eigenproblem.\n"
      12             :   "The command line options are:\n"
      13             :   "  -n <n>, where <n> = matrix dimension.\n"
      14             :   "  -seed <s>, where <s> = seed for random number generation.\n\n";
      15             : 
      16             : #include <slepceps.h>
      17             : 
      18           2 : int main(int argc,char **argv)
      19             : {
      20           2 :   EPS            eps;
      21           2 :   Mat            A;
      22           2 :   PetscRandom    rand;
      23           2 :   PetscScalar    val,c,s;
      24           2 :   PetscInt       n=30,i,seed=0x12345678;
      25           2 :   PetscMPIInt    rank;
      26             : 
      27           2 :   PetscFunctionBeginUser;
      28           2 :   PetscCall(SlepcInitialize(&argc,&argv,(char*)0,help));
      29           2 :   PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD,&rank));
      30           2 :   PetscCall(PetscOptionsGetInt(NULL,NULL,"-n",&n,NULL));
      31           2 :   PetscCall(PetscPrintf(PETSC_COMM_WORLD,"\nOrthogonal eigenproblem, n=%" PetscInt_FMT "\n\n",n));
      32             : 
      33             :   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
      34             :         Generate the matrix
      35             :      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
      36             : 
      37           2 :   PetscCall(PetscRandomCreate(PETSC_COMM_WORLD,&rand));
      38           2 :   PetscCall(PetscRandomSetFromOptions(rand));
      39           2 :   PetscCall(PetscOptionsGetInt(NULL,NULL,"-seed",&seed,NULL));
      40           2 :   PetscCall(PetscRandomSetSeed(rand,seed));
      41           2 :   PetscCall(PetscRandomSeed(rand));
      42           2 :   PetscCall(PetscRandomSetInterval(rand,0,2*PETSC_PI));
      43             : 
      44           2 :   PetscCall(MatCreate(PETSC_COMM_WORLD,&A));
      45           2 :   PetscCall(MatSetSizes(A,PETSC_DECIDE,PETSC_DECIDE,n,n));
      46           2 :   PetscCall(MatSetFromOptions(A));
      47             : 
      48           2 :   if (!rank) {
      49          32 :     for (i=0;i<n/2;i++) {
      50          30 :       PetscCall(PetscRandomGetValue(rand,&val));
      51          30 :       c = PetscCosReal(PetscRealPart(val));
      52          30 :       s = PetscSinReal(PetscRealPart(val));
      53          30 :       PetscCall(MatSetValue(A,2*i,2*i,c,INSERT_VALUES));
      54          30 :       PetscCall(MatSetValue(A,2*i+1,2*i+1,c,INSERT_VALUES));
      55          30 :       PetscCall(MatSetValue(A,2*i,2*i+1,s,INSERT_VALUES));
      56          30 :       PetscCall(MatSetValue(A,2*i+1,2*i,-s,INSERT_VALUES));
      57             :     }
      58           2 :     if (n%2) PetscCall(MatSetValue(A,n-1,n-1,-1.0,INSERT_VALUES));
      59             :   }
      60           2 :   PetscCall(MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY));
      61           2 :   PetscCall(MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY));
      62             : 
      63             :   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
      64             :                 Create the eigensolver and solve the problem
      65             :      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
      66             : 
      67           2 :   PetscCall(EPSCreate(PETSC_COMM_WORLD,&eps));
      68           2 :   PetscCall(EPSSetOperators(eps,A,NULL));
      69           2 :   PetscCall(EPSSetProblemType(eps,EPS_NHEP));
      70           2 :   PetscCall(EPSSetWhichEigenpairs(eps,EPS_LARGEST_REAL));
      71           2 :   PetscCall(EPSSetFromOptions(eps));
      72           2 :   PetscCall(EPSSolve(eps));
      73             : 
      74             :   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
      75             :                     Display solution and clean up
      76             :      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
      77             : 
      78           2 :   PetscCall(EPSErrorView(eps,EPS_ERROR_RELATIVE,NULL));
      79           2 :   PetscCall(EPSDestroy(&eps));
      80           2 :   PetscCall(MatDestroy(&A));
      81           2 :   PetscCall(PetscRandomDestroy(&rand));
      82           2 :   PetscCall(SlepcFinalize());
      83             :   return 0;
      84             : }
      85             : 
      86             : /*TEST
      87             : 
      88             :    testset:
      89             :       requires: complex double
      90             :       args: -eps_type ciss -eps_all -rg_type ring -rg_ring_center 0 -rg_ring_radius 1 -rg_ring_width 0.05 -rg_ring_startangle .93 -rg_ring_endangle .07
      91             :       filter: sed -e "s/[+-]\([0-9]\.[0-9]*i\)/+-\\1/g"
      92             :       test:
      93             :          suffix: 1_ring
      94             :          args: -eps_ciss_extraction {{ritz hankel}}
      95             : 
      96             : TEST*/