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[] = "Tests solving an eigenproblem defined with MatNest. "
      12             :   "Based on ex9.\n"
      13             :   "The command line options are:\n"
      14             :   "  -n <n>, where <n> = block dimension of the 2x2 block matrix.\n"
      15             :   "  -L <L>, where <L> = bifurcation parameter.\n"
      16             :   "  -alpha <alpha>, -beta <beta>, -delta1 <delta1>,  -delta2 <delta2>,\n"
      17             :   "       where <alpha> <beta> <delta1> <delta2> = model parameters.\n\n";
      18             : 
      19             : #include <slepceps.h>
      20             : 
      21             : /*
      22             :    This example computes the eigenvalues with largest real part of the
      23             :    following matrix
      24             : 
      25             :         A = [ tau1*T+(beta-1)*I     alpha^2*I
      26             :                   -beta*I        tau2*T-alpha^2*I ],
      27             : 
      28             :    where
      29             : 
      30             :         T = tridiag{1,-2,1}
      31             :         h = 1/(n+1)
      32             :         tau1 = delta1/(h*L)^2
      33             :         tau2 = delta2/(h*L)^2
      34             :  */
      35             : 
      36           1 : int main(int argc,char **argv)
      37             : {
      38           1 :   EPS            eps;
      39           1 :   Mat            A,T1,T2,D1,D2,mats[4];
      40           1 :   PetscScalar    alpha,beta,tau1,tau2,delta1,delta2,L,h;
      41           1 :   PetscInt       N=30,i,Istart,Iend;
      42             : 
      43           1 :   PetscFunctionBeginUser;
      44           1 :   PetscCall(SlepcInitialize(&argc,&argv,NULL,help));
      45           1 :   PetscCall(PetscOptionsGetInt(NULL,NULL,"-n",&N,NULL));
      46           1 :   PetscCall(PetscPrintf(PETSC_COMM_WORLD,"\nBrusselator wave model, n=%" PetscInt_FMT "\n\n",N));
      47             : 
      48             :   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
      49             :         Generate the matrix
      50             :      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
      51             : 
      52           1 :   alpha  = 2.0;
      53           1 :   beta   = 5.45;
      54           1 :   delta1 = 0.008;
      55           1 :   delta2 = 0.004;
      56           1 :   L      = 0.51302;
      57             : 
      58           1 :   PetscCall(PetscOptionsGetScalar(NULL,NULL,"-L",&L,NULL));
      59           1 :   PetscCall(PetscOptionsGetScalar(NULL,NULL,"-alpha",&alpha,NULL));
      60           1 :   PetscCall(PetscOptionsGetScalar(NULL,NULL,"-beta",&beta,NULL));
      61           1 :   PetscCall(PetscOptionsGetScalar(NULL,NULL,"-delta1",&delta1,NULL));
      62           1 :   PetscCall(PetscOptionsGetScalar(NULL,NULL,"-delta2",&delta2,NULL));
      63             : 
      64           1 :   h    = 1.0 / (PetscReal)(N+1);
      65           1 :   tau1 = delta1 / ((h*L)*(h*L));
      66           1 :   tau2 = delta2 / ((h*L)*(h*L));
      67             : 
      68             :   /* Create matrices T1, T2 */
      69           1 :   PetscCall(MatCreate(PETSC_COMM_WORLD,&T1));
      70           1 :   PetscCall(MatSetSizes(T1,PETSC_DECIDE,PETSC_DECIDE,N,N));
      71           1 :   PetscCall(MatSetFromOptions(T1));
      72           1 :   PetscCall(MatGetOwnershipRange(T1,&Istart,&Iend));
      73          31 :   for (i=Istart;i<Iend;i++) {
      74          30 :     if (i>0) PetscCall(MatSetValue(T1,i,i-1,1.0,INSERT_VALUES));
      75          30 :     if (i<N-1) PetscCall(MatSetValue(T1,i,i+1,1.0,INSERT_VALUES));
      76          30 :     PetscCall(MatSetValue(T1,i,i,-2.0,INSERT_VALUES));
      77             :   }
      78           1 :   PetscCall(MatAssemblyBegin(T1,MAT_FINAL_ASSEMBLY));
      79           1 :   PetscCall(MatAssemblyEnd(T1,MAT_FINAL_ASSEMBLY));
      80             : 
      81           1 :   PetscCall(MatDuplicate(T1,MAT_COPY_VALUES,&T2));
      82           1 :   PetscCall(MatScale(T1,tau1));
      83           1 :   PetscCall(MatShift(T1,beta-1.0));
      84           1 :   PetscCall(MatScale(T2,tau2));
      85           1 :   PetscCall(MatShift(T2,-alpha*alpha));
      86             : 
      87             :   /* Create matrices D1, D2 */
      88           1 :   PetscCall(MatCreateConstantDiagonal(PETSC_COMM_WORLD,PETSC_DECIDE,PETSC_DECIDE,N,N,alpha*alpha,&D1));
      89           1 :   PetscCall(MatCreateConstantDiagonal(PETSC_COMM_WORLD,PETSC_DECIDE,PETSC_DECIDE,N,N,-beta,&D2));
      90             : 
      91             :   /* Create the nest matrix */
      92           1 :   mats[0] = T1;
      93           1 :   mats[1] = D1;
      94           1 :   mats[2] = D2;
      95           1 :   mats[3] = T2;
      96           1 :   PetscCall(MatCreateNest(PETSC_COMM_WORLD,2,NULL,2,NULL,mats,&A));
      97             : 
      98             :   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
      99             :                 Create the eigensolver and solve the problem
     100             :      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
     101             : 
     102           1 :   PetscCall(EPSCreate(PETSC_COMM_WORLD,&eps));
     103           1 :   PetscCall(EPSSetOperators(eps,A,NULL));
     104           1 :   PetscCall(EPSSetProblemType(eps,EPS_NHEP));
     105           1 :   PetscCall(EPSSetWhichEigenpairs(eps,EPS_LARGEST_REAL));
     106           1 :   PetscCall(EPSSetTrueResidual(eps,PETSC_FALSE));
     107           1 :   PetscCall(EPSSetFromOptions(eps));
     108           1 :   PetscCall(EPSSolve(eps));
     109             : 
     110             :   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
     111             :                     Display solution and clean up
     112             :      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
     113             : 
     114           1 :   PetscCall(EPSErrorView(eps,EPS_ERROR_RELATIVE,NULL));
     115           1 :   PetscCall(EPSDestroy(&eps));
     116           1 :   PetscCall(MatDestroy(&A));
     117           1 :   PetscCall(MatDestroy(&T1));
     118           1 :   PetscCall(MatDestroy(&T2));
     119           1 :   PetscCall(MatDestroy(&D1));
     120           1 :   PetscCall(MatDestroy(&D2));
     121           1 :   PetscCall(SlepcFinalize());
     122             :   return 0;
     123             : }
     124             : 
     125             : /*TEST
     126             : 
     127             :    test:
     128             :       suffix: 1
     129             :       args: -eps_nev 4
     130             :       requires: !single
     131             :       filter: sed -e "s/0.00019-2.13938i, 0.00019+2.13938i/0.00019+2.13938i, 0.00019-2.13938i/" | sed -e "s/-0.67192-2.52712i, -0.67192+2.52712i/-0.67192+2.52712i, -0.67192-2.52712i/"
     132             : 
     133             : TEST*/