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[] = "Computes the smallest nonzero eigenvalue of the Laplacian of a graph.\n\n"
      12             :   "This example illustrates EPSSetDeflationSpace(). The example graph corresponds to a "
      13             :   "2-D regular mesh. The command line options are:\n"
      14             :   "  -n <n>, where <n> = number of grid subdivisions in x dimension.\n"
      15             :   "  -m <m>, where <m> = number of grid subdivisions in y dimension.\n\n";
      16             : 
      17             : #include <slepceps.h>
      18             : 
      19           8 : int main (int argc,char **argv)
      20             : {
      21           8 :   EPS            eps;             /* eigenproblem solver context */
      22           8 :   Mat            A;               /* operator matrix */
      23           8 :   Vec            x;
      24           8 :   PetscInt       N,n=10,m,i,j,II,Istart,Iend,nev;
      25           8 :   PetscScalar    w;
      26           8 :   PetscBool      flag;
      27             : 
      28           8 :   PetscFunctionBeginUser;
      29           8 :   PetscCall(SlepcInitialize(&argc,&argv,NULL,help));
      30             : 
      31           8 :   PetscCall(PetscOptionsGetInt(NULL,NULL,"-n",&n,NULL));
      32           8 :   PetscCall(PetscOptionsGetInt(NULL,NULL,"-m",&m,&flag));
      33           8 :   if (!flag) m=n;
      34           8 :   N = n*m;
      35           8 :   PetscCall(PetscPrintf(PETSC_COMM_WORLD,"\nFiedler vector of a 2-D regular mesh, N=%" PetscInt_FMT " (%" PetscInt_FMT "x%" PetscInt_FMT " grid)\n\n",N,n,m));
      36             : 
      37             :   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
      38             :      Compute the operator matrix that defines the eigensystem, Ax=kx
      39             :      In this example, A = L(G), where L is the Laplacian of graph G, i.e.
      40             :      Lii = degree of node i, Lij = -1 if edge (i,j) exists in G
      41             :      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
      42             : 
      43           8 :   PetscCall(MatCreate(PETSC_COMM_WORLD,&A));
      44           8 :   PetscCall(MatSetSizes(A,PETSC_DECIDE,PETSC_DECIDE,N,N));
      45           8 :   PetscCall(MatSetFromOptions(A));
      46             : 
      47           8 :   PetscCall(MatGetOwnershipRange(A,&Istart,&Iend));
      48         888 :   for (II=Istart;II<Iend;II++) {
      49         880 :     i = II/n; j = II-i*n;
      50         880 :     w = 0.0;
      51         880 :     if (i>0) { PetscCall(MatSetValue(A,II,II-n,-1.0,INSERT_VALUES)); w=w+1.0; }
      52         880 :     if (i<m-1) { PetscCall(MatSetValue(A,II,II+n,-1.0,INSERT_VALUES)); w=w+1.0; }
      53         880 :     if (j>0) { PetscCall(MatSetValue(A,II,II-1,-1.0,INSERT_VALUES)); w=w+1.0; }
      54         880 :     if (j<n-1) { PetscCall(MatSetValue(A,II,II+1,-1.0,INSERT_VALUES)); w=w+1.0; }
      55         880 :     PetscCall(MatSetValue(A,II,II,w,INSERT_VALUES));
      56             :   }
      57             : 
      58           8 :   PetscCall(MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY));
      59           8 :   PetscCall(MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY));
      60             : 
      61             :   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
      62             :                 Create the eigensolver and set various options
      63             :      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
      64             : 
      65             :   /*
      66             :      Create eigensolver context
      67             :   */
      68           8 :   PetscCall(EPSCreate(PETSC_COMM_WORLD,&eps));
      69             : 
      70             :   /*
      71             :      Set operators. In this case, it is a standard eigenvalue problem
      72             :   */
      73           8 :   PetscCall(EPSSetOperators(eps,A,NULL));
      74           8 :   PetscCall(EPSSetProblemType(eps,EPS_HEP));
      75             : 
      76             :   /*
      77             :      Select portion of spectrum
      78             :   */
      79           8 :   PetscCall(EPSSetWhichEigenpairs(eps,EPS_SMALLEST_REAL));
      80             : 
      81             :   /*
      82             :      Set solver parameters at runtime
      83             :   */
      84           8 :   PetscCall(EPSSetFromOptions(eps));
      85             : 
      86             :   /*
      87             :      Attach deflation space: in this case, the matrix has a constant
      88             :      nullspace, [1 1 ... 1]^T is the eigenvector of the zero eigenvalue
      89             :   */
      90           8 :   PetscCall(MatCreateVecs(A,&x,NULL));
      91           8 :   PetscCall(VecSet(x,1.0));
      92           8 :   PetscCall(EPSSetDeflationSpace(eps,1,&x));
      93           8 :   PetscCall(VecDestroy(&x));
      94             : 
      95             :   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
      96             :                       Solve the eigensystem
      97             :      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
      98             : 
      99           8 :   PetscCall(EPSSolve(eps));
     100           8 :   PetscCall(EPSGetDimensions(eps,&nev,NULL,NULL));
     101           8 :   PetscCall(PetscPrintf(PETSC_COMM_WORLD," Number of requested eigenvalues: %" PetscInt_FMT "\n",nev));
     102             : 
     103             :   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
     104             :                     Display solution and clean up
     105             :      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
     106             : 
     107           8 :   PetscCall(EPSErrorView(eps,EPS_ERROR_RELATIVE,NULL));
     108           8 :   PetscCall(EPSDestroy(&eps));
     109           8 :   PetscCall(MatDestroy(&A));
     110           8 :   PetscCall(SlepcFinalize());
     111             :   return 0;
     112             : }
     113             : 
     114             : /*TEST
     115             : 
     116             :    testset:
     117             :       args: -eps_nev 4 -m 11 -eps_max_it 500
     118             :       output_file: output/test10_1.out
     119             :       test:
     120             :          suffix: 1
     121             :          args: -eps_type {{krylovschur arnoldi gd jd rqcg}}
     122             :       test:
     123             :          suffix: 1_lobpcg
     124             :          args: -eps_type lobpcg -eps_lobpcg_blocksize 6
     125             :       test:
     126             :          suffix: 1_lanczos
     127             :          args: -eps_type lanczos -eps_lanczos_reorthog local
     128             :          requires: !single
     129             :       test:
     130             :          suffix: 1_gd2
     131             :          args: -eps_type gd -eps_gd_double_expansion
     132             : 
     133             : TEST*/