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 CISS.\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 : RG rg;
23 1 : KSP *ksp;
24 1 : PC pc;
25 1 : PEPCISSExtraction ext;
26 1 : PetscInt N,n=10,m,Istart,Iend,II,i,j,nsolve;
27 1 : PetscBool flg;
28 :
29 1 : PetscFunctionBeginUser;
30 1 : PetscCall(SlepcInitialize(&argc,&argv,NULL,help));
31 1 : PetscCall(PetscOptionsGetInt(NULL,NULL,"-n",&n,NULL));
32 1 : PetscCall(PetscOptionsGetInt(NULL,NULL,"-m",&m,&flg));
33 1 : if (!flg) m=n;
34 1 : N = n*m;
35 1 : PetscCall(PetscPrintf(PETSC_COMM_WORLD,"\nQuadratic Eigenproblem, N=%" PetscInt_FMT " (%" PetscInt_FMT "x%" PetscInt_FMT " grid)\n\n",N,n,m));
36 :
37 : /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
38 : Compute the matrices that define the eigensystem, (k^2*M+k*C+K)x=0
39 : - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
40 :
41 : /* K is the 2-D Laplacian */
42 1 : PetscCall(MatCreate(PETSC_COMM_WORLD,&K));
43 1 : PetscCall(MatSetSizes(K,PETSC_DECIDE,PETSC_DECIDE,N,N));
44 1 : PetscCall(MatSetFromOptions(K));
45 1 : PetscCall(MatGetOwnershipRange(K,&Istart,&Iend));
46 101 : for (II=Istart;II<Iend;II++) {
47 100 : i = II/n; j = II-i*n;
48 100 : if (i>0) PetscCall(MatSetValue(K,II,II-n,-1.0,INSERT_VALUES));
49 100 : if (i<m-1) PetscCall(MatSetValue(K,II,II+n,-1.0,INSERT_VALUES));
50 100 : if (j>0) PetscCall(MatSetValue(K,II,II-1,-1.0,INSERT_VALUES));
51 100 : if (j<n-1) PetscCall(MatSetValue(K,II,II+1,-1.0,INSERT_VALUES));
52 100 : PetscCall(MatSetValue(K,II,II,4.0,INSERT_VALUES));
53 : }
54 1 : PetscCall(MatAssemblyBegin(K,MAT_FINAL_ASSEMBLY));
55 1 : PetscCall(MatAssemblyEnd(K,MAT_FINAL_ASSEMBLY));
56 :
57 : /* C is the 1-D Laplacian on horizontal lines */
58 1 : PetscCall(MatCreate(PETSC_COMM_WORLD,&C));
59 1 : PetscCall(MatSetSizes(C,PETSC_DECIDE,PETSC_DECIDE,N,N));
60 1 : PetscCall(MatSetFromOptions(C));
61 1 : PetscCall(MatGetOwnershipRange(C,&Istart,&Iend));
62 101 : for (II=Istart;II<Iend;II++) {
63 100 : i = II/n; j = II-i*n;
64 100 : if (j>0) PetscCall(MatSetValue(C,II,II-1,-1.0,INSERT_VALUES));
65 100 : if (j<n-1) PetscCall(MatSetValue(C,II,II+1,-1.0,INSERT_VALUES));
66 100 : PetscCall(MatSetValue(C,II,II,2.0,INSERT_VALUES));
67 : }
68 1 : PetscCall(MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY));
69 1 : PetscCall(MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY));
70 :
71 : /* M is a diagonal matrix */
72 1 : PetscCall(MatCreate(PETSC_COMM_WORLD,&M));
73 1 : PetscCall(MatSetSizes(M,PETSC_DECIDE,PETSC_DECIDE,N,N));
74 1 : PetscCall(MatSetFromOptions(M));
75 1 : PetscCall(MatGetOwnershipRange(M,&Istart,&Iend));
76 101 : for (II=Istart;II<Iend;II++) PetscCall(MatSetValue(M,II,II,(PetscReal)(II+1),INSERT_VALUES));
77 1 : PetscCall(MatAssemblyBegin(M,MAT_FINAL_ASSEMBLY));
78 1 : PetscCall(MatAssemblyEnd(M,MAT_FINAL_ASSEMBLY));
79 :
80 : /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
81 : Create the eigensolver and solve the eigensystem
82 : - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
83 :
84 1 : PetscCall(PEPCreate(PETSC_COMM_WORLD,&pep));
85 1 : A[0] = K; A[1] = C; A[2] = M;
86 1 : PetscCall(PEPSetOperators(pep,3,A));
87 1 : PetscCall(PEPSetProblemType(pep,PEP_GENERAL));
88 :
89 : /* customize polynomial eigensolver; set runtime options */
90 1 : PetscCall(PEPSetType(pep,PEPCISS));
91 1 : PetscCall(PEPGetRG(pep,&rg));
92 1 : PetscCall(RGSetType(rg,RGELLIPSE));
93 1 : PetscCall(RGEllipseSetParameters(rg,PetscCMPLX(-0.1,0.3),0.1,0.25));
94 1 : PetscCall(PEPCISSSetSizes(pep,24,PETSC_DETERMINE,PETSC_DETERMINE,1,PETSC_DETERMINE,PETSC_TRUE));
95 1 : PetscCall(PEPCISSGetKSPs(pep,&nsolve,&ksp));
96 25 : for (i=0;i<nsolve;i++) {
97 24 : PetscCall(KSPSetTolerances(ksp[i],1e-12,PETSC_CURRENT,PETSC_CURRENT,PETSC_CURRENT));
98 24 : PetscCall(KSPSetType(ksp[i],KSPPREONLY));
99 24 : PetscCall(KSPGetPC(ksp[i],&pc));
100 24 : PetscCall(PCSetType(pc,PCLU));
101 : }
102 1 : PetscCall(PEPSetFromOptions(pep));
103 :
104 : /* solve */
105 1 : PetscCall(PetscObjectTypeCompare((PetscObject)pep,PEPCISS,&flg));
106 1 : if (flg) {
107 1 : PetscCall(PEPCISSGetExtraction(pep,&ext));
108 1 : PetscCall(PetscPrintf(PETSC_COMM_WORLD," Running CISS with %" PetscInt_FMT " KSP solvers (%s extraction)\n",nsolve,PEPCISSExtractions[ext]));
109 : }
110 1 : PetscCall(PEPSolve(pep));
111 :
112 : /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
113 : Display solution and clean up
114 : - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
115 :
116 1 : PetscCall(PEPErrorView(pep,PEP_ERROR_BACKWARD,NULL));
117 1 : PetscCall(PEPDestroy(&pep));
118 1 : PetscCall(MatDestroy(&M));
119 1 : PetscCall(MatDestroy(&C));
120 1 : PetscCall(MatDestroy(&K));
121 1 : PetscCall(SlepcFinalize());
122 : return 0;
123 : }
124 :
125 : /*TEST
126 :
127 : build:
128 : requires: complex
129 :
130 : test:
131 : suffix: 1
132 : requires: complex
133 :
134 : TEST*/
|
