Actual source code: ex44.c
slepc-3.16.2 2022-02-01
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[] = "Compute rightmost eigenvalues with Lyapunov inverse iteration.\n\n"
12: "Loads matrix from a file or builds the same problem as ex36.c (with fixed parameters).\n\n"
13: "The command line options are:\n"
14: " -file <filename>, where <filename> = matrix file in PETSc binary form.\n"
15: " -shift <sigma>, shift to make the matrix stable.\n"
16: " -n <n>, block dimension of the 2x2 block matrix (if matrix is generated).\n\n";
18: #include <slepceps.h>
20: int main(int argc,char **argv)
21: {
22: Mat A; /* operator matrix */
23: EPS eps; /* eigenproblem solver context */
24: EPSType type;
25: PetscScalar alpha,beta,tau1,tau2,delta1,delta2,L,h,sigma=0.0;
26: PetscInt n=30,i,Istart,Iend,nev;
27: char filename[PETSC_MAX_PATH_LEN];
28: PetscViewer viewer;
29: PetscBool flg,terse;
32: SlepcInitialize(&argc,&argv,(char*)0,help);if (ierr) return ierr;
34: PetscOptionsGetString(NULL,NULL,"-file",filename,sizeof(filename),&flg);
35: if (flg) {
37: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
38: Load the matrix from file
39: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
40: PetscPrintf(PETSC_COMM_WORLD,"\nEigenproblem stored in file.\n\n");
41: #if defined(PETSC_USE_COMPLEX)
42: PetscPrintf(PETSC_COMM_WORLD," Reading COMPLEX matrix from a binary file...\n");
43: #else
44: PetscPrintf(PETSC_COMM_WORLD," Reading REAL matrix from a binary file...\n");
45: #endif
46: PetscViewerBinaryOpen(PETSC_COMM_WORLD,filename,FILE_MODE_READ,&viewer);
47: MatCreate(PETSC_COMM_WORLD,&A);
48: MatSetFromOptions(A);
49: MatLoad(A,viewer);
50: PetscViewerDestroy(&viewer);
52: } else {
54: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
55: Generate Brusselator matrix
56: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
57: PetscOptionsGetInt(NULL,NULL,"-n",&n,NULL);
58: PetscPrintf(PETSC_COMM_WORLD,"\nBrusselator wave model, n=%D\n\n",n);
60: alpha = 2.0;
61: beta = 5.45;
62: delta1 = 0.008;
63: delta2 = 0.004;
64: L = 0.51302;
66: h = 1.0 / (PetscReal)(n+1);
67: tau1 = delta1 / ((h*L)*(h*L));
68: tau2 = delta2 / ((h*L)*(h*L));
70: MatCreate(PETSC_COMM_WORLD,&A);
71: MatSetSizes(A,PETSC_DECIDE,PETSC_DECIDE,2*n,2*n);
72: MatSetFromOptions(A);
73: MatSetUp(A);
75: MatGetOwnershipRange(A,&Istart,&Iend);
76: for (i=Istart;i<Iend;i++) {
77: if (i<n) { /* upper blocks */
78: if (i>0) { MatSetValue(A,i,i-1,tau1,INSERT_VALUES); }
79: if (i<n-1) { MatSetValue(A,i,i+1,tau1,INSERT_VALUES); }
80: MatSetValue(A,i,i,-2.0*tau1+beta-1.0,INSERT_VALUES);
81: MatSetValue(A,i,i+n,alpha*alpha,INSERT_VALUES);
82: } else { /* lower blocks */
83: if (i>n) { MatSetValue(A,i,i-1,tau2,INSERT_VALUES); }
84: if (i<2*n-1) { MatSetValue(A,i,i+1,tau2,INSERT_VALUES); }
85: MatSetValue(A,i,i,-2.0*tau2-alpha*alpha,INSERT_VALUES);
86: MatSetValue(A,i,i-n,-beta,INSERT_VALUES);
87: }
88: }
89: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
90: MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
91: }
93: /* Shift the matrix to make it stable, A-sigma*I */
94: PetscOptionsGetScalar(NULL,NULL,"-shift",&sigma,NULL);
95: MatShift(A,-sigma);
97: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
98: Create the eigensolver and set various options
99: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
101: EPSCreate(PETSC_COMM_WORLD,&eps);
102: EPSSetOperators(eps,A,NULL);
103: EPSSetProblemType(eps,EPS_NHEP);
104: EPSSetWhichEigenpairs(eps,EPS_LARGEST_REAL);
105: EPSSetFromOptions(eps);
107: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
108: Solve the eigensystem
109: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
111: EPSSolve(eps);
112: EPSGetType(eps,&type);
113: PetscPrintf(PETSC_COMM_WORLD," Solution method: %s\n\n",type);
114: EPSGetDimensions(eps,&nev,NULL,NULL);
115: PetscPrintf(PETSC_COMM_WORLD," Number of requested eigenvalues: %D\n",nev);
117: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
118: Display solution and clean up
119: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
121: /* show detailed info unless -terse option is given by user */
122: PetscOptionsHasName(NULL,NULL,"-terse",&terse);
123: if (terse) {
124: EPSErrorView(eps,EPS_ERROR_RELATIVE,NULL);
125: } else {
126: PetscViewerPushFormat(PETSC_VIEWER_STDOUT_WORLD,PETSC_VIEWER_ASCII_INFO_DETAIL);
127: EPSConvergedReasonView(eps,PETSC_VIEWER_STDOUT_WORLD);
128: EPSErrorView(eps,EPS_ERROR_RELATIVE,PETSC_VIEWER_STDOUT_WORLD);
129: PetscViewerPopFormat(PETSC_VIEWER_STDOUT_WORLD);
130: }
131: EPSDestroy(&eps);
132: MatDestroy(&A);
133: SlepcFinalize();
134: return ierr;
135: }
137: /*TEST
139: testset:
140: args: -eps_nev 6 -shift 0.1 -eps_type {{krylovschur lyapii}} -eps_tol 1e-7 -terse
141: requires: double
142: filter: grep -v method | sed -e "s/-0.09981-2.13938i, -0.09981+2.13938i/-0.09981+2.13938i, -0.09981-2.13938i/" | sed -e "s/-0.77192-2.52712i, -0.77192+2.52712i/-0.77192+2.52712i, -0.77192-2.52712i/" | sed -e "s/-1.88445-3.02666i, -1.88445+3.02666i/-1.88445+3.02666i, -1.88445-3.02666i/"
143: output_file: output/ex44_1.out
144: test:
145: suffix: 1
146: test:
147: suffix: 2
148: args: -eps_lyapii_ranks 8,20 -options_left no
150: test:
151: suffix: 2_feast
152: args: -eps_type feast -eps_interval -103,-90 -terse
153: requires: feast !single
155: TEST*/