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[] = "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";
17 :
18 : #include <slepceps.h>
19 :
20 4 : int main(int argc,char **argv)
21 : {
22 4 : Mat A; /* operator matrix */
23 4 : EPS eps; /* eigenproblem solver context */
24 4 : EPSType type;
25 4 : PetscScalar alpha,beta,tau1,tau2,delta1,delta2,L,h,sigma=0.0;
26 4 : PetscInt n=30,i,Istart,Iend,nev;
27 4 : char filename[PETSC_MAX_PATH_LEN];
28 4 : PetscViewer viewer;
29 4 : PetscBool flg,terse;
30 :
31 4 : PetscFunctionBeginUser;
32 4 : PetscCall(SlepcInitialize(&argc,&argv,NULL,help));
33 :
34 4 : PetscCall(PetscOptionsGetString(NULL,NULL,"-file",filename,sizeof(filename),&flg));
35 4 : if (flg) {
36 :
37 : /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
38 : Load the matrix from file
39 : - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
40 0 : PetscCall(PetscPrintf(PETSC_COMM_WORLD,"\nEigenproblem stored in file.\n\n"));
41 : #if defined(PETSC_USE_COMPLEX)
42 0 : 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 0 : PetscCall(PetscViewerBinaryOpen(PETSC_COMM_WORLD,filename,FILE_MODE_READ,&viewer));
47 0 : PetscCall(MatCreate(PETSC_COMM_WORLD,&A));
48 0 : PetscCall(MatSetFromOptions(A));
49 0 : PetscCall(MatLoad(A,viewer));
50 0 : PetscCall(PetscViewerDestroy(&viewer));
51 :
52 : } else {
53 :
54 : /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
55 : Generate Brusselator matrix
56 : - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
57 4 : PetscCall(PetscOptionsGetInt(NULL,NULL,"-n",&n,NULL));
58 4 : PetscCall(PetscPrintf(PETSC_COMM_WORLD,"\nBrusselator wave model, n=%" PetscInt_FMT "\n\n",n));
59 :
60 4 : alpha = 2.0;
61 4 : beta = 5.45;
62 4 : delta1 = 0.008;
63 4 : delta2 = 0.004;
64 4 : L = 0.51302;
65 :
66 4 : h = 1.0 / (PetscReal)(n+1);
67 4 : tau1 = delta1 / ((h*L)*(h*L));
68 4 : tau2 = delta2 / ((h*L)*(h*L));
69 :
70 4 : PetscCall(MatCreate(PETSC_COMM_WORLD,&A));
71 4 : PetscCall(MatSetSizes(A,PETSC_DECIDE,PETSC_DECIDE,2*n,2*n));
72 4 : PetscCall(MatSetFromOptions(A));
73 :
74 4 : PetscCall(MatGetOwnershipRange(A,&Istart,&Iend));
75 244 : for (i=Istart;i<Iend;i++) {
76 240 : if (i<n) { /* upper blocks */
77 120 : if (i>0) PetscCall(MatSetValue(A,i,i-1,tau1,INSERT_VALUES));
78 120 : if (i<n-1) PetscCall(MatSetValue(A,i,i+1,tau1,INSERT_VALUES));
79 120 : PetscCall(MatSetValue(A,i,i,-2.0*tau1+beta-1.0,INSERT_VALUES));
80 120 : PetscCall(MatSetValue(A,i,i+n,alpha*alpha,INSERT_VALUES));
81 : } else { /* lower blocks */
82 120 : if (i>n) PetscCall(MatSetValue(A,i,i-1,tau2,INSERT_VALUES));
83 120 : if (i<2*n-1) PetscCall(MatSetValue(A,i,i+1,tau2,INSERT_VALUES));
84 120 : PetscCall(MatSetValue(A,i,i,-2.0*tau2-alpha*alpha,INSERT_VALUES));
85 240 : PetscCall(MatSetValue(A,i,i-n,-beta,INSERT_VALUES));
86 : }
87 : }
88 4 : PetscCall(MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY));
89 4 : PetscCall(MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY));
90 : }
91 :
92 : /* Shift the matrix to make it stable, A-sigma*I */
93 4 : PetscCall(PetscOptionsGetScalar(NULL,NULL,"-shift",&sigma,NULL));
94 4 : PetscCall(MatShift(A,-sigma));
95 :
96 : /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
97 : Create the eigensolver and set various options
98 : - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
99 :
100 4 : PetscCall(EPSCreate(PETSC_COMM_WORLD,&eps));
101 4 : PetscCall(EPSSetOperators(eps,A,NULL));
102 4 : PetscCall(EPSSetProblemType(eps,EPS_NHEP));
103 4 : PetscCall(EPSSetWhichEigenpairs(eps,EPS_LARGEST_REAL));
104 4 : PetscCall(EPSSetFromOptions(eps));
105 :
106 : /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
107 : Solve the eigensystem
108 : - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
109 :
110 4 : PetscCall(EPSSolve(eps));
111 4 : PetscCall(EPSGetType(eps,&type));
112 4 : PetscCall(PetscPrintf(PETSC_COMM_WORLD," Solution method: %s\n\n",type));
113 4 : PetscCall(EPSGetDimensions(eps,&nev,NULL,NULL));
114 4 : PetscCall(PetscPrintf(PETSC_COMM_WORLD," Number of requested eigenvalues: %" PetscInt_FMT "\n",nev));
115 :
116 : /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
117 : Display solution and clean up
118 : - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
119 :
120 : /* show detailed info unless -terse option is given by user */
121 4 : PetscCall(PetscOptionsHasName(NULL,NULL,"-terse",&terse));
122 4 : if (terse) PetscCall(EPSErrorView(eps,EPS_ERROR_RELATIVE,NULL));
123 : else {
124 0 : PetscCall(PetscViewerPushFormat(PETSC_VIEWER_STDOUT_WORLD,PETSC_VIEWER_ASCII_INFO_DETAIL));
125 0 : PetscCall(EPSConvergedReasonView(eps,PETSC_VIEWER_STDOUT_WORLD));
126 0 : PetscCall(EPSErrorView(eps,EPS_ERROR_RELATIVE,PETSC_VIEWER_STDOUT_WORLD));
127 0 : PetscCall(PetscViewerPopFormat(PETSC_VIEWER_STDOUT_WORLD));
128 : }
129 4 : PetscCall(EPSDestroy(&eps));
130 4 : PetscCall(MatDestroy(&A));
131 4 : PetscCall(SlepcFinalize());
132 : return 0;
133 : }
134 :
135 : /*TEST
136 :
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
147 :
148 : TEST*/
|
