Actual source code: test37.c
slepc-3.16.0 2021-09-30
1: /*
2: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
3: SLEPc - Scalable Library for Eigenvalue Problem Computations
4: Copyright (c) 2002-2021, 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;
44: SlepcInitialize(&argc,&argv,(char*)0,help);if (ierr) return ierr;
45: PetscOptionsGetInt(NULL,NULL,"-n",&N,NULL);
46: PetscPrintf(PETSC_COMM_WORLD,"\nBrusselator wave model, n=%D\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: PetscOptionsGetScalar(NULL,NULL,"-L",&L,NULL);
59: PetscOptionsGetScalar(NULL,NULL,"-alpha",&alpha,NULL);
60: PetscOptionsGetScalar(NULL,NULL,"-beta",&beta,NULL);
61: PetscOptionsGetScalar(NULL,NULL,"-delta1",&delta1,NULL);
62: 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: MatCreate(PETSC_COMM_WORLD,&T1);
70: MatSetSizes(T1,PETSC_DECIDE,PETSC_DECIDE,N,N);
71: MatSetFromOptions(T1);
72: MatSetUp(T1);
73: MatGetOwnershipRange(T1,&Istart,&Iend);
74: for (i=Istart;i<Iend;i++) {
75: if (i>0) { MatSetValue(T1,i,i-1,1.0,INSERT_VALUES); }
76: if (i<N-1) { MatSetValue(T1,i,i+1,1.0,INSERT_VALUES); }
77: MatSetValue(T1,i,i,-2.0,INSERT_VALUES);
78: }
79: MatAssemblyBegin(T1,MAT_FINAL_ASSEMBLY);
80: MatAssemblyEnd(T1,MAT_FINAL_ASSEMBLY);
82: MatDuplicate(T1,MAT_COPY_VALUES,&T2);
83: MatScale(T1,tau1);
84: MatShift(T1,beta-1.0);
85: MatScale(T2,tau2);
86: MatShift(T2,-alpha*alpha);
88: /* Create matrices D1, D2 */
89: MatCreate(PETSC_COMM_WORLD,&D1);
90: MatSetSizes(D1,PETSC_DECIDE,PETSC_DECIDE,N,N);
91: MatSetFromOptions(D1);
92: MatSetUp(D1);
93: MatAssemblyBegin(D1,MAT_FINAL_ASSEMBLY);
94: MatAssemblyEnd(D1,MAT_FINAL_ASSEMBLY);
95: MatDuplicate(D1,MAT_COPY_VALUES,&D2);
96: MatShift(D1,alpha*alpha);
97: MatShift(D2,-beta);
99: /* Create the nest matrix */
100: mats[0] = T1;
101: mats[1] = D1;
102: mats[2] = D2;
103: mats[3] = T2;
104: MatCreateNest(PETSC_COMM_WORLD,2,NULL,2,NULL,mats,&A);
106: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
107: Create the eigensolver and solve the problem
108: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
110: EPSCreate(PETSC_COMM_WORLD,&eps);
111: EPSSetOperators(eps,A,NULL);
112: EPSSetProblemType(eps,EPS_NHEP);
113: EPSSetWhichEigenpairs(eps,EPS_LARGEST_REAL);
114: EPSSetTrueResidual(eps,PETSC_FALSE);
115: EPSSetFromOptions(eps);
116: EPSSolve(eps);
118: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
119: Display solution and clean up
120: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
122: EPSErrorView(eps,EPS_ERROR_RELATIVE,NULL);
123: EPSDestroy(&eps);
124: MatDestroy(&A);
125: MatDestroy(&T1);
126: MatDestroy(&T2);
127: MatDestroy(&D1);
128: MatDestroy(&D2);
129: SlepcFinalize();
130: return ierr;
131: }
133: /*TEST
135: test:
136: suffix: 1
137: args: -eps_nev 4
138: requires: !single
139: 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/"
141: TEST*/