1: /*
  2:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
  3:    SLEPc - Scalable Library for Eigenvalue Problem Computations
  4:    Copyright (c) 2002-, Universitat Politecnica de Valencia, Spain

  6:    This file is part of SLEPc.
  7:    SLEPc is distributed under a 2-clause BSD license (see LICENSE).
  8:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
  9: */

 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";

 19: #include <slepceps.h>

 21: /*
 22:    This example computes the eigenvalues with largest real part of the
 23:    following matrix

 25:         A = [ tau1*T+(beta-1)*I     alpha^2*I
 26:                   -beta*I        tau2*T-alpha^2*I ],

 28:    where

 30:         T = tridiag{1,-2,1}
 31:         h = 1/(n+1)
 32:         tau1 = delta1/(h*L)^2
 33:         tau2 = delta2/(h*L)^2
 34:  */

 36: int main(int argc,char **argv)
 37: {
 38:   EPS            eps;
 39:   Mat            A,T1,T2,D1,D2,mats[4];
 40:   PetscScalar    alpha,beta,tau1,tau2,delta1,delta2,L,h;
 41:   PetscInt       N=30,i,Istart,Iend;

 43:   PetscFunctionBeginUser;
 44:   PetscCall(SlepcInitialize(&argc,&argv,NULL,help));
 45:   PetscCall(PetscOptionsGetInt(NULL,NULL,"-n",&N,NULL));
 46:   PetscCall(PetscPrintf(PETSC_COMM_WORLD,"\nBrusselator wave model, n=%" PetscInt_FMT "\n\n",N));

 48:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 49:         Generate the matrix
 50:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

 52:   alpha  = 2.0;
 53:   beta   = 5.45;
 54:   delta1 = 0.008;
 55:   delta2 = 0.004;
 56:   L      = 0.51302;

 58:   PetscCall(PetscOptionsGetScalar(NULL,NULL,"-L",&L,NULL));
 59:   PetscCall(PetscOptionsGetScalar(NULL,NULL,"-alpha",&alpha,NULL));
 60:   PetscCall(PetscOptionsGetScalar(NULL,NULL,"-beta",&beta,NULL));
 61:   PetscCall(PetscOptionsGetScalar(NULL,NULL,"-delta1",&delta1,NULL));
 62:   PetscCall(PetscOptionsGetScalar(NULL,NULL,"-delta2",&delta2,NULL));

 64:   h    = 1.0 / (PetscReal)(N+1);
 65:   tau1 = delta1 / ((h*L)*(h*L));
 66:   tau2 = delta2 / ((h*L)*(h*L));

 68:   /* Create matrices T1, T2 */
 69:   PetscCall(MatCreate(PETSC_COMM_WORLD,&T1));
 70:   PetscCall(MatSetSizes(T1,PETSC_DECIDE,PETSC_DECIDE,N,N));
 71:   PetscCall(MatSetFromOptions(T1));
 72:   PetscCall(MatGetOwnershipRange(T1,&Istart,&Iend));
 73:   for (i=Istart;i<Iend;i++) {
 74:     if (i>0) PetscCall(MatSetValue(T1,i,i-1,1.0,INSERT_VALUES));
 75:     if (i<N-1) PetscCall(MatSetValue(T1,i,i+1,1.0,INSERT_VALUES));
 76:     PetscCall(MatSetValue(T1,i,i,-2.0,INSERT_VALUES));
 77:   }
 78:   PetscCall(MatAssemblyBegin(T1,MAT_FINAL_ASSEMBLY));
 79:   PetscCall(MatAssemblyEnd(T1,MAT_FINAL_ASSEMBLY));

 81:   PetscCall(MatDuplicate(T1,MAT_COPY_VALUES,&T2));
 82:   PetscCall(MatScale(T1,tau1));
 83:   PetscCall(MatShift(T1,beta-1.0));
 84:   PetscCall(MatScale(T2,tau2));
 85:   PetscCall(MatShift(T2,-alpha*alpha));

 87:   /* Create matrices D1, D2 */
 88:   PetscCall(MatCreateConstantDiagonal(PETSC_COMM_WORLD,PETSC_DECIDE,PETSC_DECIDE,N,N,alpha*alpha,&D1));
 89:   PetscCall(MatCreateConstantDiagonal(PETSC_COMM_WORLD,PETSC_DECIDE,PETSC_DECIDE,N,N,-beta,&D2));

 91:   /* Create the nest matrix */
 92:   mats[0] = T1;
 93:   mats[1] = D1;
 94:   mats[2] = D2;
 95:   mats[3] = T2;
 96:   PetscCall(MatCreateNest(PETSC_COMM_WORLD,2,NULL,2,NULL,mats,&A));

 98:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 99:                 Create the eigensolver and solve the problem
100:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

102:   PetscCall(EPSCreate(PETSC_COMM_WORLD,&eps));
103:   PetscCall(EPSSetOperators(eps,A,NULL));
104:   PetscCall(EPSSetProblemType(eps,EPS_NHEP));
105:   PetscCall(EPSSetWhichEigenpairs(eps,EPS_LARGEST_REAL));
106:   PetscCall(EPSSetTrueResidual(eps,PETSC_FALSE));
107:   PetscCall(EPSSetFromOptions(eps));
108:   PetscCall(EPSSolve(eps));

110:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
111:                     Display solution and clean up
112:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

114:   PetscCall(EPSErrorView(eps,EPS_ERROR_RELATIVE,NULL));
115:   PetscCall(EPSDestroy(&eps));
116:   PetscCall(MatDestroy(&A));
117:   PetscCall(MatDestroy(&T1));
118:   PetscCall(MatDestroy(&T2));
119:   PetscCall(MatDestroy(&D1));
120:   PetscCall(MatDestroy(&D2));
121:   PetscCall(SlepcFinalize());
122:   return 0;
123: }

125: /*TEST

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/"

133: TEST*/