Actual source code: butterfly.c
slepc-3.11.2 2019-07-30
1: /*
2: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
3: SLEPc - Scalable Library for Eigenvalue Problem Computations
4: Copyright (c) 2002-2019, 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: */
10: /*
11: This example implements one of the problems found at
12: NLEVP: A Collection of Nonlinear Eigenvalue Problems,
13: The University of Manchester.
14: The details of the collection can be found at:
15: [1] T. Betcke et al., "NLEVP: A Collection of Nonlinear Eigenvalue
16: Problems", ACM Trans. Math. Software 39(2), Article 7, 2013.
18: The butterfly problem is a quartic PEP with T-even structure.
19: */
21: static char help[] = "Quartic polynomial eigenproblem with T-even structure.\n\n"
22: "The command line options are:\n"
23: " -m <m>, grid size, the dimension of the matrices is n=m*m.\n"
24: " -c <array>, problem parameters, must be 10 comma-separated real values.\n\n";
26: #include <slepcpep.h>
28: #define NMAT 5
30: int main(int argc,char **argv)
31: {
32: Mat A[NMAT]; /* problem matrices */
33: PEP pep; /* polynomial eigenproblem solver context */
34: PetscInt n,m=8,k,II,Istart,Iend,i,j;
35: PetscReal c[10] = { 0.6, 1.3, 1.3, 0.1, 0.1, 1.2, 1.0, 1.0, 1.2, 1.0 };
36: PetscBool flg;
37: PetscBool terse;
40: SlepcInitialize(&argc,&argv,(char*)0,help);if (ierr) return ierr;
42: PetscOptionsGetInt(NULL,NULL,"-m",&m,NULL);
43: n = m*m;
44: k = 10;
45: PetscOptionsGetRealArray(NULL,NULL,"-c",c,&k,&flg);
46: if (flg && k!=10) SETERRQ1(PETSC_COMM_WORLD,1,"The number of parameters -c should be 10, you provided %D",k);
47: PetscPrintf(PETSC_COMM_WORLD,"\nButterfly problem, n=%D (m=%D)\n\n",n,m);
49: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
50: Compute the polynomial matrices
51: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
53: /* initialize matrices */
54: for (i=0;i<NMAT;i++) {
55: MatCreate(PETSC_COMM_WORLD,&A[i]);
56: MatSetSizes(A[i],PETSC_DECIDE,PETSC_DECIDE,n,n);
57: MatSetFromOptions(A[i]);
58: MatSetUp(A[i]);
59: }
60: MatGetOwnershipRange(A[0],&Istart,&Iend);
62: /* A0 */
63: for (II=Istart;II<Iend;II++) {
64: i = II/m; j = II-i*m;
65: MatSetValue(A[0],II,II,4.0*c[0]/6.0+4.0*c[1]/6.0,INSERT_VALUES);
66: if (j>0) { MatSetValue(A[0],II,II-1,c[0]/6.0,INSERT_VALUES); }
67: if (j<m-1) { MatSetValue(A[0],II,II+1,c[0]/6.0,INSERT_VALUES); }
68: if (i>0) { MatSetValue(A[0],II,II-m,c[1]/6.0,INSERT_VALUES); }
69: if (i<m-1) { MatSetValue(A[0],II,II+m,c[1]/6.0,INSERT_VALUES); }
70: }
72: /* A1 */
73: for (II=Istart;II<Iend;II++) {
74: i = II/m; j = II-i*m;
75: if (j>0) { MatSetValue(A[1],II,II-1,c[2],INSERT_VALUES); }
76: if (j<m-1) { MatSetValue(A[1],II,II+1,-c[2],INSERT_VALUES); }
77: if (i>0) { MatSetValue(A[1],II,II-m,c[3],INSERT_VALUES); }
78: if (i<m-1) { MatSetValue(A[1],II,II+m,-c[3],INSERT_VALUES); }
79: }
81: /* A2 */
82: for (II=Istart;II<Iend;II++) {
83: i = II/m; j = II-i*m;
84: MatSetValue(A[2],II,II,-2.0*c[4]-2.0*c[5],INSERT_VALUES);
85: if (j>0) { MatSetValue(A[2],II,II-1,c[4],INSERT_VALUES); }
86: if (j<m-1) { MatSetValue(A[2],II,II+1,c[4],INSERT_VALUES); }
87: if (i>0) { MatSetValue(A[2],II,II-m,c[5],INSERT_VALUES); }
88: if (i<m-1) { MatSetValue(A[2],II,II+m,c[5],INSERT_VALUES); }
89: }
91: /* A3 */
92: for (II=Istart;II<Iend;II++) {
93: i = II/m; j = II-i*m;
94: if (j>0) { MatSetValue(A[3],II,II-1,c[6],INSERT_VALUES); }
95: if (j<m-1) { MatSetValue(A[3],II,II+1,-c[6],INSERT_VALUES); }
96: if (i>0) { MatSetValue(A[3],II,II-m,c[7],INSERT_VALUES); }
97: if (i<m-1) { MatSetValue(A[3],II,II+m,-c[7],INSERT_VALUES); }
98: }
100: /* A4 */
101: for (II=Istart;II<Iend;II++) {
102: i = II/m; j = II-i*m;
103: MatSetValue(A[4],II,II,2.0*c[8]+2.0*c[9],INSERT_VALUES);
104: if (j>0) { MatSetValue(A[4],II,II-1,-c[8],INSERT_VALUES); }
105: if (j<m-1) { MatSetValue(A[4],II,II+1,-c[8],INSERT_VALUES); }
106: if (i>0) { MatSetValue(A[4],II,II-m,-c[9],INSERT_VALUES); }
107: if (i<m-1) { MatSetValue(A[4],II,II+m,-c[9],INSERT_VALUES); }
108: }
110: /* assemble matrices */
111: for (i=0;i<NMAT;i++) {
112: MatAssemblyBegin(A[i],MAT_FINAL_ASSEMBLY);
113: }
114: for (i=0;i<NMAT;i++) {
115: MatAssemblyEnd(A[i],MAT_FINAL_ASSEMBLY);
116: }
118: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
119: Create the eigensolver and solve the problem
120: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
122: PEPCreate(PETSC_COMM_WORLD,&pep);
123: PEPSetOperators(pep,NMAT,A);
124: PEPSetFromOptions(pep);
125: PEPSolve(pep);
127: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
128: Display solution and clean up
129: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
131: /* show detailed info unless -terse option is given by user */
132: PetscOptionsHasName(NULL,NULL,"-terse",&terse);
133: if (terse) {
134: PEPErrorView(pep,PEP_ERROR_BACKWARD,NULL);
135: } else {
136: PetscViewerPushFormat(PETSC_VIEWER_STDOUT_WORLD,PETSC_VIEWER_ASCII_INFO_DETAIL);
137: PEPReasonView(pep,PETSC_VIEWER_STDOUT_WORLD);
138: PEPErrorView(pep,PEP_ERROR_BACKWARD,PETSC_VIEWER_STDOUT_WORLD);
139: PetscViewerPopFormat(PETSC_VIEWER_STDOUT_WORLD);
140: }
141: PEPDestroy(&pep);
142: for (i=0;i<NMAT;i++) {
143: MatDestroy(&A[i]);
144: }
145: SlepcFinalize();
146: return ierr;
147: }
149: /*TEST
151: testset:
152: args: -pep_nev 4 -st_type sinvert -pep_target 0.01 -terse
153: requires: !complex
154: output_file: output/butterfly_1.out
155: test:
156: suffix: 1_toar
157: args: -pep_type toar -pep_toar_restart 0.3
158: test:
159: suffix: 1_linear
160: args: -pep_type linear
162: test:
163: suffix: 2
164: args: -pep_type {{toar linear}} -pep_nev 4 -terse
165: requires: double !complex !define(PETSC_USE_64BIT_INDICES)
167: TEST*/