Stokhos Package Browser (Single Doxygen Collection) Version of the Day
Loading...
Searching...
No Matches
Stokhos_GSPreconditioner.hpp
Go to the documentation of this file.
1/*
2// @HEADER
3// ***********************************************************************
4//
5// Stokhos Package
6// Copyright (2009) Sandia Corporation
7//
8// Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
9// license for use of this work by or on behalf of the U.S. Government.
10//
11// Redistribution and use in source and binary forms, with or without
12// modification, are permitted provided that the following conditions are
13// met:
14//
15// 1. Redistributions of source code must retain the above copyright
16// notice, this list of conditions and the following disclaimer.
17//
18// 2. Redistributions in binary form must reproduce the above copyright
19// notice, this list of conditions and the following disclaimer in the
20// documentation and/or other materials provided with the distribution.
21//
22// 3. Neither the name of the Corporation nor the names of the
23// contributors may be used to endorse or promote products derived from
24// this software without specific prior written permission.
25//
26// THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
27// EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
28// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
29// PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
30// CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
31// EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
32// PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
33// PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
34// LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
35// NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
36// SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
37//
38// Questions? Contact Eric T. Phipps (etphipp@sandia.gov).
39//
40// ***********************************************************************
41// @HEADER
42*/
43
44
45#ifndef STOKHOS_GSPRECONDITIONER_HPP
46#define STOKHOS_GSPRECONDITIONER_HPP
47
48#include "Teuchos_RCP.hpp"
49#include "Stokhos_Operator.hpp"
50#include "Teuchos_SerialDenseMatrix.hpp"
51#include "Teuchos_LAPACK.hpp"
52
53namespace Stokhos {
54
55 template <typename ordinal_type, typename value_type>
57 public Stokhos::Operator<ordinal_type,value_type> {
58 public:
59
62 const Teuchos::SerialDenseMatrix<ordinal_type,value_type> & A_,
63 const ordinal_type sym_) : A(A_), sym(sym_) {}
64
66 virtual ~GSPreconditioner() {}
67
68 // m is the number of GS iterations to solve Mz=r
69 // If sym=0 then do symmetric Gauss Seidel
70 virtual ordinal_type ApplyInverse(
71 const Teuchos::SerialDenseMatrix<ordinal_type, value_type>& Input,
72 Teuchos::SerialDenseMatrix<ordinal_type, value_type>& Result,
73 ordinal_type m) const {
74 Result.assign(Input);
75 ordinal_type n=A.numRows();
76 ordinal_type info;
78
79 //Get lower triangular part of A
80 Teuchos::SerialDenseMatrix<ordinal_type, value_type> L(A);
81 for (ordinal_type i=0; i<n; i++){
82 for (ordinal_type j=0; j<n; j++){
83 if (j>i)
84 L(i,j)=0;
85 }
86 }
87
88 if (sym==0){
89 //Get inv(diag(A))=D
90 Teuchos::SerialDenseMatrix<ordinal_type, value_type> D(n,n);
91 for (ordinal_type i=0; i<n; i++){
92 for (ordinal_type j=0; j<n; j++){
93 D(i,i)=1/A(i,i);
94 }
95 }
96
97 //Get upper triangular part of A
98 Teuchos::SerialDenseMatrix<ordinal_type, value_type> U(A);
99 for (ordinal_type i=0; i<n; i++){
100 for (ordinal_type j=0; j<n; j++){
101 if (i>j)
102 U(i,j)=0;
103 }
104 }
105
106 Result.assign(Input);
107
108 //compute M=(L+D)inv(diagA)
109 Teuchos::SerialDenseMatrix<ordinal_type, value_type> M(n,n);
110 M.multiply(Teuchos::NO_TRANS,Teuchos::NO_TRANS,1.0, L, D, 0.0);
111
112 //Forward solve Lz=r
113 lapack.TRTRS('L', 'N', 'N', M.numRows(), 1, M.values(), M.stride(), Result.values(), Result.stride(),&info);
114
115 //Backward solve Uw=z
116 lapack.TRTRS('U', 'N', 'N', U.numRows(), 1, U.values(), U.stride(), Result.values(), Result.stride(),&info);
117
118 }
119 else{
120 lapack.TRTRS('L', 'N', 'N', L.numRows(), 1, L.values(), L.stride(), Result.values(), Result.stride(),&info);
121 }
122
123 return 0;
124 }
125
126 protected:
127 const Teuchos::SerialDenseMatrix<ordinal_type,value_type> & A;
128
129 const ordinal_type sym;
130 }; // class GSPreconditioner
131
132} // namespace Stokhos
133
134#endif // STOKHOS_GSPRECONDITIONER_HPP
135
virtual ~GSPreconditioner()
Destructor.
GSPreconditioner(const Teuchos::SerialDenseMatrix< ordinal_type, value_type > &A_, const ordinal_type sym_)
Constructor.
const Teuchos::SerialDenseMatrix< ordinal_type, value_type > & A
virtual ordinal_type ApplyInverse(const Teuchos::SerialDenseMatrix< ordinal_type, value_type > &Input, Teuchos::SerialDenseMatrix< ordinal_type, value_type > &Result, ordinal_type m) const
Returns the result of a Operator inverse applied to a Teuchos::SerialDenseMatrix Input in Result.
Specialization for Sacado::UQ::PCE< Storage<...> >
Top-level namespace for Stokhos classes and functions.