Actual source code: ex44.c
slepc-3.21.0 2024-03-30
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[] = "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;
31: PetscFunctionBeginUser;
32: PetscCall(SlepcInitialize(&argc,&argv,(char*)0,help));
34: PetscCall(PetscOptionsGetString(NULL,NULL,"-file",filename,sizeof(filename),&flg));
35: if (flg) {
37: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
38: Load the matrix from file
39: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
40: PetscCall(PetscPrintf(PETSC_COMM_WORLD,"\nEigenproblem stored in file.\n\n"));
41: #if defined(PETSC_USE_COMPLEX)
42: PetscCall(PetscPrintf(PETSC_COMM_WORLD," Reading COMPLEX matrix from a binary file...\n"));
43: #else
44: PetscCall(PetscPrintf(PETSC_COMM_WORLD," Reading REAL matrix from a binary file...\n"));
45: #endif
46: PetscCall(PetscViewerBinaryOpen(PETSC_COMM_WORLD,filename,FILE_MODE_READ,&viewer));
47: PetscCall(MatCreate(PETSC_COMM_WORLD,&A));
48: PetscCall(MatSetFromOptions(A));
49: PetscCall(MatLoad(A,viewer));
50: PetscCall(PetscViewerDestroy(&viewer));
52: } else {
54: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
55: Generate Brusselator matrix
56: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
57: PetscCall(PetscOptionsGetInt(NULL,NULL,"-n",&n,NULL));
58: PetscCall(PetscPrintf(PETSC_COMM_WORLD,"\nBrusselator wave model, n=%" PetscInt_FMT "\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: PetscCall(MatCreate(PETSC_COMM_WORLD,&A));
71: PetscCall(MatSetSizes(A,PETSC_DECIDE,PETSC_DECIDE,2*n,2*n));
72: PetscCall(MatSetFromOptions(A));
74: PetscCall(MatGetOwnershipRange(A,&Istart,&Iend));
75: for (i=Istart;i<Iend;i++) {
76: if (i<n) { /* upper blocks */
77: if (i>0) PetscCall(MatSetValue(A,i,i-1,tau1,INSERT_VALUES));
78: if (i<n-1) PetscCall(MatSetValue(A,i,i+1,tau1,INSERT_VALUES));
79: PetscCall(MatSetValue(A,i,i,-2.0*tau1+beta-1.0,INSERT_VALUES));
80: PetscCall(MatSetValue(A,i,i+n,alpha*alpha,INSERT_VALUES));
81: } else { /* lower blocks */
82: if (i>n) PetscCall(MatSetValue(A,i,i-1,tau2,INSERT_VALUES));
83: if (i<2*n-1) PetscCall(MatSetValue(A,i,i+1,tau2,INSERT_VALUES));
84: PetscCall(MatSetValue(A,i,i,-2.0*tau2-alpha*alpha,INSERT_VALUES));
85: PetscCall(MatSetValue(A,i,i-n,-beta,INSERT_VALUES));
86: }
87: }
88: PetscCall(MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY));
89: PetscCall(MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY));
90: }
92: /* Shift the matrix to make it stable, A-sigma*I */
93: PetscCall(PetscOptionsGetScalar(NULL,NULL,"-shift",&sigma,NULL));
94: PetscCall(MatShift(A,-sigma));
96: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
97: Create the eigensolver and set various options
98: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
100: PetscCall(EPSCreate(PETSC_COMM_WORLD,&eps));
101: PetscCall(EPSSetOperators(eps,A,NULL));
102: PetscCall(EPSSetProblemType(eps,EPS_NHEP));
103: PetscCall(EPSSetWhichEigenpairs(eps,EPS_LARGEST_REAL));
104: PetscCall(EPSSetFromOptions(eps));
106: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
107: Solve the eigensystem
108: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
110: PetscCall(EPSSolve(eps));
111: PetscCall(EPSGetType(eps,&type));
112: PetscCall(PetscPrintf(PETSC_COMM_WORLD," Solution method: %s\n\n",type));
113: PetscCall(EPSGetDimensions(eps,&nev,NULL,NULL));
114: PetscCall(PetscPrintf(PETSC_COMM_WORLD," Number of requested eigenvalues: %" PetscInt_FMT "\n",nev));
116: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
117: Display solution and clean up
118: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
120: /* show detailed info unless -terse option is given by user */
121: PetscCall(PetscOptionsHasName(NULL,NULL,"-terse",&terse));
122: if (terse) PetscCall(EPSErrorView(eps,EPS_ERROR_RELATIVE,NULL));
123: else {
124: PetscCall(PetscViewerPushFormat(PETSC_VIEWER_STDOUT_WORLD,PETSC_VIEWER_ASCII_INFO_DETAIL));
125: PetscCall(EPSConvergedReasonView(eps,PETSC_VIEWER_STDOUT_WORLD));
126: PetscCall(EPSErrorView(eps,EPS_ERROR_RELATIVE,PETSC_VIEWER_STDOUT_WORLD));
127: PetscCall(PetscViewerPopFormat(PETSC_VIEWER_STDOUT_WORLD));
128: }
129: PetscCall(EPSDestroy(&eps));
130: PetscCall(MatDestroy(&A));
131: PetscCall(SlepcFinalize());
132: return 0;
133: }
135: /*TEST
137: testset:
138: args: -eps_nev 6 -shift 0.1 -eps_type {{krylovschur lyapii}} -eps_tol 1e-7 -terse
139: requires: double
140: 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/"
141: output_file: output/ex44_1.out
142: test:
143: suffix: 1
144: test:
145: suffix: 2
146: args: -eps_lyapii_ranks 8,20 -options_left no
148: TEST*/