OOFEM  2.4
OOFEM.org - Object Oriented Finite Element Solver
slepcsolver.C
Go to the documentation of this file.
1 /*
2  *
3  * ##### ##### ###### ###### ### ###
4  * ## ## ## ## ## ## ## ### ##
5  * ## ## ## ## #### #### ## # ##
6  * ## ## ## ## ## ## ## ##
7  * ## ## ## ## ## ## ## ##
8  * ##### ##### ## ###### ## ##
9  *
10  *
11  * OOFEM : Object Oriented Finite Element Code
12  *
13  * Copyright (C) 1993 - 2013 Borek Patzak
14  *
15  *
16  *
17  * Czech Technical University, Faculty of Civil Engineering,
18  * Department of Structural Mechanics, 166 29 Prague, Czech Republic
19  *
20  * This library is free software; you can redistribute it and/or
21  * modify it under the terms of the GNU Lesser General Public
22  * License as published by the Free Software Foundation; either
23  * version 2.1 of the License, or (at your option) any later version.
24  *
25  * This program is distributed in the hope that it will be useful,
26  * but WITHOUT ANY WARRANTY; without even the implied warranty of
27  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
28  * Lesser General Public License for more details.
29  *
30  * You should have received a copy of the GNU Lesser General Public
31  * License along with this library; if not, write to the Free Software
32  * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
33  */
34 
35 #include "slepcsolver.h"
36 #include "petscsparsemtrx.h"
37 #include "engngm.h"
38 #include "floatarray.h"
39 #include "verbose.h"
40 #include "domain.h"
41 #include "classfactory.h"
42 
43 #define TIME_REPORT
44 
45 #ifdef TIME_REPORT
46  #include "timer.h"
47 #endif
48 
49 namespace oofem {
51 
52 
54 {
55  A = B = NULL;
56  epsInit = false;
57 }
58 
59 
61 {
62  if ( epsInit ) {
63  EPSDestroy(& eps);
64  }
65 }
66 
68 SLEPcSolver :: solve(SparseMtrx &a, SparseMtrx &b, FloatArray &_eigv, FloatMatrix &_r, double rtol, int nroot)
69 {
70  PetscErrorCode ierr;
71  int size;
72  ST st;
73 
74  // first check whether Lhs is defined
75 
76  if ( a.giveNumberOfRows() != a.giveNumberOfColumns() ||
79  OOFEM_ERROR("matrices size mismatch");
80  }
81 
82  A = dynamic_cast< PetscSparseMtrx * >(&a);
83  B = dynamic_cast< PetscSparseMtrx * >(&b);
84 
85  if ( !A || !B ) {
86  OOFEM_ERROR("PetscSparseMtrx Expected");
87  }
88 
89  size = engngModel->giveParallelContext( A->giveDomainIndex() )->giveNumberOfNaturalEqs(); // A->giveLeqs();
90 
91  _r.resize(size, nroot);
92  _eigv.resize(nroot);
93 
94 
95  /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
96  * Create the eigensolver and set various options
97  * - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
98  int nconv, nite;
99  EPSConvergedReason reason;
100 
101 #ifdef TIME_REPORT
102  Timer timer;
103  timer.startTimer();
104 #endif
105 
106  if ( !epsInit ) {
107  /*
108  * Create eigensolver context
109  */
110 #ifdef __PARALLEL_MODE
111  MPI_Comm comm = engngModel->giveParallelComm();
112 #else
113  MPI_Comm comm = PETSC_COMM_SELF;
114 #endif
115  ierr = EPSCreate(comm, & eps);
116  CHKERRQ(ierr);
117  epsInit = true;
118  }
119 
120  /*
121  * Set operators. In this case, it is a generalized eigenvalue problem
122  */
123 
124  ierr = EPSSetOperators( eps, * A->giveMtrx(), * B->giveMtrx() );
125  CHKERRQ(ierr);
126  ierr = EPSSetProblemType(eps, EPS_GHEP);
127  CHKERRQ(ierr);
128  ierr = EPSGetST(eps, & st);
129  CHKERRQ(ierr);
130  ierr = STSetType(st, STSINVERT);
131  CHKERRQ(ierr);
132  ierr = STSetMatStructure(st, SAME_NONZERO_PATTERN);
133  CHKERRQ(ierr);
134  ierr = EPSSetTolerances(eps, ( PetscReal ) rtol, PETSC_DECIDE);
135  CHKERRQ(ierr);
136  ierr = EPSSetDimensions(eps, ( PetscInt ) nroot, PETSC_DECIDE, PETSC_DECIDE);
137  CHKERRQ(ierr);
138  ierr = EPSSetWhichEigenpairs(eps, EPS_SMALLEST_MAGNITUDE);
139  CHKERRQ(ierr);
140 
141  /*
142  * Set solver parameters at runtime
143  */
144 
145  ierr = EPSSetFromOptions(eps);
146  CHKERRQ(ierr);
147 
148  /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
149  * Solve the eigensystem
150  * - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
151 
152  ierr = EPSSolve(eps);
153  CHKERRQ(ierr);
154 
155  ierr = EPSGetConvergedReason(eps, & reason);
156  CHKERRQ(ierr);
157  ierr = EPSGetIterationNumber(eps, & nite);
158  CHKERRQ(ierr);
159  OOFEM_LOG_INFO("SLEPcSolver::solve EPSConvergedReason: %d, number of iterations: %d\n", reason, nite);
160 
161  ierr = EPSGetConverged(eps, & nconv);
162  CHKERRQ(ierr);
163 
164  if ( nconv > 0 ) {
165  OOFEM_LOG_INFO("SLEPcSolver :: solveYourselfAt: Convergence reached for RTOL=%20.15f", rtol);
166  PetscScalar kr;
167  Vec Vr;
168 
169  ierr = MatGetVecs(* B->giveMtrx(), PETSC_NULL, & Vr);
170  CHKERRQ(ierr);
171 
172  FloatArray Vr_loc;
173 
174  for ( int i = 0; i < nconv && i < nroot; i++ ) {
175  ierr = EPSGetEigenpair(eps, nconv - i - 1, & kr, PETSC_NULL, Vr, PETSC_NULL);
176  CHKERRQ(ierr);
177 
178  //Store the eigenvalue
179  _eigv.at(i + 1) = kr;
180 
181  //Store the eigenvector
182  A->scatterG2L(Vr, Vr_loc);
183  for ( int j = 0; j < size; j++ ) {
184  _r.at(j + 1, i + 1) = Vr_loc.at(j + 1);
185  }
186  }
187 
188  ierr = VecDestroy(& Vr);
189  CHKERRQ(ierr);
190  } else {
191  OOFEM_ERROR("No converged eigenpairs");
192  }
193 
194 #ifdef TIME_REPORT
195  timer.stopTimer();
196  OOFEM_LOG_INFO( "SLEPcSolver info: user time consumed by solution: %.2fs\n", timer.getUtime() );
197 #endif
198 
199  return NM_Success;
200 }
201 } // end namespace oofem
bool epsInit
Flag if context initialized.
Definition: slepcsolver.h:60
#define NM_Success
Numerical method exited with success.
Definition: nmstatus.h:47
EPS eps
Eigenvalue solver context.
Definition: slepcsolver.h:58
Class and object Domain.
Definition: domain.h:115
Base class for all matrices stored in sparse format.
Definition: sparsemtrx.h:60
SLEPcSolver(Domain *d, EngngModel *m)
Definition: slepcsolver.C:53
double & at(int i)
Coefficient access function.
Definition: floatarray.h:131
unsigned long NM_Status
Mask defining NumMetod Status; which can be asked after finishing computation by Numerical Method...
Definition: nmstatus.h:44
virtual ~SLEPcSolver()
Definition: slepcsolver.C:60
REGISTER_GeneralizedEigenValueSolver(InverseIteration, GES_InverseIt)
PetscSparseMtrx * B
Definition: slepcsolver.h:56
int giveNumberOfRows() const
Returns number of rows of receiver.
Definition: sparsemtrx.h:114
#define OOFEM_LOG_INFO(...)
Definition: logger.h:127
#define OOFEM_ERROR(...)
Definition: error.h:61
virtual ParallelContext * giveParallelContext(int n)
Returns the parallel context corresponding to given domain (n) and unknown type Default implementatio...
Definition: engngm.C:1745
PetscSparseMtrx * A
Definition: slepcsolver.h:55
double at(int i, int j) const
Coefficient access function.
Definition: floatmatrix.h:176
MPI_Comm giveParallelComm()
Returns the communication object of reciever.
Definition: engngm.h:549
virtual NM_Status solve(SparseMtrx &a, SparseMtrx &b, FloatArray &v, FloatMatrix &x, double rtol, int nroot)
Solves the given sparse generalized eigen value system of equations .
Definition: slepcsolver.C:68
Initializes the variable VERBOSE, in order to get a few intermediate messages on screen: beginning an...
Class representing vector of real numbers.
Definition: floatarray.h:82
Implementation of matrix containing floating point numbers.
Definition: floatmatrix.h:94
void resize(int rows, int cols)
Checks size of receiver towards requested bounds.
Definition: floatmatrix.C:1358
Class implementing single timer, providing wall clock and user time capabilities. ...
Definition: timer.h:46
int giveNumberOfColumns() const
Returns number of columns of receiver.
Definition: sparsemtrx.h:116
Abstract base class representing the "problem" under consideration.
Definition: engngm.h:181
the oofem namespace is to define a context or scope in which all oofem names are defined.
int scatterG2L(Vec src, FloatArray &dest) const
Scatters global vector to natural one.
This base class is an abstraction for all numerical methods solving sparse linear system of equations...
EngngModel * engngModel
Pointer to engineering model.
Definition: nummet.h:86
void startTimer()
Definition: timer.C:69
This class provides an sparse matrix interface to PETSc Matrices.
double getUtime()
Returns total user time elapsed in seconds.
Definition: timer.C:105
void stopTimer()
Definition: timer.C:77
void resize(int s)
Resizes receiver towards requested size.
Definition: floatarray.C:631

This page is part of the OOFEM documentation. Copyright (c) 2011 Borek Patzak
Project e-mail: info@oofem.org
Generated at Tue Jan 2 2018 20:07:31 for OOFEM by doxygen 1.8.11 written by Dimitri van Heesch, © 1997-2011