Stokhos Package Browser (Single Doxygen Collection) Version of the Day
Loading...
Searching...
No Matches
Stokhos_MLPrecOp.cpp
Go to the documentation of this file.
1// $Id$
2// $Source$
3// @HEADER
4// ***********************************************************************
5//
6// Stokhos Package
7// Copyright (2009) Sandia Corporation
8//
9// Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
10// license for use of this work by or on behalf of the U.S. Government.
11//
12// Redistribution and use in source and binary forms, with or without
13// modification, are permitted provided that the following conditions are
14// met:
15//
16// 1. Redistributions of source code must retain the above copyright
17// notice, this list of conditions and the following disclaimer.
18//
19// 2. Redistributions in binary form must reproduce the above copyright
20// notice, this list of conditions and the following disclaimer in the
21// documentation and/or other materials provided with the distribution.
22//
23// 3. Neither the name of the Corporation nor the names of the
24// contributors may be used to endorse or promote products derived from
25// this software without specific prior written permission.
26//
27// THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
28// EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
29// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
30// PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
31// CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
32// EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
33// PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
34// PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
35// LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
36// NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
37// SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
38//
39// Questions? Contact Eric T. Phipps (etphipp@sandia.gov).
40//
41// ***********************************************************************
42// @HEADER
43
44#include "Stokhos_MLPrecOp.hpp"
45
46#ifdef HAVE_STOKHOS_ML
47
48//==============================================================================
49// constructor -- it's presumed that the user has constructed the ML
50// object somewhere else. Uses AMG to invert the mean stiffness matrix.
51Stokhos::MLPrecOp::MLPrecOp(const Epetra_CrsMatrix& mean_op,const Teuchos::Array<double>& norms, const Epetra_Comm& Comm, const Epetra_Map& DMap,const Epetra_Map& RMap)
52 : //Epetra_Operator(),
53 Comm_(Comm),
54 Label_(0),
55 DomainMap_(DMap),
56 RangeMap_(RMap),
57 norms_(norms)
58 {
59 Label_ = "Stochos::MLPrecOp";
60
61
62 // create a parameter list for ML options
63 Teuchos::ParameterList MLList;
64
65 // Sets default parameters for classic smoothed aggregation. After this
66 // call, MLList contains the default values for the ML parameters,
67 // as required by typical smoothed aggregation for symmetric systems.
68 // Other sets of parameters are available for non-symmetric systems
69 // ("DD" and "DD-ML"), and for the Maxwell equations ("maxwell").
70 ML_Epetra::SetDefaults("SA",MLList);
71
72 // overwrite some parameters. Please refer to the user's guide
73 // for more information
74 // some of the parameters do not differ from their default value,
75 // and they are here reported for the sake of clarity
76
77 // output level, 0 being silent and 10 verbose
78 MLList.set("ML output", 10);
79 // maximum number of levels
80 MLList.set("max levels",5);
81 // set finest level to 0
82 MLList.set("increasing or decreasing","increasing");
83
84 // use Uncoupled scheme to create the aggregate
85 MLList.set("aggregation: type", "Uncoupled");
86
87 // smoother is Chebyshev. Example file
88 // `ml/examples/TwoLevelDD/ml_2level_DD.cpp' shows how to use
89 // AZTEC's preconditioners as smoothers
90
91 MLList.set("smoother: type","Chebyshev");
92 MLList.set("smoother: sweeps",3);
93
94 // use both pre and post smoothing
95 MLList.set("smoother: pre or post", "both");
96 //MLList.set("coarse: max size", );
97
98#ifdef HAVE_ML_AMESOS
99 // solve with serial direct solver KLU
100 MLList.set("coarse: type","Amesos-KLU");
101#else
102 // this is for testing purposes only, you should have
103 // a direct solver for the coarse problem (either Amesos, or the SuperLU/
104 // SuperLU_DIST interface of ML)
105 MLList.set("coarse: type","Jacobi");
106#endif
107
108 // Creates the preconditioning object. We suggest to use `new' and
109 // `delete' because the destructor contains some calls to MPI (as
110 // required by ML and possibly Amesos). This is an issue only if the
111 // destructor is called **after** MPI_Finalize().
112
113 MLPrec = new ML_Epetra::MultiLevelPreconditioner(mean_op, MLList);
114
115 // verify unused parameters on process 0 (put -1 to print on all
116 // processes)
117 MLPrec->PrintUnused(0);
118
119
120}
121
122//==============================================================================
123//Applies the preconditioner implicitly to the global stiffness matrix.
124//==============================================================================
125
126int Stokhos::MLPrecOp::ApplyInverse(const Epetra_MultiVector& X, Epetra_MultiVector& Y) const{
127
128 if (!X.Map().SameAs(OperatorDomainMap()))
129 std::cout << "!X.Map().SameAs(OperatorDomainMap())\n";
130 if (!Y.Map().SameAs(OperatorRangeMap()))
131 std::cout<< "!Y.Map().SameAs(OperatorRangeMap())\n";
132 if (Y.NumVectors()!=X.NumVectors())
133 std::cout<< "Y.NumVectors()!=X.NumVectors()\n";
134
135
136
137 for(int mm = 0; mm< X.NumVectors(); mm++){
138
139
140 int N_xi = norms_.size();
141 int N_x = X.MyLength()/N_xi;
142
143 // Construct a Map with NumElements and index base of 0
144 //Epetra_Map::Epetra_Map Map(N_x, 0, Comm_);
145 Epetra_Map Map(MLPrec->OperatorDomainMap());
146
147 // Form x and y into block vectors.
148 Epetra_MultiVector xBlock(Map,N_xi);
149 Epetra_MultiVector yBlock(MLPrec->OperatorRangeMap(),N_xi);
150
151 // get the local size of the vector
152 int MyLength = xBlock.MyLength();
153 for( int c=0; c<N_xi ; c++){
154 for( int i=0; i<MyLength; i++){
155 xBlock[c][i] = (X)[mm][c*N_x + i];
156 yBlock[c][i] = 0;
157 }
158 }
159
160 Epetra_MultiVector blockProducts(MLPrec->OperatorRangeMap(),N_xi);
161 MLPrec->ApplyInverse(xBlock, blockProducts);
162
163
164 for(int j = 0; j<N_xi; j++){
165 (*yBlock(j)).Update(1/norms_[j],*blockProducts(j), 1.0);
166 }
167
168 for( int c=0; c<N_xi ; c++){
169 for( int i=0; i<MyLength; i++){
170 (Y)[mm][c*N_x + i] = yBlock[c][i];
171 }
172 }
173 }
174 return 1;
175}
176
177#endif
178
179
180
181
const Epetra_BlockMap & Map() const
int NumVectors() const
int MyLength() const