1 // @HEADER
2 // ***********************************************************************
3 //
4 // Teuchos: Common Tools Package
5 // Copyright (2004) Sandia Corporation
6 //
7 // Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
8 // license for use of this work by or on behalf of the U.S. Government.
9 //
10 // Redistribution and use in source and binary forms, with or without
11 // modification, are permitted provided that the following conditions are
12 // met:
13 //
14 // 1. Redistributions of source code must retain the above copyright
15 // notice, this list of conditions and the following disclaimer.
16 //
17 // 2. Redistributions in binary form must reproduce the above copyright
18 // notice, this list of conditions and the following disclaimer in the
19 // documentation and/or other materials provided with the distribution.
20 //
21 // 3. Neither the name of the Corporation nor the names of the
22 // contributors may be used to endorse or promote products derived from
23 // this software without specific prior written permission.
24 //
25 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
26 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
27 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
28 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
29 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
30 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
31 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
32 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
33 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
34 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
35 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
36 //
37 // Questions? Contact Michael A. Heroux (maherou@sandia.gov)
38 //
39 // ***********************************************************************
40 // @HEADER
41
42 #ifndef _TEUCHOS_SERIALBANDDENSESOLVER_HPP_
43 #define _TEUCHOS_SERIALBANDDENSESOLVER_HPP_
44 /// \file Teuchos_SerialBandDenseSolver.hpp
45 ///
46 /// Declaration and definition of SerialBandDenseSolver,
47 /// a templated class for solving banded dense linear systems.
48
49 #include "Teuchos_CompObject.hpp"
50 #include "Teuchos_BLAS.hpp"
51 #include "Teuchos_LAPACK.hpp"
52 #include "Teuchos_RCP.hpp"
53 #include "Teuchos_ConfigDefs.hpp"
54 #include "Teuchos_SerialDenseMatrix.hpp"
55 #include "Teuchos_SerialBandDenseMatrix.hpp"
56 #include "Teuchos_SerialDenseSolver.hpp"
57 #include "Teuchos_ScalarTraits.hpp"
58
59 #ifdef HAVE_TEUCHOSNUMERICS_EIGEN
60 #include "Eigen/Dense"
61 #endif
62
63 namespace Teuchos {
64
65 /*! \class SerialBandDenseSolver
66 \brief A class for representing and solving banded dense linear systems.
67 \tparam OrdinalType The index type used by the linear algebra implementation.
68 This should always be \c int.
69 \tparam ScalarType The type of entries in the matrix.
70
71 \section Teuchos_SerialBandDenseSolver_Intro Introduction
72
73 This class defines a banded dense matrix, which may have any number
74 of rows or columns (not necessarily equal). It's called "serial"
75 because the matrix lives in a single memory space. Thus, it's the
76 kind of matrix that one might give to the BLAS or LAPACK, not a
77 distributed matrix like one would give to ScaLAPACK.
78
79 This class also has methods for computing the (banded) LU
80 factorization of the matrix, and solving linear systems with the
81 matrix. We use instances of SerialDenseVector to represent the
82 right-hand side b or the solution x in the linear system \f$Ax=b\f$.
83 The instance of this class used store the banded matrix must contain
84 KL extra superdiagonals to store the L and U factors (see details
85 below).
86
87 Users have the option to do equilibration before factoring the matrix.
88 This may improve accuracy when solving ill-conditioned problems.
89
90 \section Teuchos_SerialBandDenseSolver_LAPACK SerialBandDenseSolver and LAPACK
91
92 Teuchos' LAPACK class wraps LAPACK's LU factorizations, including
93 the banded factorizations. It thus gives access to most of the same
94 functionality as SerialBandDenseSolver. The main difference is that
95 SerialBandDenseSolver offers a higher level of abstraction. It
96 hides the details of which LAPACK routines to call. Furthermore, if
97 you have built Teuchos with support for the third-party library
98 Eigen, SerialBandDenseSolver lets you solve linear systems for
99 <tt>ScalarType</tt> other than the four supported by the LAPACK
100 library.
101
102 \section Teuchos_SerialBandDenseSolver_Construct Constructing SerialBandDenseSolver objects
103
104 There is a single Teuchos::SerialBandDenseSolver constructor.
105 However, the matrix, right hand side and solution vectors must be
106 set prior to executing most methods in this class.
107
108 \subsection Teuchos_SerialBandDenseSolver_Construct_Setting Setting vectors used for linear solves
109
110 The matrix A, the left hand side X and the right hand side B (when
111 solving AX = B, for X), can be set by appropriate set methods. Each
112 of these three objects must be a SerialDenseMatrix or a
113 SerialDenseVector object. The set methods are as follows:
114 - setMatrix(): Sets the matrix
115 - setVectors() - Sets the left and right hand side vector(s)
116
117 \subsection Teuchos_SerialBandDenseSolver_Construct_Format Format of the matrix A
118
119 The SerialBandDenseMatrix must contain KL extra superdiagonals to store the L and U factors, where KL
120 is the upper bandwidth. Consider using the non-member conversion routines generalToBanded and bandedToGeneral if the
121 full SerialDenseMatrix is already in storage. However, it is more efficient simply to construct the
122 SerialBandDenseMatrix with the desired parameters and use the provided matrix access operators so
123 that the full rectangular matrix need not be stored. The conversion routine generalToBanded has a flag to store
124 the input Teuchos::SerialDenseMatrix in banded format with KL extra superdiagonals so this class can use it. Again,
125 it is more efficient to simply construct a Teuchos::SerialBandDenseMatrix object with KL extra superdiagonals than are
126 needed for the matrix data and fill the matrix using the matrix access operators.
127
128 See the documentation of Teuchos::SerialBandDenseMatrix for further details on the storage format.
129
130 \section Teuchos_SerialBandDenseSolver_Util Vector and Utility Functions
131
132 Once a Teuchos::SerialBandDenseSolver is constructed, several mathematical functions can be applied to
133 the object. Specifically:
134 <ul>
135 <li> Conversion between storage formats
136 <li> Factorizations
137 <li> Solves
138 <li> Condition estimates
139 <li> Equilibration
140 <li> Norms
141 </ul>
142
143 \section Teuchos_SerialBandDenseSolver_Strategies Strategies for Solving Linear Systems
144
145 In many cases, linear systems can be accurately solved by simply computing the LU factorization
146 of the matrix and then performing a forward back solve with a given set of right hand side vectors. However,
147 in some instances, the factorization may be very poorly conditioned and this simple approach may not work. In
148 these situations, equilibration and iterative refinement may improve the accuracy, or prevent a breakdown in
149 the factorization.
150
151 SerialBandDenseSolver will use equilibration with the factorization if, once the object
152 is constructed and \e before it is factored, you call the function factorWithEquilibration(true) to force
153 equilibration to be used. If you are uncertain if equilibration should be used, you may call the function
154 shouldEquilibrate() which will return true if equilibration could possibly help. shouldEquilibrate() uses
155 guidelines specified in the LAPACK User Guide, namely if SCOND < 0.1 and AMAX < Underflow or AMAX > Overflow, to
156 determine if equilibration \e might be useful.
157
158 SerialBandDenseSolver will use iterative refinement after a forward/back solve if you call
159 solveToRefinedSolution(true). It will also compute forward and backward error estimates if you call
160 estimateSolutionErrors(true). Access to the forward (back) error estimates is available via FERR() (BERR()).
161
162 Examples using SerialBandDenseSolver can be found in the Teuchos test directories.
163 */
164
165 template<typename OrdinalType, typename ScalarType>
166 class SerialBandDenseSolver : public CompObject, public Object, public BLAS<OrdinalType, ScalarType>,
167 public LAPACK<OrdinalType, ScalarType>
168 {
169
170 public:
171
172 typedef typename Teuchos::ScalarTraits<ScalarType>::magnitudeType MagnitudeType;
173 #ifdef HAVE_TEUCHOSNUMERICS_EIGEN
174 typedef typename Eigen::Matrix<ScalarType,Eigen::Dynamic,Eigen::Dynamic,Eigen::ColMajor> EigenMatrix;
175 typedef typename Eigen::internal::BandMatrix<ScalarType,Eigen::Dynamic,Eigen::Dynamic,Eigen::ColMajor> EigenBandMatrix;
176 typedef typename Eigen::Matrix<ScalarType,Eigen::Dynamic,1> EigenVector;
177 typedef typename Eigen::InnerStride<Eigen::Dynamic> EigenInnerStride;
178 typedef typename Eigen::OuterStride<Eigen::Dynamic> EigenOuterStride;
179 typedef typename Eigen::Map<EigenVector,0,EigenInnerStride > EigenVectorMap;
180 typedef typename Eigen::Map<const EigenVector,0,EigenInnerStride > EigenConstVectorMap;
181 typedef typename Eigen::Map<EigenMatrix,0,EigenOuterStride > EigenMatrixMap;
182 typedef typename Eigen::Map<const EigenMatrix,0,EigenOuterStride > EigenConstMatrixMap;
183 typedef typename Eigen::PermutationMatrix<Eigen::Dynamic,Eigen::Dynamic,OrdinalType> EigenPermutationMatrix;
184 typedef typename Eigen::Array<OrdinalType,Eigen::Dynamic,1> EigenOrdinalArray;
185 typedef typename Eigen::Map<EigenOrdinalArray> EigenOrdinalArrayMap;
186 typedef typename Eigen::Array<ScalarType,Eigen::Dynamic,1> EigenScalarArray;
187 typedef typename Eigen::Map<EigenScalarArray> EigenScalarArrayMap;
188 #endif
189
190 //! @name Constructor/Destructor Methods
191 //@{
192
193 //! Default constructor; matrix should be set using setMatrix(), LHS and RHS set with setVectors().
194 SerialBandDenseSolver();
195
196 //! SerialBandDenseSolver destructor.
197 virtual ~SerialBandDenseSolver();
198
199 //@}
200 //! @name Set Methods
201 //@{
202
203 //! Sets the pointer for coefficient matrix
204 int setMatrix(const RCP<SerialBandDenseMatrix<OrdinalType, ScalarType> >& AB);
205
206 //! Sets the pointers for left and right hand side vector(s).
207 /*! Row dimension of X must match column dimension of matrix A, row dimension of B
208 must match row dimension of A. X and B must have the same dimensions.
209 */
210 int setVectors(const RCP<SerialDenseMatrix<OrdinalType, ScalarType> >& X,
211 const RCP<SerialDenseMatrix<OrdinalType, ScalarType> >& B);
212
213 //@}
214 //! @name Strategy Modifying Methods
215 //@{
216
217 /// Set whether or not to equilibrate just before the matrix factorization.
218 /// This function must be called before the factorization is performed.
factorWithEquilibration(bool flag)219 void factorWithEquilibration(bool flag) {equilibrate_ = flag; return;}
220
221 /// If \c flag is true, causes all subsequent function calls to work with the transpose of \e this matrix, otherwise not.
222 ///
223 /// \note This interface will not work correctly for complex-valued linear systems, use solveWithTransposeFlag().
solveWithTranspose(bool flag)224 void solveWithTranspose(bool flag) {transpose_ = flag; if (flag) TRANS_ = Teuchos::TRANS; else TRANS_ = Teuchos::NO_TRANS; return;}
225
226 /// All subsequent function calls will work with the transpose-type set by this method (\c Teuchos::NO_TRANS, \c Teuchos::TRANS, and \c Teuchos::CONJ_TRANS).
227 /// \note This interface will allow correct behavior for complex-valued linear systems, solveWithTranspose() will not.
solveWithTransposeFlag(Teuchos::ETransp trans)228 void solveWithTransposeFlag(Teuchos::ETransp trans) {TRANS_ = trans; transpose_ = (trans != Teuchos::NO_TRANS) ? true : false;}
229
230 //! Set whether or not to use iterative refinement to improve solutions to linear systems.
solveToRefinedSolution(bool flag)231 void solveToRefinedSolution(bool flag) {refineSolution_ = flag; return;}
232
233 //! Causes all solves to estimate the forward and backward solution error.
234 /*! Error estimates will be in the arrays FERR and BERR, resp, after the solve step is complete.
235 These arrays are accessible via the FERR() and BERR() access functions.
236 */
237 void estimateSolutionErrors(bool flag);
238 //@}
239
240 //! @name Factor/Solve Methods
241 //@{
242
243 //! Computes the in-place LU factorization of the matrix using the LAPACK routine \e _GBTRF.
244 /*!
245 \return Integer error code, set to 0 if successful.
246 */
247 int factor();
248
249 //! Computes the solution X to AX = B for the \e this matrix and the B provided to SetVectors()..
250 /*!
251 \return Integer error code, set to 0 if successful.
252 */
253 int solve();
254
255 //! Computes the scaling vector S(i) = 1/sqrt(A(i,i)) of the \e this matrix.
256 /*!
257 \return Integer error code, set to 0 if successful. Otherwise returns the LAPACK error code INFO.
258 */
259 int computeEquilibrateScaling();
260
261 //! Equilibrates the \e this matrix.
262 /*!
263 \return Integer error code, set to 0 if successful. Otherwise returns the LAPACK error code INFO.
264 */
265 int equilibrateMatrix();
266
267 //! Equilibrates the current RHS.
268 /*!
269 \return Integer error code, set to 0 if successful. Otherwise returns the LAPACK error code INFO.
270 */
271 int equilibrateRHS();
272
273 //! Apply Iterative Refinement.
274 /*!
275 \return Integer error code, set to 0 if successful. Otherwise returns the LAPACK error code INFO.
276 */
277 int applyRefinement();
278
279 //! Unscales the solution vectors if equilibration was used to solve the system.
280 /*!
281 \return Integer error code, set to 0 if successful. Otherwise returns the LAPACK error code INFO.
282 */
283 int unequilibrateLHS();
284
285 //! Returns the reciprocal of the 1-norm condition number of the \e this matrix.
286 /*!
287 \param Value Out
288 On return contains the reciprocal of the 1-norm condition number of the \e this matrix.
289
290 \return Integer error code, set to 0 if successful. Otherwise returns the LAPACK error code INFO.
291 */
292 int reciprocalConditionEstimate(MagnitudeType & Value);
293 //@}
294
295 //! @name Query methods
296 //@{
297
298 //! Returns true if transpose of \e this matrix has and will be used.
transpose()299 bool transpose() {return(transpose_);}
300
301 //! Returns true if matrix is factored (factor available via AF() and LDAF()).
factored()302 bool factored() {return(factored_);}
303
304 //! Returns true if factor is equilibrated (factor available via AF() and LDAF()).
equilibratedA()305 bool equilibratedA() {return(equilibratedA_);}
306
307 //! Returns true if RHS is equilibrated (RHS available via B() and LDB()).
equilibratedB()308 bool equilibratedB() {return(equilibratedB_);}
309
310 //! Returns true if the LAPACK general rules for equilibration suggest you should equilibrate the system.
shouldEquilibrate()311 bool shouldEquilibrate() {computeEquilibrateScaling(); return(shouldEquilibrate_);}
312
313 //! Returns true if forward and backward error estimated have been computed (available via FERR() and BERR()).
solutionErrorsEstimated()314 bool solutionErrorsEstimated() {return(solutionErrorsEstimated_);}
315
316 //! Returns true if the condition number of the \e this matrix has been computed (value available via ReciprocalConditionEstimate()).
reciprocalConditionEstimated()317 bool reciprocalConditionEstimated() {return(reciprocalConditionEstimated_);}
318
319 //! Returns true if the current set of vectors has been solved.
solved()320 bool solved() {return(solved_);}
321
322 //! Returns true if the current set of vectors has been refined.
solutionRefined()323 bool solutionRefined() {return(solutionRefined_);}
324 //@}
325
326 //! @name Data Accessor methods
327 //@{
328
329 //! Returns pointer to current matrix.
getMatrix() const330 RCP<SerialBandDenseMatrix<OrdinalType, ScalarType> > getMatrix() const {return(Matrix_);}
331
332 //! Returns pointer to factored matrix (assuming factorization has been performed).
getFactoredMatrix() const333 RCP<SerialBandDenseMatrix<OrdinalType, ScalarType> > getFactoredMatrix() const {return(Factor_);}
334
335 //! Returns pointer to current LHS.
getLHS() const336 RCP<SerialDenseMatrix<OrdinalType, ScalarType> > getLHS() const {return(LHS_);}
337
338 //! Returns pointer to current RHS.
getRHS() const339 RCP<SerialDenseMatrix<OrdinalType, ScalarType> > getRHS() const {return(RHS_);}
340
341 //! Returns row dimension of system.
numRows() const342 OrdinalType numRows() const {return(M_);}
343
344 //! Returns column dimension of system.
numCols() const345 OrdinalType numCols() const {return(N_);}
346
347 //! Returns lower bandwidth of system.
numLower() const348 OrdinalType numLower() const {return(KL_);}
349
350 //! Returns upper bandwidth of system.
numUpper() const351 OrdinalType numUpper() const {return(KU_);}
352
353 //! Returns pointer to pivot vector (if factorization has been computed), zero otherwise.
IPIV() const354 std::vector<OrdinalType> IPIV() const {return(IPIV_);}
355
356 //! Returns the 1-Norm of the \e this matrix (returns -1 if not yet computed).
ANORM() const357 MagnitudeType ANORM() const {return(ANORM_);}
358
359 //! Returns the reciprocal of the condition number of the \e this matrix (returns -1 if not yet computed).
RCOND() const360 MagnitudeType RCOND() const {return(RCOND_);}
361
362 //! Ratio of smallest to largest row scale factors for the \e this matrix (returns -1 if not yet computed).
363 /*! If ROWCND() is >= 0.1 and AMAX() is not close to overflow or underflow, then equilibration is not needed.
364 */
ROWCND() const365 MagnitudeType ROWCND() const {return(ROWCND_);}
366
367 //! Ratio of smallest to largest column scale factors for the \e this matrix (returns -1 if not yet computed).
368 /*! If COLCND() is >= 0.1 then equilibration is not needed.
369 */
COLCND() const370 MagnitudeType COLCND() const {return(COLCND_);}
371
372 //! Returns the absolute value of the largest entry of the \e this matrix (returns -1 if not yet computed).
AMAX() const373 MagnitudeType AMAX() const {return(AMAX_);}
374
375 //! Returns a pointer to the forward error estimates computed by LAPACK.
FERR() const376 std::vector<MagnitudeType> FERR() const {return(FERR_);}
377
378 //! Returns a pointer to the backward error estimates computed by LAPACK.
BERR() const379 std::vector<MagnitudeType> BERR() const {return(BERR_);}
380
381 //! Returns a pointer to the row scaling vector used for equilibration.
R() const382 std::vector<MagnitudeType> R() const {return(R_);}
383
384 //! Returns a pointer to the column scale vector used for equilibration.
C() const385 std::vector<MagnitudeType> C() const {return(C_);}
386 //@}
387
388 //! @name I/O methods
389 //@{
390 //! Print service methods; defines behavior of ostream << operator.
391 void Print(std::ostream& os) const;
392 //@}
393 protected:
394
allocateWORK()395 void allocateWORK() { LWORK_ = 3*N_; WORK_.resize( LWORK_ ); return;}
396 void resetMatrix();
397 void resetVectors();
398
399
400 bool equilibrate_;
401 bool shouldEquilibrate_;
402 bool equilibratedA_;
403 bool equilibratedB_;
404 bool transpose_;
405 bool factored_;
406 bool estimateSolutionErrors_;
407 bool solutionErrorsEstimated_;
408 bool solved_;
409 bool reciprocalConditionEstimated_;
410 bool refineSolution_;
411 bool solutionRefined_;
412
413 Teuchos::ETransp TRANS_;
414
415 OrdinalType M_;
416 OrdinalType N_;
417 OrdinalType KL_;
418 OrdinalType KU_;
419 OrdinalType Min_MN_;
420 OrdinalType LDA_;
421 OrdinalType LDAF_;
422 OrdinalType INFO_;
423 OrdinalType LWORK_;
424
425 std::vector<OrdinalType> IPIV_;
426 std::vector<int> IWORK_;
427
428 MagnitudeType ANORM_;
429 MagnitudeType RCOND_;
430 MagnitudeType ROWCND_;
431 MagnitudeType COLCND_;
432 MagnitudeType AMAX_;
433
434 RCP<SerialBandDenseMatrix<OrdinalType, ScalarType> > Matrix_;
435 RCP<SerialDenseMatrix<OrdinalType, ScalarType> > LHS_;
436 RCP<SerialDenseMatrix<OrdinalType, ScalarType> > RHS_;
437 RCP<SerialBandDenseMatrix<OrdinalType, ScalarType> > Factor_;
438
439 ScalarType * A_;
440 ScalarType * AF_;
441 std::vector<MagnitudeType> FERR_;
442 std::vector<MagnitudeType> BERR_;
443 std::vector<ScalarType> WORK_;
444 std::vector<MagnitudeType> R_;
445 std::vector<MagnitudeType> C_;
446 #ifdef HAVE_TEUCHOSNUMERICS_EIGEN
447 Eigen::PartialPivLU<EigenMatrix> lu_;
448 #endif
449
450 private:
451
452 // SerialBandDenseSolver copy constructor (put here because we don't want user access)
453 SerialBandDenseSolver(const SerialBandDenseSolver<OrdinalType, ScalarType>& Source);
454 SerialBandDenseSolver & operator=(const SerialBandDenseSolver<OrdinalType, ScalarType>& Source);
455
456 };
457
458 // Helper traits to distinguish work arrays for real and complex-valued datatypes.
459 using namespace details;
460
461 //=============================================================================
462
463 template<typename OrdinalType, typename ScalarType>
SerialBandDenseSolver()464 SerialBandDenseSolver<OrdinalType,ScalarType>::SerialBandDenseSolver()
465 : CompObject(),
466 equilibrate_(false),
467 shouldEquilibrate_(false),
468 equilibratedA_(false),
469 equilibratedB_(false),
470 transpose_(false),
471 factored_(false),
472 estimateSolutionErrors_(false),
473 solutionErrorsEstimated_(false),
474 solved_(false),
475 reciprocalConditionEstimated_(false),
476 refineSolution_(false),
477 solutionRefined_(false),
478 TRANS_(Teuchos::NO_TRANS),
479 M_(0),
480 N_(0),
481 KL_(0),
482 KU_(0),
483 Min_MN_(0),
484 LDA_(0),
485 LDAF_(0),
486 INFO_(0),
487 LWORK_(0),
488 ANORM_(0.0),
489 RCOND_(ScalarTraits<MagnitudeType>::zero()),
490 ROWCND_(ScalarTraits<MagnitudeType>::zero()),
491 COLCND_(ScalarTraits<MagnitudeType>::zero()),
492 AMAX_(ScalarTraits<MagnitudeType>::zero()),
493 A_(NULL),
494 AF_(NULL)
495 {
496 resetMatrix();
497 }
498 //=============================================================================
499
500 template<typename OrdinalType, typename ScalarType>
~SerialBandDenseSolver()501 SerialBandDenseSolver<OrdinalType,ScalarType>::~SerialBandDenseSolver()
502 {}
503
504 //=============================================================================
505
506 template<typename OrdinalType, typename ScalarType>
resetVectors()507 void SerialBandDenseSolver<OrdinalType,ScalarType>::resetVectors()
508 {
509 LHS_ = Teuchos::null;
510 RHS_ = Teuchos::null;
511 reciprocalConditionEstimated_ = false;
512 solutionRefined_ = false;
513 solved_ = false;
514 solutionErrorsEstimated_ = false;
515 equilibratedB_ = false;
516 }
517 //=============================================================================
518
519 template<typename OrdinalType, typename ScalarType>
resetMatrix()520 void SerialBandDenseSolver<OrdinalType,ScalarType>::resetMatrix()
521 {
522 resetVectors();
523 equilibratedA_ = false;
524 factored_ = false;
525 M_ = 0;
526 N_ = 0;
527 KL_ = 0;
528 KU_ = 0;
529 Min_MN_ = 0;
530 LDA_ = 0;
531 LDAF_ = 0;
532 RCOND_ = -ScalarTraits<MagnitudeType>::one();
533 ROWCND_ = -ScalarTraits<MagnitudeType>::one();
534 COLCND_ = -ScalarTraits<MagnitudeType>::one();
535 AMAX_ = -ScalarTraits<MagnitudeType>::one();
536 A_ = 0;
537 AF_ = 0;
538 INFO_ = 0;
539 LWORK_ = 0;
540 R_.resize(0);
541 C_.resize(0);
542 }
543 //=============================================================================
544
545 template<typename OrdinalType, typename ScalarType>
setMatrix(const RCP<SerialBandDenseMatrix<OrdinalType,ScalarType>> & AB)546 int SerialBandDenseSolver<OrdinalType,ScalarType>::setMatrix(const RCP<SerialBandDenseMatrix<OrdinalType,ScalarType> >& AB)
547 {
548
549 OrdinalType m = AB->numRows();
550 OrdinalType n = AB->numCols();
551 OrdinalType kl = AB->lowerBandwidth();
552 OrdinalType ku = AB->upperBandwidth();
553
554 // Check that the new matrix is consistent.
555 TEUCHOS_TEST_FOR_EXCEPTION(AB->values()==0, std::invalid_argument,
556 "SerialBandDenseSolver<T>::setMatrix: A is an empty SerialBandDenseMatrix<T>!");
557
558 resetMatrix();
559 Matrix_ = AB;
560 Factor_ = AB;
561 M_ = m;
562 N_ = n;
563 KL_ = kl;
564 KU_ = ku-kl;
565 Min_MN_ = TEUCHOS_MIN(M_,N_);
566 LDA_ = AB->stride();
567 LDAF_ = LDA_;
568 A_ = AB->values();
569 AF_ = AB->values();
570
571 return(0);
572 }
573 //=============================================================================
574
575 template<typename OrdinalType, typename ScalarType>
setVectors(const RCP<SerialDenseMatrix<OrdinalType,ScalarType>> & X,const RCP<SerialDenseMatrix<OrdinalType,ScalarType>> & B)576 int SerialBandDenseSolver<OrdinalType,ScalarType>::setVectors(const RCP<SerialDenseMatrix<OrdinalType,ScalarType> >& X,
577 const RCP<SerialDenseMatrix<OrdinalType,ScalarType> >& B)
578 {
579 // Check that these new vectors are consistent.
580 TEUCHOS_TEST_FOR_EXCEPTION(B->numRows()!=X->numRows() || B->numCols() != X->numCols(), std::invalid_argument,
581 "SerialBandDenseSolver<T>::setVectors: X and B are not the same size!");
582 TEUCHOS_TEST_FOR_EXCEPTION(B->values()==0, std::invalid_argument,
583 "SerialBandDenseSolver<T>::setVectors: B is an empty SerialDenseMatrix<T>!");
584 TEUCHOS_TEST_FOR_EXCEPTION(X->values()==0, std::invalid_argument,
585 "SerialBandDenseSolver<T>::setVectors: X is an empty SerialDenseMatrix<T>!");
586 TEUCHOS_TEST_FOR_EXCEPTION(B->stride()<1, std::invalid_argument,
587 "SerialBandDenseSolver<T>::setVectors: B has an invalid stride!");
588 TEUCHOS_TEST_FOR_EXCEPTION(X->stride()<1, std::invalid_argument,
589 "SerialBandDenseSolver<T>::setVectors: X has an invalid stride!");
590
591 resetVectors();
592 LHS_ = X;
593 RHS_ = B;
594 return(0);
595 }
596 //=============================================================================
597
598 template<typename OrdinalType, typename ScalarType>
estimateSolutionErrors(bool flag)599 void SerialBandDenseSolver<OrdinalType,ScalarType>::estimateSolutionErrors(bool flag)
600 {
601 estimateSolutionErrors_ = flag;
602
603 // If the errors are estimated, this implies that the solution must be refined
604 refineSolution_ = refineSolution_ || flag;
605 }
606 //=============================================================================
607
608 template<typename OrdinalType, typename ScalarType>
factor()609 int SerialBandDenseSolver<OrdinalType,ScalarType>::factor() {
610
611 if (factored()) return(0); // Already factored
612
613 ANORM_ = Matrix_->normOne(); // Compute 1-Norm of A
614
615 // If we want to refine the solution, then the factor must
616 // be stored separatedly from the original matrix
617
618 if (A_ == AF_)
619 if (refineSolution_ ) {
620 Factor_ = rcp( new SerialBandDenseMatrix<OrdinalType,ScalarType>(*Matrix_) );
621 AF_ = Factor_->values();
622 LDAF_ = Factor_->stride();
623 }
624
625 int ierr = 0;
626 if (equilibrate_) ierr = equilibrateMatrix();
627
628 if (ierr!=0) return(ierr);
629
630 if (IPIV_.size()==0) IPIV_.resize( N_ ); // Allocated Pivot vector if not already done.
631
632 INFO_ = 0;
633
634 #ifdef HAVE_TEUCHOSNUMERICS_EIGEN
635 EigenMatrixMap aMap( AF_+KL_, KL_+KU_+1, N_, EigenOuterStride(LDAF_) );
636 EigenMatrix fullMatrix(M_, N_);
637 for (OrdinalType j=0; j<N_; j++) {
638 for (OrdinalType i=TEUCHOS_MAX(0,j-KU_); i<=TEUCHOS_MIN(M_-1, j+KL_); i++) {
639 fullMatrix(i,j) = aMap(KU_+i-j, j);
640 }
641 }
642
643 EigenPermutationMatrix p;
644 EigenOrdinalArray ind;
645 EigenOrdinalArrayMap ipivMap( &IPIV_[0], Min_MN_ );
646
647 lu_.compute(fullMatrix);
648 p = lu_.permutationP();
649 ind = p.indices();
650
651 for (OrdinalType i=0; i<ipivMap.innerSize(); i++) {
652 ipivMap(i) = ind(i);
653 }
654
655 #else
656 this->GBTRF(M_, N_, KL_, KU_, AF_, LDAF_, &IPIV_[0], &INFO_);
657 #endif
658
659 factored_ = true;
660
661 return(INFO_);
662 }
663
664 //=============================================================================
665
666 template<typename OrdinalType, typename ScalarType>
solve()667 int SerialBandDenseSolver<OrdinalType,ScalarType>::solve() {
668
669 // If the user want the matrix to be equilibrated or wants a refined solution, we will
670 // call the X interface.
671 // Otherwise, if the matrix is already factored we will call the TRS interface.
672 // Otherwise, if the matrix is unfactored we will call the SV interface.
673
674 int ierr = 0;
675 if (equilibrate_) {
676 ierr = equilibrateRHS();
677 }
678 if (ierr != 0) return(ierr); // Can't equilibrate B, so return.
679
680 TEUCHOS_TEST_FOR_EXCEPTION( RHS_==Teuchos::null, std::invalid_argument,
681 "SerialBandDenseSolver<T>::solve: No right-hand side vector (RHS) has been set for the linear system!");
682 TEUCHOS_TEST_FOR_EXCEPTION( LHS_==Teuchos::null, std::invalid_argument,
683 "SerialBandDenseSolver<T>::solve: No solution vector (LHS) has been set for the linear system!");
684
685 if (!factored()) factor(); // Matrix must be factored
686
687 TEUCHOS_TEST_FOR_EXCEPTION( (equilibratedA_ && !equilibratedB_) || (!equilibratedA_ && equilibratedB_) ,
688 std::logic_error, "SerialBandDenseSolver<T>::solve: Matrix and vectors must be similarly scaled!");
689
690 if (shouldEquilibrate() && !equilibratedA_)
691 std::cout << "WARNING! SerialBandDenseSolver<T>::solve: System should be equilibrated!" << std::endl;
692
693 if (RHS_->values()!=LHS_->values()) {
694 (*LHS_) = (*RHS_); // Copy B to X if needed
695 }
696 INFO_ = 0;
697
698 #ifdef HAVE_TEUCHOSNUMERICS_EIGEN
699 EigenMatrixMap rhsMap( RHS_->values(), RHS_->numRows(), RHS_->numCols(), EigenOuterStride(RHS_->stride()) );
700 EigenMatrixMap lhsMap( LHS_->values(), LHS_->numRows(), LHS_->numCols(), EigenOuterStride(LHS_->stride()) );
701 if ( TRANS_ == Teuchos::NO_TRANS ) {
702 lhsMap = lu_.solve(rhsMap);
703 } else if ( TRANS_ == Teuchos::TRANS ) {
704 lu_.matrixLU().template triangularView<Eigen::Upper>().transpose().solveInPlace(lhsMap);
705 lu_.matrixLU().template triangularView<Eigen::UnitLower>().transpose().solveInPlace(lhsMap);
706 EigenMatrix x = lu_.permutationP().transpose() * lhsMap;
707 lhsMap = x;
708 } else {
709 lu_.matrixLU().template triangularView<Eigen::Upper>().adjoint().solveInPlace(lhsMap);
710 lu_.matrixLU().template triangularView<Eigen::UnitLower>().adjoint().solveInPlace(lhsMap);
711 EigenMatrix x = lu_.permutationP().transpose() * lhsMap;
712 lhsMap = x;
713 }
714 #else
715 this->GBTRS(ETranspChar[TRANS_], N_, KL_, KU_, RHS_->numCols(), AF_, LDAF_, &IPIV_[0], LHS_->values(), LHS_->stride(), &INFO_);
716 #endif
717
718 if (INFO_!=0) return(INFO_);
719 solved_ = true;
720
721 int ierr1=0;
722 if (refineSolution_) ierr1 = applyRefinement();
723 if (ierr1!=0)
724 return(ierr1);
725
726 if (equilibrate_) ierr1 = unequilibrateLHS();
727
728 return(ierr1);
729 }
730 //=============================================================================
731
732 template<typename OrdinalType, typename ScalarType>
applyRefinement()733 int SerialBandDenseSolver<OrdinalType,ScalarType>::applyRefinement()
734 {
735 TEUCHOS_TEST_FOR_EXCEPTION(!solved(), std::logic_error,
736 "SerialBandDenseSolver<T>::applyRefinement: Must have an existing solution!");
737 TEUCHOS_TEST_FOR_EXCEPTION(A_==AF_, std::logic_error,
738 "SerialBandDenseSolver<T>::applyRefinement: Cannot apply refinement if no original copy of A!");
739
740 #ifdef HAVE_TEUCHOSNUMERICS_EIGEN
741 // Implement templated GERFS or use Eigen.
742 return (-1);
743 #else
744
745 OrdinalType NRHS = RHS_->numCols();
746 FERR_.resize( NRHS );
747 BERR_.resize( NRHS );
748 allocateWORK();
749
750 INFO_ = 0;
751 std::vector<typename details::lapack_traits<ScalarType>::iwork_type> GBRFS_WORK( N_ );
752 this->GBRFS(ETranspChar[TRANS_], N_, KL_, KU_, NRHS, A_+KL_, LDA_, AF_, LDAF_, &IPIV_[0],
753 RHS_->values(), RHS_->stride(), LHS_->values(), LHS_->stride(),
754 &FERR_[0], &BERR_[0], &WORK_[0], &GBRFS_WORK[0], &INFO_);
755
756 solutionErrorsEstimated_ = true;
757 reciprocalConditionEstimated_ = true;
758 solutionRefined_ = true;
759
760 return(INFO_);
761 #endif
762 }
763
764 //=============================================================================
765
766 template<typename OrdinalType, typename ScalarType>
computeEquilibrateScaling()767 int SerialBandDenseSolver<OrdinalType,ScalarType>::computeEquilibrateScaling()
768 {
769 if (R_.size()!=0) return(0); // Already computed
770
771 R_.resize( M_ );
772 C_.resize( N_ );
773
774 INFO_ = 0;
775 this->GBEQU (M_, N_, KL_, KU_, AF_+KL_, LDAF_, &R_[0], &C_[0], &ROWCND_, &COLCND_, &AMAX_, &INFO_);
776
777 if (COLCND_<0.1*ScalarTraits<MagnitudeType>::one() ||
778 ROWCND_<0.1*ScalarTraits<MagnitudeType>::one() ||
779 AMAX_ < ScalarTraits<ScalarType>::rmin() || AMAX_ > ScalarTraits<ScalarType>::rmax() )
780 shouldEquilibrate_ = true;
781
782 return(INFO_);
783 }
784
785 //=============================================================================
786
787 template<typename OrdinalType, typename ScalarType>
equilibrateMatrix()788 int SerialBandDenseSolver<OrdinalType,ScalarType>::equilibrateMatrix()
789 {
790 OrdinalType i, j;
791 int ierr = 0;
792
793 if (equilibratedA_) return(0); // Already done.
794 if (R_.size()==0) ierr = computeEquilibrateScaling(); // Compute R and C if needed.
795 if (ierr!=0) return(ierr); // If return value is not zero, then can't equilibrate.
796
797 if (A_==AF_) {
798
799 ScalarType * ptr;
800 for (j=0; j<N_; j++) {
801 ptr = A_ + KL_ + j*LDA_ + TEUCHOS_MAX(KU_-j, 0);
802 ScalarType s1 = C_[j];
803 for (i=TEUCHOS_MAX(0,j-KU_); i<=TEUCHOS_MIN(M_-1,j+KL_); i++) {
804 *ptr = *ptr*s1*R_[i];
805 ptr++;
806 }
807 }
808 } else {
809
810 ScalarType * ptr;
811 ScalarType * ptrL;
812 ScalarType * ptrU;
813 for (j=0; j<N_; j++) {
814 ptr = A_ + KL_ + j*LDA_ + TEUCHOS_MAX(KU_-j, 0);
815 ScalarType s1 = C_[j];
816 for (i=TEUCHOS_MAX(0,j-KU_); i<=TEUCHOS_MIN(M_-1,j+KL_); i++) {
817 *ptr = *ptr*s1*R_[i];
818 ptr++;
819 }
820 }
821 for (j=0; j<N_; j++) {
822 ptrU = AF_ + j*LDAF_ + TEUCHOS_MAX(KL_+KU_-j, 0);
823 ptrL = AF_ + KL_ + KU_ + 1 + j*LDAF_;
824 ScalarType s1 = C_[j];
825 for (i=TEUCHOS_MAX(0,j-(KL_+KU_)); i<=TEUCHOS_MIN(M_-1,j); i++) {
826 *ptrU = *ptrU*s1*R_[i];
827 ptrU++;
828 }
829 for (i=TEUCHOS_MAX(0,j); i<=TEUCHOS_MIN(M_-1,j+KL_); i++) {
830 *ptrL = *ptrL*s1*R_[i];
831 ptrL++;
832 }
833 }
834 }
835
836 equilibratedA_ = true;
837
838 return(0);
839 }
840
841 //=============================================================================
842
843 template<typename OrdinalType, typename ScalarType>
equilibrateRHS()844 int SerialBandDenseSolver<OrdinalType,ScalarType>::equilibrateRHS()
845 {
846 OrdinalType i, j;
847 int ierr = 0;
848
849 if (equilibratedB_) return(0); // Already done.
850 if (R_.size()==0) ierr = computeEquilibrateScaling(); // Compute R and C if needed.
851 if (ierr!=0) return(ierr); // Can't count on R and C being computed.
852
853 MagnitudeType * R_tmp = &R_[0];
854 if (transpose_) R_tmp = &C_[0];
855
856 OrdinalType LDB = RHS_->stride(), NRHS = RHS_->numCols();
857 ScalarType * B = RHS_->values();
858 ScalarType * ptr;
859 for (j=0; j<NRHS; j++) {
860 ptr = B + j*LDB;
861 for (i=0; i<M_; i++) {
862 *ptr = *ptr*R_tmp[i];
863 ptr++;
864 }
865 }
866
867 equilibratedB_ = true;
868
869 return(0);
870 }
871
872
873 //=============================================================================
874
875 template<typename OrdinalType, typename ScalarType>
unequilibrateLHS()876 int SerialBandDenseSolver<OrdinalType,ScalarType>::unequilibrateLHS()
877 {
878 OrdinalType i, j;
879
880 if (!equilibratedB_) return(0); // Nothing to do
881
882 MagnitudeType * C_tmp = &C_[0];
883 if (transpose_) C_tmp = &R_[0];
884
885 OrdinalType LDX = LHS_->stride(), NLHS = LHS_->numCols();
886 ScalarType * X = LHS_->values();
887 ScalarType * ptr;
888 for (j=0; j<NLHS; j++) {
889 ptr = X + j*LDX;
890 for (i=0; i<N_; i++) {
891 *ptr = *ptr*C_tmp[i];
892 ptr++;
893 }
894 }
895
896 return(0);
897 }
898
899 //=============================================================================
900
901 template<typename OrdinalType, typename ScalarType>
reciprocalConditionEstimate(MagnitudeType & Value)902 int SerialBandDenseSolver<OrdinalType,ScalarType>::reciprocalConditionEstimate(MagnitudeType & Value)
903 {
904 #ifdef HAVE_TEUCHOSNUMERICS_EIGEN
905 // Implement templated GECON or use Eigen. Eigen currently doesn't have a condition estimation function.
906 return (-1);
907 #else
908
909 if (reciprocalConditionEstimated()) {
910 Value = RCOND_;
911 return(0); // Already computed, just return it.
912 }
913
914 if ( ANORM_<ScalarTraits<MagnitudeType>::zero() ) ANORM_ = Matrix_->normOne();
915
916 int ierr = 0;
917 if (!factored()) ierr = factor(); // Need matrix factored.
918 if (ierr!=0) return(ierr);
919
920 allocateWORK();
921
922 // We will assume a one-norm condition number
923 INFO_ = 0;
924 std::vector<typename details::lapack_traits<ScalarType>::iwork_type> GBCON_WORK( N_ );
925 this->GBCON( '1', N_, KL_, KU_, AF_, LDAF_, &IPIV_[0], ANORM_, &RCOND_, &WORK_[0], &GBCON_WORK[0], &INFO_);
926 reciprocalConditionEstimated_ = true;
927 Value = RCOND_;
928
929 return(INFO_);
930 #endif
931 }
932 //=============================================================================
933
934 template<typename OrdinalType, typename ScalarType>
Print(std::ostream & os) const935 void SerialBandDenseSolver<OrdinalType,ScalarType>::Print(std::ostream& os) const {
936
937 if (Matrix_!=Teuchos::null) os << "Solver Matrix" << std::endl << *Matrix_ << std::endl;
938 if (Factor_!=Teuchos::null) os << "Solver Factored Matrix" << std::endl << *Factor_ << std::endl;
939 if (LHS_ !=Teuchos::null) os << "Solver LHS" << std::endl << *LHS_ << std::endl;
940 if (RHS_ !=Teuchos::null) os << "Solver RHS" << std::endl << *RHS_ << std::endl;
941
942 }
943
944 //=============================================================================
945
946
947 //=============================================================================
948
949
950 } // namespace Teuchos
951
952 #endif /* _TEUCHOS_SERIALBANDDENSESOLVER_HPP_ */
953