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[] = "Solve a quadratic problem with PEPLINEAR with a user-provided EPS.\n\n"
12 : "The command line options are:\n"
13 : " -n <n>, where <n> = number of grid subdivisions in x dimension.\n"
14 : " -m <m>, where <m> = number of grid subdivisions in y dimension.\n\n";
15 :
16 : #include <slepcpep.h>
17 :
18 1 : int main(int argc,char **argv)
19 : {
20 1 : Mat M,C,K,A[3];
21 1 : PEP pep;
22 1 : PetscInt N,n=10,m,Istart,Iend,II,i,j;
23 1 : PetscBool flag,expmat;
24 1 : PetscReal alpha,beta;
25 1 : EPS eps;
26 1 : ST st;
27 1 : KSP ksp;
28 1 : PC pc;
29 :
30 1 : PetscFunctionBeginUser;
31 1 : PetscCall(SlepcInitialize(&argc,&argv,NULL,help));
32 1 : PetscCall(PetscOptionsGetInt(NULL,NULL,"-n",&n,NULL));
33 1 : PetscCall(PetscOptionsGetInt(NULL,NULL,"-m",&m,&flag));
34 1 : if (!flag) m=n;
35 1 : N = n*m;
36 1 : PetscCall(PetscPrintf(PETSC_COMM_WORLD,"\nQuadratic Eigenproblem, N=%" PetscInt_FMT " (%" PetscInt_FMT "x%" PetscInt_FMT " grid)\n\n",N,n,m));
37 :
38 : /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
39 : Compute the matrices that define the eigensystem, (k^2*M+k*C+K)x=0
40 : - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
41 :
42 : /* K is the 2-D Laplacian */
43 1 : PetscCall(MatCreate(PETSC_COMM_WORLD,&K));
44 1 : PetscCall(MatSetSizes(K,PETSC_DECIDE,PETSC_DECIDE,N,N));
45 1 : PetscCall(MatSetFromOptions(K));
46 1 : PetscCall(MatGetOwnershipRange(K,&Istart,&Iend));
47 101 : for (II=Istart;II<Iend;II++) {
48 100 : i = II/n; j = II-i*n;
49 100 : if (i>0) PetscCall(MatSetValue(K,II,II-n,-1.0,INSERT_VALUES));
50 100 : if (i<m-1) PetscCall(MatSetValue(K,II,II+n,-1.0,INSERT_VALUES));
51 100 : if (j>0) PetscCall(MatSetValue(K,II,II-1,-1.0,INSERT_VALUES));
52 100 : if (j<n-1) PetscCall(MatSetValue(K,II,II+1,-1.0,INSERT_VALUES));
53 100 : PetscCall(MatSetValue(K,II,II,4.0,INSERT_VALUES));
54 : }
55 1 : PetscCall(MatAssemblyBegin(K,MAT_FINAL_ASSEMBLY));
56 1 : PetscCall(MatAssemblyEnd(K,MAT_FINAL_ASSEMBLY));
57 :
58 : /* C is the 1-D Laplacian on horizontal lines */
59 1 : PetscCall(MatCreate(PETSC_COMM_WORLD,&C));
60 1 : PetscCall(MatSetSizes(C,PETSC_DECIDE,PETSC_DECIDE,N,N));
61 1 : PetscCall(MatSetFromOptions(C));
62 1 : PetscCall(MatGetOwnershipRange(C,&Istart,&Iend));
63 101 : for (II=Istart;II<Iend;II++) {
64 100 : i = II/n; j = II-i*n;
65 100 : if (j>0) PetscCall(MatSetValue(C,II,II-1,-1.0,INSERT_VALUES));
66 100 : if (j<n-1) PetscCall(MatSetValue(C,II,II+1,-1.0,INSERT_VALUES));
67 100 : PetscCall(MatSetValue(C,II,II,2.0,INSERT_VALUES));
68 : }
69 1 : PetscCall(MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY));
70 1 : PetscCall(MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY));
71 :
72 : /* M is a diagonal matrix */
73 1 : PetscCall(MatCreate(PETSC_COMM_WORLD,&M));
74 1 : PetscCall(MatSetSizes(M,PETSC_DECIDE,PETSC_DECIDE,N,N));
75 1 : PetscCall(MatSetFromOptions(M));
76 1 : PetscCall(MatGetOwnershipRange(M,&Istart,&Iend));
77 101 : for (II=Istart;II<Iend;II++) PetscCall(MatSetValue(M,II,II,(PetscReal)(II+1),INSERT_VALUES));
78 1 : PetscCall(MatAssemblyBegin(M,MAT_FINAL_ASSEMBLY));
79 1 : PetscCall(MatAssemblyEnd(M,MAT_FINAL_ASSEMBLY));
80 :
81 : /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
82 : Create a standalone EPS with appropriate settings
83 : - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
84 :
85 1 : PetscCall(EPSCreate(PETSC_COMM_WORLD,&eps));
86 1 : PetscCall(EPSSetWhichEigenpairs(eps,EPS_TARGET_MAGNITUDE));
87 : #if defined(PETSC_USE_COMPLEX)
88 : PetscCall(EPSSetTarget(eps,0.01*PETSC_i));
89 : #endif
90 1 : PetscCall(EPSGetST(eps,&st));
91 1 : PetscCall(STSetType(st,STSINVERT));
92 1 : PetscCall(STGetKSP(st,&ksp));
93 1 : PetscCall(KSPSetType(ksp,KSPBCGS));
94 1 : PetscCall(KSPGetPC(ksp,&pc));
95 1 : PetscCall(PCSetType(pc,PCJACOBI));
96 1 : PetscCall(EPSSetFromOptions(eps));
97 :
98 : /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
99 : Create the eigensolver and solve the eigensystem
100 : - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
101 :
102 1 : PetscCall(PEPCreate(PETSC_COMM_WORLD,&pep));
103 1 : PetscCall(PetscObjectSetName((PetscObject)pep,"PEP_solver"));
104 1 : A[0] = K; A[1] = C; A[2] = M;
105 1 : PetscCall(PEPSetOperators(pep,3,A));
106 1 : PetscCall(PEPSetType(pep,PEPLINEAR));
107 1 : PetscCall(PEPSetProblemType(pep,PEP_GENERAL));
108 1 : PetscCall(PEPLinearSetEPS(pep,eps));
109 1 : PetscCall(PEPSetFromOptions(pep));
110 1 : PetscCall(PEPSolve(pep));
111 1 : PetscCall(PEPLinearGetLinearization(pep,&alpha,&beta));
112 1 : PetscCall(PetscPrintf(PETSC_COMM_WORLD," Linearization with alpha=%g, beta=%g",(double)alpha,(double)beta));
113 1 : PetscCall(PEPLinearGetExplicitMatrix(pep,&expmat));
114 1 : if (expmat) PetscCall(PetscPrintf(PETSC_COMM_WORLD," with explicit matrix"));
115 1 : PetscCall(PetscPrintf(PETSC_COMM_WORLD,"\n"));
116 :
117 : /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
118 : Display solution and clean up
119 : - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
120 :
121 1 : PetscCall(PEPErrorView(pep,PEP_ERROR_BACKWARD,NULL));
122 1 : PetscCall(PEPDestroy(&pep));
123 1 : PetscCall(EPSDestroy(&eps));
124 1 : PetscCall(MatDestroy(&M));
125 1 : PetscCall(MatDestroy(&C));
126 1 : PetscCall(MatDestroy(&K));
127 1 : PetscCall(SlepcFinalize());
128 : return 0;
129 : }
130 :
131 : /*TEST
132 :
133 : testset:
134 : args: -pep_linear_explicitmatrix -pep_view_vectors ::ascii_info
135 : test:
136 : suffix: 1_real
137 : requires: !single !complex
138 : test:
139 : suffix: 1
140 : requires: !single complex
141 :
142 : TEST*/
|
