ShyLU  Version of the Day
shylu_local_schur_operator.cpp
1 
2 //@HEADER
3 // ************************************************************************
4 //
5 // ShyLU: Hybrid preconditioner package
6 // Copyright 2012 Sandia Corporation
7 //
8 // Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
9 // the U.S. Government retains certain rights in this software.
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 Michael A. Heroux (maherou@sandia.gov)
39 //
40 // ************************************************************************
41 //@HEADER
42 
43 
44 #include <Epetra_Operator.h>
45 #include <Epetra_CrsMatrix.h>
46 #include <Epetra_LinearProblem.h>
47 #include <Amesos_BaseSolver.h>
48 #include <Epetra_MultiVector.h>
49 #include <Teuchos_Time.hpp>
50 #include "shylu_local_schur_operator.h"
51 #include "shylu_util.h"
52 
53 
54 ShyLU_Local_Schur_Operator::ShyLU_Local_Schur_Operator(
55  shylu_config *config,
56  shylu_symbolic *ssym, // symbolic structure
57  Epetra_CrsMatrix *G,
58  Epetra_CrsMatrix *R,
59  Epetra_LinearProblem *LP, Amesos_BaseSolver *solver,
60  Ifpack_Preconditioner *ifSolver, Epetra_CrsMatrix *C,
61  Epetra_Map *localDRowMap, int nvectors)
62 {
63  ssym_ = ssym;
64  config_ = config;
65  G_ = G;
66  R_ = R;
67  LP_ = LP;
68  solver_ = solver;
69  ifSolver_ = ifSolver;
70  C_ = C;
71  localDRowMap_ = localDRowMap;
72  orig_lhs_ = ssym_->OrigLP->GetLHS();
73  orig_rhs_ = ssym_->OrigLP->GetRHS();
74  ResetTempVectors(nvectors);
75  nvectors_ = nvectors;
76 
77  cntApply = 0;
78 
79 #ifdef TIMING_OUTPUT
80  using Teuchos::RCP;
81  using Teuchos::Time;
82  matvec_time_ = RCP<Time>(new Time("First 2 Matvec time in prober: "));
83  localize_time_ = RCP<Time>(new Time("Localize time in prober: "));
84  trisolve_time_ = RCP<Time>(new Time("Triangular solve time in prober: "));
85  dist_time_ = RCP<Time>(new Time("Distribute time in prober: "));
86  matvec2_time_ = RCP<Time>(new Time("Second 1 matvec time in prober: "));
87  apply_time_ = RCP<Time>(new Time("Apply time in prober: "));
88  update_time_ = RCP<Time>(new Time("Update time in prober: "));
89 #endif
90 }
91 
92 int ShyLU_Local_Schur_Operator::SetUseTranspose(bool useTranspose)
93 {
94  return 0;
95 }
96 
97 int ShyLU_Local_Schur_Operator::Apply(const Epetra_MultiVector &X,
98  Epetra_MultiVector &Y) const
99 {
100 #ifdef TIMING_OUTPUT
101  apply_time_->start();
102 #endif
103 
104  //int nvectors = X.NumVectors();
105  // TODO: For local_scur_operation this will always be local !!!
106  // Remove these cases or merge with schur operator.
107  //bool local = (G_->Comm().NumProc() == 1);
108  int err;
109 
110 #ifdef TIMING_OUTPUT
111  matvec_time_->start();
112 #endif
113 
114  err = G_->Multiply(false, X, *temp2);
115  assert((err == 0));
116  err = C_->Multiply(false, X, *temp);
117  assert((err == 0));
118 
119 #ifdef TIMING_OUTPUT
120  matvec_time_->stop();
121 #endif
122 
123 
124 #ifdef DEBUG
125  int nrows = C_->RowMap().NumMyElements();
126  cout << "DEBUG MODE" << endl;
127  ASSERT((nrows == localDRowMap_->NumGlobalElements()));
128 
129  int gids[nrows], gids1[nrows];
130  C_->RowMap().MyGlobalElements(gids);
131  localDRowMap_->MyGlobalElements(gids1);
132 
133  for (int i = 0; i < nrows; i++)
134  {
135  ASSERT((gids[i] == gids1[i]));
136  }
137 #endif
138 
139 #ifdef TIMING_OUTPUT
140  localize_time_->start();
141 #endif
142 
143 #ifdef TIMING_OUTPUT
144  localize_time_->stop();
145  trisolve_time_->start();
146 #endif
147 
148  // The lhs and rhs is set in ResetTempVectors()
149  if (config_->amesosForDiagonal)
150  solver_->Solve();
151  else
152  ifSolver_->ApplyInverse(*temp, *localX);
153 
154 #ifdef TIMING_OUTPUT
155  trisolve_time_->stop();
156  dist_time_->start();
157 #endif
158 
159 #ifdef TIMING_OUTPUT
160  dist_time_->stop();
161  matvec2_time_->start();
162 #endif
163 
164  R_->Multiply(false, *localX, Y);
165 
166 #ifdef TIMING_OUTPUT
167  matvec2_time_->stop();
168  update_time_->start();
169 #endif
170 
171  err = Y.Update(1.0, *temp2, -1.0);
172  //cout << Y.MyLength() << " " << temp2.MyLength() << endl;
173  assert((err == 0));
174 
175 #ifdef TIMING_OUTPUT
176  update_time_->stop();
177  apply_time_->stop();
178 #endif
179  //cout << "Out of local schur's Apply" << endl;
180  cntApply++;
181  return 0;
182 }
183 
184 void ShyLU_Local_Schur_Operator::PrintTimingInfo()
185 {
186 #ifdef TIMING_OUTPUT
187  cout << matvec_time_->name()<< matvec_time_->totalElapsedTime() << endl;
188  cout << localize_time_->name()<< localize_time_->totalElapsedTime() << endl;
189  cout << trisolve_time_->name()<< trisolve_time_->totalElapsedTime() << endl;
190  cout << dist_time_->name() <<dist_time_->totalElapsedTime() << endl;
191  cout << matvec2_time_->name() <<matvec2_time_->totalElapsedTime() << endl;
192  cout << update_time_->name() <<update_time_->totalElapsedTime() << endl;
193  cout << apply_time_->name() <<apply_time_->totalElapsedTime() << endl;
194 #endif
195 }
196 
197 void ShyLU_Local_Schur_Operator::ResetTempVectors(int nvectors)
198 {
199  using Teuchos::RCP;
200  nvectors_ = nvectors;
201  // If vectors were created already, they will be freed.
203  (new Epetra_MultiVector(View, *(ssym_->Drhs), 0, nvectors));
205  (new Epetra_MultiVector(View, *(ssym_->Dlhs), 0, nvectors));
207  (new Epetra_MultiVector(View, *(ssym_->Gvec), 0, nvectors));
208  ssym_->OrigLP->SetLHS(localX.getRawPtr());
209  ssym_->OrigLP->SetRHS(temp.getRawPtr());
210  ssym_->ReIdx_LP->fwd();
211 
212 }
213 
214 int ShyLU_Local_Schur_Operator::ApplyInverse(const Epetra_MultiVector &X,
215  Epetra_MultiVector &Y) const
216 {
217  return 0;
218 }
219 
220 double ShyLU_Local_Schur_Operator::NormInf() const
221 {
222  return -1;
223 }
224 
225 const char* ShyLU_Local_Schur_Operator::Label() const
226 {
227  return "Shylu Local Schur Operator (Wide Separator)";
228 }
229 
230 bool ShyLU_Local_Schur_Operator::UseTranspose() const
231 {
232  return false;
233 }
234 
235 bool ShyLU_Local_Schur_Operator::HasNormInf() const
236 {
237  return false;
238 }
239 
240 const Epetra_Comm& ShyLU_Local_Schur_Operator::Comm() const
241 {
242  return G_->Comm();
243 }
244 
245 const Epetra_Map& ShyLU_Local_Schur_Operator::OperatorDomainMap() const
246 {
247  //return C_->ColMap();
248  return G_->DomainMap();
249 }
250 
251 const Epetra_Map& ShyLU_Local_Schur_Operator::OperatorRangeMap() const
252 {
253  //return R_->RowMap();
254  return G_->RangeMap();
255 }
Utilities for ShyLU.
T * getRawPtr() const