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[] = "Illustrates use of PEPSetEigenvalueComparison().\n\n"
12 : "Based on butterfly.c. The command line options are:\n"
13 : " -m <m>, grid size, the dimension of the matrices is n=m*m.\n"
14 : " -c <array>, problem parameters, must be 10 comma-separated real values.\n\n";
15 :
16 : #include <slepcpep.h>
17 :
18 : #define NMAT 5
19 :
20 : /*
21 : Function for user-defined eigenvalue ordering criterion.
22 :
23 : Given two eigenvalues ar+i*ai and br+i*bi, the subroutine must choose
24 : one of them as the preferred one according to the criterion.
25 : In this example, eigenvalues are sorted by magnitude but those with
26 : positive real part are preferred.
27 : */
28 335 : PetscErrorCode MyEigenSort(PetscScalar ar,PetscScalar ai,PetscScalar br,PetscScalar bi,PetscInt *r,void *ctx)
29 : {
30 335 : PetscReal rea,reb;
31 :
32 335 : PetscFunctionBeginUser;
33 : #if defined(PETSC_USE_COMPLEX)
34 : rea = PetscRealPart(ar); reb = PetscRealPart(br);
35 : #else
36 335 : rea = ar; reb = br;
37 : #endif
38 335 : *r = rea<0.0? 1: (reb<0.0? -1: PetscSign(SlepcAbsEigenvalue(br,bi)-SlepcAbsEigenvalue(ar,ai)));
39 335 : PetscFunctionReturn(PETSC_SUCCESS);
40 : }
41 :
42 1 : int main(int argc,char **argv)
43 : {
44 1 : Mat A[NMAT]; /* problem matrices */
45 1 : PEP pep; /* polynomial eigenproblem solver context */
46 1 : PetscInt n,m=8,k,II,Istart,Iend,i,j;
47 1 : PetscReal c[10] = { 0.6, 1.3, 1.3, 0.1, 0.1, 1.2, 1.0, 1.0, 1.2, 1.0 };
48 1 : PetscBool flg;
49 1 : PetscBool terse;
50 1 : const char *prefix;
51 :
52 1 : PetscFunctionBeginUser;
53 1 : PetscCall(SlepcInitialize(&argc,&argv,(char*)0,help));
54 :
55 1 : PetscCall(PetscOptionsGetInt(NULL,NULL,"-m",&m,NULL));
56 1 : n = m*m;
57 1 : k = 10;
58 1 : PetscCall(PetscOptionsGetRealArray(NULL,NULL,"-c",c,&k,&flg));
59 1 : PetscCheck(!flg || k==10,PETSC_COMM_WORLD,PETSC_ERR_USER,"The number of parameters -c should be 10, you provided %" PetscInt_FMT,k);
60 1 : PetscCall(PetscPrintf(PETSC_COMM_WORLD,"\nButterfly problem, n=%" PetscInt_FMT " (m=%" PetscInt_FMT ")\n\n",n,m));
61 :
62 : /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
63 : Compute the polynomial matrices
64 : - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
65 :
66 : /* initialize matrices */
67 6 : for (i=0;i<NMAT;i++) {
68 5 : PetscCall(MatCreate(PETSC_COMM_WORLD,&A[i]));
69 5 : PetscCall(MatSetSizes(A[i],PETSC_DECIDE,PETSC_DECIDE,n,n));
70 5 : PetscCall(MatSetFromOptions(A[i]));
71 : }
72 1 : PetscCall(MatGetOwnershipRange(A[0],&Istart,&Iend));
73 :
74 : /* A0 */
75 65 : for (II=Istart;II<Iend;II++) {
76 64 : i = II/m; j = II-i*m;
77 64 : PetscCall(MatSetValue(A[0],II,II,4.0*c[0]/6.0+4.0*c[1]/6.0,INSERT_VALUES));
78 64 : if (j>0) PetscCall(MatSetValue(A[0],II,II-1,c[0]/6.0,INSERT_VALUES));
79 64 : if (j<m-1) PetscCall(MatSetValue(A[0],II,II+1,c[0]/6.0,INSERT_VALUES));
80 64 : if (i>0) PetscCall(MatSetValue(A[0],II,II-m,c[1]/6.0,INSERT_VALUES));
81 64 : if (i<m-1) PetscCall(MatSetValue(A[0],II,II+m,c[1]/6.0,INSERT_VALUES));
82 : }
83 :
84 : /* A1 */
85 65 : for (II=Istart;II<Iend;II++) {
86 64 : i = II/m; j = II-i*m;
87 64 : if (j>0) PetscCall(MatSetValue(A[1],II,II-1,c[2],INSERT_VALUES));
88 64 : if (j<m-1) PetscCall(MatSetValue(A[1],II,II+1,-c[2],INSERT_VALUES));
89 64 : if (i>0) PetscCall(MatSetValue(A[1],II,II-m,c[3],INSERT_VALUES));
90 64 : if (i<m-1) PetscCall(MatSetValue(A[1],II,II+m,-c[3],INSERT_VALUES));
91 : }
92 :
93 : /* A2 */
94 65 : for (II=Istart;II<Iend;II++) {
95 64 : i = II/m; j = II-i*m;
96 64 : PetscCall(MatSetValue(A[2],II,II,-2.0*c[4]-2.0*c[5],INSERT_VALUES));
97 64 : if (j>0) PetscCall(MatSetValue(A[2],II,II-1,c[4],INSERT_VALUES));
98 64 : if (j<m-1) PetscCall(MatSetValue(A[2],II,II+1,c[4],INSERT_VALUES));
99 64 : if (i>0) PetscCall(MatSetValue(A[2],II,II-m,c[5],INSERT_VALUES));
100 64 : if (i<m-1) PetscCall(MatSetValue(A[2],II,II+m,c[5],INSERT_VALUES));
101 : }
102 :
103 : /* A3 */
104 65 : for (II=Istart;II<Iend;II++) {
105 64 : i = II/m; j = II-i*m;
106 64 : if (j>0) PetscCall(MatSetValue(A[3],II,II-1,c[6],INSERT_VALUES));
107 64 : if (j<m-1) PetscCall(MatSetValue(A[3],II,II+1,-c[6],INSERT_VALUES));
108 64 : if (i>0) PetscCall(MatSetValue(A[3],II,II-m,c[7],INSERT_VALUES));
109 64 : if (i<m-1) PetscCall(MatSetValue(A[3],II,II+m,-c[7],INSERT_VALUES));
110 : }
111 :
112 : /* A4 */
113 65 : for (II=Istart;II<Iend;II++) {
114 64 : i = II/m; j = II-i*m;
115 64 : PetscCall(MatSetValue(A[4],II,II,2.0*c[8]+2.0*c[9],INSERT_VALUES));
116 64 : if (j>0) PetscCall(MatSetValue(A[4],II,II-1,-c[8],INSERT_VALUES));
117 64 : if (j<m-1) PetscCall(MatSetValue(A[4],II,II+1,-c[8],INSERT_VALUES));
118 64 : if (i>0) PetscCall(MatSetValue(A[4],II,II-m,-c[9],INSERT_VALUES));
119 64 : if (i<m-1) PetscCall(MatSetValue(A[4],II,II+m,-c[9],INSERT_VALUES));
120 : }
121 :
122 : /* assemble matrices */
123 6 : for (i=0;i<NMAT;i++) PetscCall(MatAssemblyBegin(A[i],MAT_FINAL_ASSEMBLY));
124 6 : for (i=0;i<NMAT;i++) PetscCall(MatAssemblyEnd(A[i],MAT_FINAL_ASSEMBLY));
125 :
126 : /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
127 : Create the eigensolver and solve the problem
128 : - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
129 :
130 1 : PetscCall(PEPCreate(PETSC_COMM_WORLD,&pep));
131 1 : PetscCall(PEPSetOptionsPrefix(pep,"check_"));
132 1 : PetscCall(PEPAppendOptionsPrefix(pep,"myprefix_"));
133 1 : PetscCall(PEPGetOptionsPrefix(pep,&prefix));
134 1 : PetscCall(PetscPrintf(PETSC_COMM_WORLD,"PEP prefix is currently: %s\n\n",prefix));
135 :
136 1 : PetscCall(PEPSetOperators(pep,NMAT,A));
137 1 : PetscCall(PEPSetEigenvalueComparison(pep,MyEigenSort,NULL));
138 1 : PetscCall(PEPSetFromOptions(pep));
139 1 : PetscCall(PEPSolve(pep));
140 :
141 : /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
142 : Display solution and clean up
143 : - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
144 :
145 : /* show detailed info unless -terse option is given by user */
146 1 : PetscCall(PetscOptionsHasName(NULL,NULL,"-terse",&terse));
147 1 : if (terse) PetscCall(PEPErrorView(pep,PEP_ERROR_BACKWARD,NULL));
148 : else {
149 0 : PetscCall(PetscViewerPushFormat(PETSC_VIEWER_STDOUT_WORLD,PETSC_VIEWER_ASCII_INFO_DETAIL));
150 0 : PetscCall(PEPConvergedReasonView(pep,PETSC_VIEWER_STDOUT_WORLD));
151 0 : PetscCall(PEPErrorView(pep,PEP_ERROR_BACKWARD,PETSC_VIEWER_STDOUT_WORLD));
152 0 : PetscCall(PetscViewerPopFormat(PETSC_VIEWER_STDOUT_WORLD));
153 : }
154 1 : PetscCall(PEPDestroy(&pep));
155 6 : for (i=0;i<NMAT;i++) PetscCall(MatDestroy(&A[i]));
156 1 : PetscCall(SlepcFinalize());
157 : return 0;
158 : }
159 :
160 : /*TEST
161 :
162 : test:
163 : args: -check_myprefix_pep_nev 4 -terse
164 : requires: double
165 :
166 : TEST*/
|
