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[] = "Generalized eigenproblem, illustrates setting MUMPS options.\n\n"
12 : "The problem is Ax = lambda Bx, with:\n"
13 : " A = Laplacian operator in 2-D\n"
14 : " B = diagonal matrix with all values equal to 4\n\n"
15 : "The command line options are:\n"
16 : " -n <n>, where <n> = number of grid subdivisions in x dimension.\n"
17 : " -m <m>, where <m> = number of grid subdivisions in y dimension.\n\n";
18 :
19 : #include <slepceps.h>
20 :
21 1 : int main(int argc,char **argv)
22 : {
23 1 : Mat A,B;
24 : #if defined(PETSC_HAVE_MUMPS)
25 : Mat K;
26 : #endif
27 1 : EPS eps;
28 1 : EPSType type;
29 1 : ST st;
30 1 : KSP ksp;
31 1 : PC pc;
32 1 : PetscInt N,n=10,m=12,Istart,Iend,II,nev,i,j;
33 1 : PetscBool flag,terse;
34 :
35 1 : PetscFunctionBeginUser;
36 1 : PetscCall(SlepcInitialize(&argc,&argv,(char*)0,help));
37 :
38 1 : PetscCall(PetscOptionsGetInt(NULL,NULL,"-n",&n,NULL));
39 1 : PetscCall(PetscOptionsGetInt(NULL,NULL,"-m",&m,&flag));
40 1 : N = n*m;
41 1 : PetscCall(PetscPrintf(PETSC_COMM_WORLD,"\nGeneralized Eigenproblem, N=%" PetscInt_FMT " (%" PetscInt_FMT "x%" PetscInt_FMT " grid)\n\n",N,n,m));
42 :
43 : /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
44 : Compute the matrices that define the eigensystem, Ax=kBx
45 : - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
46 :
47 1 : PetscCall(MatCreate(PETSC_COMM_WORLD,&A));
48 1 : PetscCall(MatSetSizes(A,PETSC_DECIDE,PETSC_DECIDE,N,N));
49 1 : PetscCall(MatSetFromOptions(A));
50 :
51 1 : PetscCall(MatCreate(PETSC_COMM_WORLD,&B));
52 1 : PetscCall(MatSetSizes(B,PETSC_DECIDE,PETSC_DECIDE,N,N));
53 1 : PetscCall(MatSetFromOptions(B));
54 :
55 1 : PetscCall(MatGetOwnershipRange(A,&Istart,&Iend));
56 121 : for (II=Istart;II<Iend;II++) {
57 120 : i = II/n; j = II-i*n;
58 120 : if (i>0) PetscCall(MatSetValue(A,II,II-n,-1.0,INSERT_VALUES));
59 120 : if (i<m-1) PetscCall(MatSetValue(A,II,II+n,-1.0,INSERT_VALUES));
60 120 : if (j>0) PetscCall(MatSetValue(A,II,II-1,-1.0,INSERT_VALUES));
61 120 : if (j<n-1) PetscCall(MatSetValue(A,II,II+1,-1.0,INSERT_VALUES));
62 120 : PetscCall(MatSetValue(A,II,II,4.0,INSERT_VALUES));
63 120 : PetscCall(MatSetValue(B,II,II,4.0,INSERT_VALUES));
64 : }
65 :
66 1 : PetscCall(MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY));
67 1 : PetscCall(MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY));
68 1 : PetscCall(MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY));
69 1 : PetscCall(MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY));
70 :
71 : /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
72 : Create the eigensolver and set various options
73 : - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
74 :
75 : /*
76 : Create eigensolver context
77 : */
78 1 : PetscCall(EPSCreate(PETSC_COMM_WORLD,&eps));
79 :
80 : /*
81 : Set operators. In this case, it is a generalized eigenvalue problem
82 : */
83 1 : PetscCall(EPSSetOperators(eps,A,B));
84 1 : PetscCall(EPSSetProblemType(eps,EPS_GNHEP));
85 :
86 : /*
87 : Set some solver options
88 : */
89 1 : PetscCall(EPSSetTarget(eps,1.3));
90 1 : PetscCall(EPSSetDimensions(eps,2,PETSC_DEFAULT,PETSC_DEFAULT));
91 1 : PetscCall(EPSGetST(eps,&st));
92 1 : PetscCall(STSetType(st,STSINVERT));
93 :
94 1 : PetscCall(STGetKSP(st,&ksp));
95 1 : PetscCall(KSPSetType(ksp,KSPPREONLY));
96 1 : PetscCall(KSPGetPC(ksp,&pc));
97 1 : PetscCall(PCSetType(pc,PCLU));
98 :
99 : /*
100 : Set MUMPS options if available
101 : */
102 : #if defined(PETSC_HAVE_MUMPS)
103 : PetscCall(PCFactorSetMatSolverType(pc,MATSOLVERMUMPS));
104 : /* the next line is required to force the creation of the ST operator and its passing to KSP */
105 : PetscCall(STGetOperator(st,NULL));
106 : PetscCall(PCFactorSetUpMatSolverType(pc));
107 : PetscCall(PCFactorGetMatrix(pc,&K));
108 : PetscCall(MatMumpsSetIcntl(K,14,50));
109 : PetscCall(MatMumpsSetCntl(K,3,1e-12));
110 : #endif
111 :
112 : /*
113 : Let the user change settings at runtime
114 : */
115 1 : PetscCall(EPSSetFromOptions(eps));
116 :
117 : /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
118 : Solve the eigensystem
119 : - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
120 :
121 1 : PetscCall(EPSSolve(eps));
122 :
123 : /*
124 : Optional: Get some information from the solver and display it
125 : */
126 1 : PetscCall(EPSGetType(eps,&type));
127 1 : PetscCall(PetscPrintf(PETSC_COMM_WORLD," Solution method: %s\n\n",type));
128 1 : PetscCall(EPSGetDimensions(eps,&nev,NULL,NULL));
129 1 : PetscCall(PetscPrintf(PETSC_COMM_WORLD," Number of requested eigenvalues: %" PetscInt_FMT "\n",nev));
130 :
131 : /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
132 : Display solution and clean up
133 : - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
134 :
135 : /* show detailed info unless -terse option is given by user */
136 1 : PetscCall(PetscOptionsHasName(NULL,NULL,"-terse",&terse));
137 1 : if (terse) PetscCall(EPSErrorView(eps,EPS_ERROR_RELATIVE,NULL));
138 : else {
139 0 : PetscCall(PetscViewerPushFormat(PETSC_VIEWER_STDOUT_WORLD,PETSC_VIEWER_ASCII_INFO_DETAIL));
140 0 : PetscCall(EPSConvergedReasonView(eps,PETSC_VIEWER_STDOUT_WORLD));
141 0 : PetscCall(EPSErrorView(eps,EPS_ERROR_RELATIVE,PETSC_VIEWER_STDOUT_WORLD));
142 0 : PetscCall(PetscViewerPopFormat(PETSC_VIEWER_STDOUT_WORLD));
143 : }
144 1 : PetscCall(EPSDestroy(&eps));
145 1 : PetscCall(MatDestroy(&A));
146 1 : PetscCall(MatDestroy(&B));
147 1 : PetscCall(SlepcFinalize());
148 : return 0;
149 : }
150 :
151 : /*TEST
152 :
153 : testset:
154 : args: -terse
155 : output_file: output/ex43_1.out
156 : test:
157 : suffix: 1
158 : test:
159 : suffix: 2
160 : nsize: 2
161 : args: -st_pc_factor_mat_solver_type mumps
162 : requires: mumps
163 :
164 : TEST*/
|
