ShyLU  Version of the Day
shylu_probing_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 <Ifpack_Preconditioner.h>
49 #include <Epetra_MultiVector.h>
50 #include <Teuchos_Time.hpp>
51 #include "shylu_probing_operator.h"
52 
53 // TODO: 1. ltemp is not needed in the all local case.
54 
55 ShyLU_Probing_Operator::ShyLU_Probing_Operator(
56  shylu_config *config,
57  shylu_symbolic *ssym, // symbolic structure
58  Epetra_CrsMatrix *G,
59  Epetra_CrsMatrix *R,
60  Epetra_LinearProblem *LP, Amesos_BaseSolver *solver,
61  Ifpack_Preconditioner *ifSolver, Epetra_CrsMatrix *C,
62  Epetra_Map *localDRowMap, int nvectors)
63 {
64  config_ = config;
65  ssym_ = ssym;
66  G_ = G;
67  R_ = R;
68  LP_ = LP;
69  solver_ = solver;
70  ifSolver_ = ifSolver;
71  C_ = C;
72  localDRowMap_ = localDRowMap;
73  ResetTempVectors(nvectors);
74  nvectors_ = nvectors;
75 
76  cntApply = 0;
77 
78 #ifdef TIMING_OUTPUT
79  using Teuchos::RCP;
80  using Teuchos::Time;
81  matvec_time_ = RCP<Time>(new Time("First 2 Matvec time in prober: "));
82  localize_time_ = RCP<Time>(new Time("Localize time in prober: "));
83  trisolve_time_ = RCP<Time>(new Time("Triangular solve time in prober: "));
84  dist_time_ = RCP<Time>(new Time("Distribute time in prober: "));
85  matvec2_time_ = RCP<Time>(new Time("Second 1 matvec time in prober: "));
86  apply_time_ = RCP<Time>(new Time("Apply time in prober: "));
87  update_time_ = RCP<Time>(new Time("Update time in prober: "));
88 #endif
89 }
90 
91 int ShyLU_Probing_Operator::SetUseTranspose(bool useTranspose)
92 {
93  return 0;
94 }
95 
96 int ShyLU_Probing_Operator::Apply(const Epetra_MultiVector &X,
97  Epetra_MultiVector &Y) const
98 {
99 #ifdef TIMING_OUTPUT
100  apply_time_->start();
101 #endif
102 
103  int nvectors = X.NumVectors();
104  bool local = (C_->Comm().NumProc() == 1);
105  int err;
106  //cout << "No of colors after probing" << nvectors << endl;
107 
108 #ifdef TIMING_OUTPUT
109  matvec_time_->start();
110 #endif
111 
112  err = G_->Multiply(false, X, *temp2);
113  assert(err == 0);
114  if (!local)
115  err = C_->Multiply(false, X, *temp);
116  else
117  {
118  // localize X
119  double *values;
120  int mylda;
121  X.ExtractView(&values, &mylda);
122 
123  Epetra_SerialComm LComm; // Use Serial Comm for the local blocks.
124  Epetra_Map SerialMap(X.Map().NumMyElements(), X.Map().NumMyElements(),
125  X.Map().MyGlobalElements(), 0, LComm);
126  Epetra_MultiVector Xl(View, SerialMap, values, mylda, X.NumVectors());
127  err = C_->Multiply(false, Xl, *temp);
128  }
129  assert(err == 0);
130 
131 #ifdef TIMING_OUTPUT
132  matvec_time_->stop();
133 #endif
134 
135  int nrows = C_->RowMap().NumMyElements();
136 
137 #ifdef DEBUG
138  cout << "DEBUG MODE" << endl;
139  assert(nrows == localDRowMap_->NumGlobalElements());
140 
141  int gids[nrows], gids1[nrows];
142  C_->RowMap().MyGlobalElements(gids);
143  localDRowMap_->MyGlobalElements(gids1);
144 
145  for (int i = 0; i < nrows; i++)
146  {
147  assert(gids[i] == gids1[i]);
148  }
149 #endif
150 
151 #ifdef TIMING_OUTPUT
152  localize_time_->start();
153 #endif
154 
155  //int err;
156  int lda;
157  double *values;
158  if (!local)
159  {
160  err = temp->ExtractView(&values, &lda);
161  assert (err == 0);
162 
163  // copy to local vector //TODO: OMP parallel
164  assert(lda == nrows);
165 
166  //#pragma omp parallel for shared(nvectors, nrows, values)
167  for (int v = 0; v < nvectors; v++)
168  {
169  for (int i = 0; i < nrows; i++)
170  {
171  err = ltemp->ReplaceMyValue(i, v, values[i+v*lda]);
172  assert (err == 0);
173  }
174  }
175  }
176 
177 #ifdef TIMING_OUTPUT
178  localize_time_->stop();
179  trisolve_time_->start();
180 #endif
181 
182  if (!local)
183  {
184  LP_->SetRHS(ltemp.getRawPtr());
185  }
186  else
187  {
188  //LP_->SetRHS(temp.getRawPtr());
189  }
190  //LP_->SetLHS(localX.getRawPtr());
191 
192  //TODO: Why not just in Reset(). Check the distr path.
193  if (config_->amesosForDiagonal)
194  {
195  ssym_->OrigLP->SetLHS(localX.getRawPtr());
196  ssym_->OrigLP->SetRHS(temp.getRawPtr());
197  ssym_->ReIdx_LP->fwd();
198  solver_->Solve();
199  }
200  else
201  ifSolver_->ApplyInverse(*temp, *localX);
202 
203 #ifdef TIMING_OUTPUT
204  trisolve_time_->stop();
205  dist_time_->start();
206 #endif
207 
208  if (!local)
209  {
210  err = localX->ExtractView(&values, &lda);
211  assert (err == 0);
212 
213  //Copy back to dist vector //TODO: OMP parallel
214  //#pragma omp parallel for
215  for (int v = 0; v < nvectors; v++)
216  {
217  for (int i = 0; i < nrows; i++)
218  {
219  err = temp->ReplaceMyValue(i, v, values[i+v*lda]);
220  assert (err == 0);
221  }
222  }
223  }
224 
225 #ifdef TIMING_OUTPUT
226  dist_time_->stop();
227  matvec2_time_->start();
228 #endif
229 
230  if (!local)
231  {
232  R_->Multiply(false, *temp, Y);
233  }
234  else
235  {
236  // Should Y be localY in Multiply and then exported to Y ?? TODO:
237  // Use view mode ?
238  double *values;
239  int mylda;
240  Y.ExtractView(&values, &mylda);
241 
242  Epetra_SerialComm LComm; // Use Serial Comm for the local blocks.
243  Epetra_Map SerialMap(Y.Map().NumMyElements(), Y.Map().NumMyElements(),
244  Y.Map().MyGlobalElements(), 0, LComm);
245  Epetra_MultiVector Yl(View, SerialMap, values, mylda, Y.NumVectors());
246  R_->Multiply(false, *localX, Yl);
247  }
248 
249 #ifdef TIMING_OUTPUT
250  matvec2_time_->stop();
251  update_time_->start();
252 #endif
253 
254  err = Y.Update(1.0, *temp2, -1.0);
255  //cout << Y.MyLength() << " " << temp2.MyLength() << endl;
256  assert(err == 0);
257 
258 #ifdef TIMING_OUTPUT
259  update_time_->stop();
260  apply_time_->stop();
261 #endif
262  cntApply++;
263  return 0;
264 }
265 
266 void ShyLU_Probing_Operator::PrintTimingInfo()
267 {
268 #ifdef TIMING_OUTPUT
269  cout << matvec_time_->name()<< matvec_time_->totalElapsedTime() << endl;
270  cout << localize_time_->name()<< localize_time_->totalElapsedTime() << endl;
271  cout << trisolve_time_->name()<< trisolve_time_->totalElapsedTime() << endl;
272  cout << dist_time_->name() <<dist_time_->totalElapsedTime() << endl;
273  cout << matvec2_time_->name() <<matvec2_time_->totalElapsedTime() << endl;
274  cout << update_time_->name() <<update_time_->totalElapsedTime() << endl;
275  cout << apply_time_->name() <<apply_time_->totalElapsedTime() << endl;
276 #endif
277 }
278 
279 void ShyLU_Probing_Operator::ResetTempVectors(int nvectors)
280 {
281  using Teuchos::RCP;
282  nvectors_ = nvectors;
283  if (nvectors_ == 0) return;
284 
285  if (nvectors <= ssym_->Drhs->NumVectors())
286  {
287  // If vectors were created already, they will be freed.
288  ltemp = Teuchos::RCP<Epetra_MultiVector>
289  (new Epetra_MultiVector(View, *(ssym_->Drhs), 0, nvectors));
290  localX = Teuchos::RCP<Epetra_MultiVector>
291  (new Epetra_MultiVector(View, *(ssym_->Dlhs), 0, nvectors));
292  temp2 = Teuchos::RCP<Epetra_MultiVector>
293  (new Epetra_MultiVector(View, *(ssym_->Gvec), 0, nvectors));
294  }
295  else
296  {
297  ltemp = Teuchos::RCP<Epetra_MultiVector>
298  (new Epetra_MultiVector((*(ssym_->Drhs)).Map(), nvectors));
299  localX = Teuchos::RCP<Epetra_MultiVector>
300  (new Epetra_MultiVector((*(ssym_->Dlhs)).Map(), nvectors));
301  temp2 = Teuchos::RCP<Epetra_MultiVector>
302  (new Epetra_MultiVector((*(ssym_->Gvec)).Map(), nvectors));
303  }
304  temp = RCP<Epetra_MultiVector>(new Epetra_MultiVector(C_->RowMap(),
305  nvectors));
306  ssym_->OrigLP->SetLHS(localX.getRawPtr());
307  ssym_->OrigLP->SetRHS(temp.getRawPtr());
308  ssym_->ReIdx_LP->fwd();
309 }
310 
311 int ShyLU_Probing_Operator::ApplyInverse(const Epetra_MultiVector &X,
312  Epetra_MultiVector &Y) const
313 {
314  return 0;
315 }
316 
317 double ShyLU_Probing_Operator::NormInf() const
318 {
319  return -1;
320 }
321 
322 const char* ShyLU_Probing_Operator::Label() const
323 {
324  return "Shylu probing";
325 }
326 
327 bool ShyLU_Probing_Operator::UseTranspose() const
328 {
329  return false;
330 }
331 
332 bool ShyLU_Probing_Operator::HasNormInf() const
333 {
334  return false;
335 }
336 
337 const Epetra_Comm& ShyLU_Probing_Operator::Comm() const
338 {
339  return G_->Comm();
340 }
341 
342 const Epetra_Map& ShyLU_Probing_Operator::OperatorDomainMap() const
343 {
344  //return C_->ColMap();
345  return G_->DomainMap();
346 }
347 
348 const Epetra_Map& ShyLU_Probing_Operator::OperatorRangeMap() const
349 {
350  //return R_->RowMap();
351  return G_->RangeMap();
352 }