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 interface to external package BLOPEX.\n\n"
      12             :   "This is based on ex12.c. The command line options are:\n"
      13             :   "  -n <n>, where <n> = number of grid subdivisions in x dimension.\n"
      14             :   "  -m <m>, where <m> = number of grid subdivisions in y dimension.\n\n";
      15             : 
      16             : #include <slepceps.h>
      17             : 
      18           1 : int main(int argc,char **argv)
      19             : {
      20           1 :   Mat            A,B;         /* matrices */
      21           1 :   EPS            eps;         /* eigenproblem solver context */
      22           1 :   ST             st;          /* spectral transformation context */
      23           1 :   KSP            ksp;
      24           1 :   PC             pc;
      25           1 :   PetscInt       N,n=35,m,Istart,Iend,II,i,j,bs;
      26           1 :   PetscBool      flag;
      27             : 
      28           1 :   PetscFunctionBeginUser;
      29           1 :   PetscCall(SlepcInitialize(&argc,&argv,(char*)0,help));
      30             : 
      31           1 :   PetscCall(PetscOptionsGetInt(NULL,NULL,"-n",&n,NULL));
      32           1 :   PetscCall(PetscOptionsGetInt(NULL,NULL,"-m",&m,&flag));
      33           1 :   if (!flag) m=n;
      34           1 :   N = n*m;
      35           1 :   PetscCall(PetscPrintf(PETSC_COMM_WORLD,"\nGeneralized eigenproblem with BLOPEX, N=%" PetscInt_FMT " (%" PetscInt_FMT "x%" PetscInt_FMT " grid)\n\n",N,n,m));
      36             : 
      37             :   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
      38             :      Compute the matrices that define the eigensystem, Ax=kBx
      39             :      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
      40             : 
      41           1 :   PetscCall(MatCreate(PETSC_COMM_WORLD,&A));
      42           1 :   PetscCall(MatSetSizes(A,PETSC_DECIDE,PETSC_DECIDE,N,N));
      43           1 :   PetscCall(MatSetFromOptions(A));
      44             : 
      45           1 :   PetscCall(MatCreate(PETSC_COMM_WORLD,&B));
      46           1 :   PetscCall(MatSetSizes(B,PETSC_DECIDE,PETSC_DECIDE,N,N));
      47           1 :   PetscCall(MatSetFromOptions(B));
      48             : 
      49           1 :   PetscCall(MatGetOwnershipRange(A,&Istart,&Iend));
      50        1226 :   for (II=Istart;II<Iend;II++) {
      51        1225 :     i = II/n; j = II-i*n;
      52        1225 :     if (i>0) PetscCall(MatSetValue(A,II,II-n,-1.0,INSERT_VALUES));
      53        1225 :     if (i<m-1) PetscCall(MatSetValue(A,II,II+n,-1.0,INSERT_VALUES));
      54        1225 :     if (j>0) PetscCall(MatSetValue(A,II,II-1,-1.0,INSERT_VALUES));
      55        1225 :     if (j<n-1) PetscCall(MatSetValue(A,II,II+1,-1.0,INSERT_VALUES));
      56        1225 :     PetscCall(MatSetValue(A,II,II,4.0,INSERT_VALUES));
      57        1225 :     PetscCall(MatSetValue(B,II,II,2.0,INSERT_VALUES));
      58             :   }
      59           1 :   if (Istart==0) {
      60           1 :     PetscCall(MatSetValue(B,0,0,6.0,INSERT_VALUES));
      61           1 :     PetscCall(MatSetValue(B,0,1,-1.0,INSERT_VALUES));
      62           1 :     PetscCall(MatSetValue(B,1,0,-1.0,INSERT_VALUES));
      63           1 :     PetscCall(MatSetValue(B,1,1,1.0,INSERT_VALUES));
      64             :   }
      65             : 
      66           1 :   PetscCall(MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY));
      67           1 :   PetscCall(MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY));
      68           1 :   PetscCall(MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY));
      69           1 :   PetscCall(MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY));
      70             : 
      71             :   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
      72             :                 Create the eigensolver and set various options
      73             :      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
      74             : 
      75           1 :   PetscCall(EPSCreate(PETSC_COMM_WORLD,&eps));
      76           1 :   PetscCall(EPSSetOperators(eps,A,B));
      77           1 :   PetscCall(EPSSetProblemType(eps,EPS_GHEP));
      78           1 :   PetscCall(EPSSetType(eps,EPSBLOPEX));
      79             : 
      80             :   /*
      81             :      Set several options
      82             :   */
      83           1 :   PetscCall(EPSSetWhichEigenpairs(eps,EPS_SMALLEST_REAL));
      84           1 :   PetscCall(EPSGetST(eps,&st));
      85           1 :   PetscCall(STSetType(st,STPRECOND));
      86           1 :   PetscCall(STGetKSP(st,&ksp));
      87           1 :   PetscCall(KSPGetPC(ksp,&pc));
      88           1 :   PetscCall(KSPSetType(ksp,KSPPREONLY));
      89           1 :   PetscCall(PCSetType(pc,PCICC));
      90             : 
      91           1 :   PetscCall(EPSBLOPEXSetBlockSize(eps,4));
      92           1 :   PetscCall(EPSSetFromOptions(eps));
      93             : 
      94             :   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
      95             :            Compute all eigenvalues in interval and display info
      96             :      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
      97             : 
      98           1 :   PetscCall(EPSSolve(eps));
      99           1 :   PetscCall(EPSBLOPEXGetBlockSize(eps,&bs));
     100           1 :   PetscCall(PetscPrintf(PETSC_COMM_WORLD," BLOPEX: using block size %" PetscInt_FMT "\n",bs));
     101             : 
     102           1 :   PetscCall(EPSErrorView(eps,EPS_ERROR_ABSOLUTE,NULL));
     103             : 
     104           1 :   PetscCall(EPSDestroy(&eps));
     105           1 :   PetscCall(MatDestroy(&A));
     106           1 :   PetscCall(MatDestroy(&B));
     107           1 :   PetscCall(SlepcFinalize());
     108             :   return 0;
     109             : }
     110             : 
     111             : /*TEST
     112             : 
     113             :    build:
     114             :       requires: blopex
     115             : 
     116             :    test:
     117             :       suffix: 1
     118             :       args: -eps_nev 8 -eps_view
     119             :       requires: blopex
     120             :       filter: grep -v tolerance | sed -e "s/hermitian/symmetric/"
     121             : 
     122             : TEST*/