Amesos2 - Direct Sparse Solver Interfaces  Version of the Day
Amesos2_TpetraMultiVecAdapter_def.hpp
Go to the documentation of this file.
1 // @HEADER
2 //
3 // ***********************************************************************
4 //
5 // Amesos2: Templated Direct Sparse Solver Package
6 // Copyright 2011 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 //
42 // @HEADER
43 
53 #ifndef AMESOS2_TPETRA_MULTIVEC_ADAPTER_DEF_HPP
54 #define AMESOS2_TPETRA_MULTIVEC_ADAPTER_DEF_HPP
55 
57 
58 
59 namespace Amesos2 {
60 
61  using Tpetra::MultiVector;
62 
63  template <typename Scalar, typename LocalOrdinal, typename GlobalOrdinal, class Node >
64  MultiVecAdapter<
65  MultiVector<Scalar,
66  LocalOrdinal,
67  GlobalOrdinal,
68  Node> >::MultiVecAdapter( const Teuchos::RCP<multivec_t>& m )
69  : mv_(m)
70  {}
71 
72 
73  template <typename Scalar, typename LocalOrdinal, typename GlobalOrdinal, class Node >
74  void
75  MultiVecAdapter<
76  MultiVector<Scalar,
77  LocalOrdinal,
78  GlobalOrdinal,
79  Node> >::get1dCopy(const Teuchos::ArrayView<scalar_t>& av,
80  size_t lda,
81  Teuchos::Ptr<
82  const Tpetra::Map<LocalOrdinal,
83  GlobalOrdinal,
84  Node> > distribution_map ) const
85  {
86  using Teuchos::as;
87  using Teuchos::RCP;
88  typedef Tpetra::Map<LocalOrdinal, GlobalOrdinal, Node> map_type;
89  const size_t num_vecs = getGlobalNumVectors ();
90 
91  TEUCHOS_TEST_FOR_EXCEPTION(
92  distribution_map.getRawPtr () == NULL, std::invalid_argument,
93  "Amesos2::MultiVecAdapter::get1dCopy: distribution_map argument is null.");
94  TEUCHOS_TEST_FOR_EXCEPTION(
95  mv_.is_null (), std::logic_error,
96  "Amesos2::MultiVecAdapter::get1dCopy: mv_ is null.");
97  // Check mv_ before getMap(), because the latter calls mv_->getMap().
98  TEUCHOS_TEST_FOR_EXCEPTION(
99  this->getMap ().is_null (), std::logic_error,
100  "Amesos2::MultiVecAdapter::get1dCopy: this->getMap() returns null.");
101 
102 #ifdef HAVE_AMESOS2_DEBUG
103  const size_t requested_vector_length = distribution_map->getNodeNumElements ();
104  TEUCHOS_TEST_FOR_EXCEPTION(
105  lda < requested_vector_length, std::invalid_argument,
106  "Amesos2::MultiVecAdapter::get1dCopy: On process " <<
107  distribution_map->getComm ()->getRank () << " of the distribution Map's "
108  "communicator, the given stride lda = " << lda << " is not large enough "
109  "for the local vector length " << requested_vector_length << ".");
110  TEUCHOS_TEST_FOR_EXCEPTION(
111  as<size_t> (av.size ()) < as<size_t> ((num_vecs - 1) * lda + requested_vector_length),
112  std::invalid_argument, "Amesos2::MultiVector::get1dCopy: MultiVector "
113  "storage not large enough given leading dimension and number of vectors." );
114 #endif // HAVE_AMESOS2_DEBUG
115 
116  // (Re)compute the Export object if necessary. If not, then we
117  // don't need to clone distribution_map; we can instead just get
118  // the previously cloned target Map from the Export object.
119  RCP<const map_type> distMap;
120  if (exporter_.is_null () ||
121  ! exporter_->getSourceMap ()->isSameAs (* (this->getMap ())) ||
122  ! exporter_->getTargetMap ()->isSameAs (* distribution_map)) {
123 
124  // Since we're caching the Export object, and since the Export
125  // needs to keep the distribution Map, we have to make a copy of
126  // the latter in order to ensure that it will stick around past
127  // the scope of this function call. (Ptr is not reference
128  // counted.) Map's clone() method suffices, even though it only
129  // makes a shallow copy of some of Map's data, because Map is
130  // immutable and those data are reference-counted (e.g.,
131  // ArrayRCP or RCP).
132  distMap = distribution_map->template clone<Node> (distribution_map->getNode ());
133 
134  // (Re)create the Export object.
135  exporter_ = rcp (new export_type (this->getMap (), distMap));
136  }
137  else {
138  distMap = exporter_->getTargetMap ();
139  }
140 
141  multivec_t redist_mv (distMap, num_vecs);
142 
143  // Redistribute the input (multi)vector.
144  redist_mv.doExport (*mv_, *exporter_, Tpetra::REPLACE);
145 
146  // Copy the imported (multi)vector's data into the ArrayView.
147  redist_mv.get1dCopy (av, lda);
148  }
149 
150  template <typename Scalar, typename LocalOrdinal, typename GlobalOrdinal, class Node >
151  Teuchos::ArrayRCP<Scalar>
152  MultiVecAdapter<
153  MultiVector<Scalar,
154  LocalOrdinal,
155  GlobalOrdinal,
156  Node> >::get1dViewNonConst (bool local)
157  {
158  // FIXME (mfh 22 Jan 2014) When I first found this routine, all of
159  // its code was commented out, and it didn't return anything. The
160  // latter is ESPECIALLY dangerous, given that the return value is
161  // an ArrayRCP. Thus, I added the exception throw below.
162  TEUCHOS_TEST_FOR_EXCEPTION(
163  true, std::logic_error, "Amesos2::MultiVecAdapter::get1dViewNonConst: "
164  "Not implemented.");
165 
166  // if( local ){
167  // this->localize();
168  // /* Use the global element list returned by
169  // * mv_->getMap()->getNodeElementList() to get a subCopy of mv_ which we
170  // * assign to l_l_mv_, then return get1dViewNonConst() of l_l_mv_
171  // */
172  // if(l_l_mv_.is_null() ){
173  // Teuchos::ArrayView<const GlobalOrdinal> nodeElements_go
174  // = mv_->getMap()->getNodeElementList();
175  // Teuchos::Array<size_t> nodeElements_st(nodeElements_go.size());
176 
177  // // Convert the node element to a list of size_t type objects
178  // typename Teuchos::ArrayView<const GlobalOrdinal>::iterator it_go;
179  // Teuchos::Array<size_t>::iterator it_st = nodeElements_st.begin();
180  // for( it_go = nodeElements_go.begin(); it_go != nodeElements_go.end(); ++it_go ){
181  // *(it_st++) = Teuchos::as<size_t>(*it_go);
182  // }
183 
184  // // To be consistent with the globalize() function, get a view of the local mv
185  // l_l_mv_ = l_mv_->subViewNonConst(nodeElements_st);
186 
187  // return(l_l_mv_->get1dViewNonConst());
188  // } else {
189  // // We need to re-import values to the local, since changes may have been
190  // // made to the global structure that are not reflected in the local
191  // // view.
192  // Teuchos::ArrayView<const GlobalOrdinal> nodeElements_go
193  // = mv_->getMap()->getNodeElementList();
194  // Teuchos::Array<size_t> nodeElements_st(nodeElements_go.size());
195 
196  // // Convert the node element to a list of size_t type objects
197  // typename Teuchos::ArrayView<const GlobalOrdinal>::iterator it_go;
198  // Teuchos::Array<size_t>::iterator it_st = nodeElements_st.begin();
199  // for( it_go = nodeElements_go.begin(); it_go != nodeElements_go.end(); ++it_go ){
200  // *(it_st++) = Teuchos::as<size_t>(*it_go);
201  // }
202 
203  // return l_l_mv_->get1dViewNonConst();
204  // }
205  // } else {
206  // if( mv_->isDistributed() ){
207  // this->localize();
208 
209  // return l_mv_->get1dViewNonConst();
210  // } else { // not distributed, no need to import
211  // return mv_->get1dViewNonConst();
212  // }
213  // }
214  }
215 
216 
217  template <typename Scalar, typename LocalOrdinal, typename GlobalOrdinal, class Node>
218  void
219  MultiVecAdapter<
220  MultiVector<Scalar,
221  LocalOrdinal,
222  GlobalOrdinal,
223  Node> >::put1dData(const Teuchos::ArrayView<const scalar_t>& new_data,
224  size_t lda,
225  Teuchos::Ptr<
226  const Tpetra::Map<LocalOrdinal,
227  GlobalOrdinal,
228  Node> > source_map)
229  {
230  using Teuchos::RCP;
231  typedef Tpetra::Map<LocalOrdinal, GlobalOrdinal, Node> map_type;
232 
233  TEUCHOS_TEST_FOR_EXCEPTION(
234  source_map.getRawPtr () == NULL, std::invalid_argument,
235  "Amesos2::MultiVecAdapter::put1dData: source_map argument is null.");
236  TEUCHOS_TEST_FOR_EXCEPTION(
237  mv_.is_null (), std::logic_error,
238  "Amesos2::MultiVecAdapter::put1dData: the internal MultiVector mv_ is null.");
239  // getMap() calls mv_->getMap(), so test first whether mv_ is null.
240  TEUCHOS_TEST_FOR_EXCEPTION(
241  this->getMap ().is_null (), std::logic_error,
242  "Amesos2::MultiVecAdapter::put1dData: this->getMap() returns null.");
243 
244  const size_t num_vecs = getGlobalNumVectors ();
245 
246  // (Re)compute the Import object if necessary. If not, then we
247  // don't need to clone source_map; we can instead just get the
248  // previously cloned source Map from the Import object.
249  RCP<const map_type> srcMap;
250  if (importer_.is_null () ||
251  ! importer_->getSourceMap ()->isSameAs (* source_map) ||
252  ! importer_->getTargetMap ()->isSameAs (* (this->getMap ()))) {
253 
254  // Since we're caching the Import object, and since the Import
255  // needs to keep the source Map, we have to make a copy of the
256  // latter in order to ensure that it will stick around past the
257  // scope of this function call. (Ptr is not reference counted.)
258  // Map's clone() method suffices, even though it only makes a
259  // shallow copy of some of Map's data, because Map is immutable
260  // and those data are reference-counted (e.g., ArrayRCP or RCP).
261  srcMap = source_map->template clone<Node> (source_map->getNode ());
262  importer_ = rcp (new import_type (srcMap, this->getMap ()));
263  }
264  else {
265  srcMap = importer_->getSourceMap ();
266  }
267 
268  const multivec_t source_mv (srcMap, new_data, lda, num_vecs);
269 
270  // Redistribute the output (multi)vector.
271  mv_->doImport (source_mv, *importer_, Tpetra::REPLACE);
272  }
273 
274 
275  template <typename Scalar, typename LocalOrdinal, typename GlobalOrdinal, class Node >
276  std::string
277  MultiVecAdapter<
278  MultiVector<Scalar,
279  LocalOrdinal,
280  GlobalOrdinal,
281  Node> >::description() const
282  {
283  std::ostringstream oss;
284  oss << "Amesos2 adapter wrapping: ";
285  oss << mv_->description();
286  return oss.str();
287  }
288 
289 
290  template <typename Scalar, typename LocalOrdinal, typename GlobalOrdinal, class Node >
291  void
292  MultiVecAdapter<
293  MultiVector<Scalar,
294  LocalOrdinal,
295  GlobalOrdinal,
296  Node> >::describe (Teuchos::FancyOStream& os,
297  const Teuchos::EVerbosityLevel verbLevel) const
298  {
299  mv_->describe (os, verbLevel);
300  }
301 
302 
303  template <typename Scalar, typename LocalOrdinal, typename GlobalOrdinal, class Node >
304  const char* MultiVecAdapter<
305  MultiVector<Scalar,
306  LocalOrdinal,
307  GlobalOrdinal,
308  Node> >::name = "Amesos2 adapter for Tpetra::MultiVector";
309 
310 } // end namespace Amesos2
311 
312 #endif // AMESOS2_TPETRA_MULTIVEC_ADAPTER_DEF_HPP
Definition: Amesos2_AbstractConcreteMatrixAdapter.hpp:48
Amesos2::MultiVecAdapter specialization for the Tpetra::MultiVector class.