Xpetra_MatrixUtils.hpp
Go to the documentation of this file.
1 // @HEADER
2 //
3 // ***********************************************************************
4 //
5 // Xpetra: A linear algebra interface 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
39 // Jonathan Hu (jhu@sandia.gov)
40 // Andrey Prokopenko (aprokop@sandia.gov)
41 // Ray Tuminaro (rstumin@sandia.gov)
42 // Tobias Wiesner (tawiesn@sandia.gov)
43 //
44 // ***********************************************************************
45 //
46 // @HEADER
47 #ifndef PACKAGES_XPETRA_SUP_MATRIX_UTILS_HPP_
48 #define PACKAGES_XPETRA_SUP_MATRIX_UTILS_HPP_
49 
50 #include "Xpetra_ConfigDefs.hpp"
51 
52 
53 #include "Xpetra_Map.hpp"
54 #include "Xpetra_MapUtils.hpp"
55 #include "Xpetra_StridedMap.hpp"
56 #include "Xpetra_MapFactory.hpp"
57 #include "Xpetra_MapExtractor.hpp"
59 #include "Xpetra_Matrix.hpp"
60 #include "Xpetra_MatrixFactory.hpp"
62 
63 namespace Xpetra {
64 
74 template <class Scalar,
75  class LocalOrdinal,
76  class GlobalOrdinal,
77  class Node>
78 class MatrixUtils {
79 #undef XPETRA_MATRIXUTILS_SHORT
80 #include "Xpetra_UseShortNames.hpp"
81 
82 public:
83 
84  static Teuchos::RCP<Xpetra::MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node> > xpetraGidNumbering2ThyraGidNumbering(
89  RCP<const Map> map = MapUtils::shrinkMapGIDs(*(input.getMap()),*(input.getMap()));
90  RCP<MultiVector> ret = MultiVectorFactory::Build(map, input.getNumVectors(), true);
91  for (size_t c = 0; c < input.getNumVectors(); c++) {
92  Teuchos::ArrayRCP< const Scalar > data = input.getData(c);
93  for (size_t r = 0; r < input.getLocalLength(); r++) {
94  ret->replaceLocalValue(Teuchos::as<LocalOrdinal>(r), c, data[r]);
95  }
96  }
97 
98  return ret;
99  }
100 
101  static Teuchos::RCP<Xpetra::Map<LocalOrdinal, GlobalOrdinal, Node> > findColumnSubMap(
104 
105  RCP< const Teuchos::Comm<int> > comm = input.getRowMap()->getComm();
106 
107  // build an overlapping version of mySpecialMap
108  Teuchos::Array<GlobalOrdinal> ovlUnknownStatusGids;
109  Teuchos::Array<GlobalOrdinal> ovlFoundStatusGids;
110 
111  // loop over global column map of A and find all GIDs where it is not sure, whether they are special or not
112  for(size_t i = 0; i<input.getColMap()->getNodeNumElements(); i++) {
113  GlobalOrdinal gcid = input.getColMap()->getGlobalElement(i);
114  if(domainMap.isNodeGlobalElement(gcid) == false) {
115  ovlUnknownStatusGids.push_back(gcid);
116  }
117  }
118 
119  // We need a locally replicated list of all DOF gids of the (non-overlapping) range map of A10
120  // Communicate the number of DOFs on each processor
121  std::vector<int> myUnknownDofGIDs(comm->getSize(),0);
122  std::vector<int> numUnknownDofGIDs(comm->getSize(),0);
123  myUnknownDofGIDs[comm->getRank()] = ovlUnknownStatusGids.size();
124  Teuchos::reduceAll(*comm,Teuchos::REDUCE_MAX,comm->getSize(),&myUnknownDofGIDs[0],&numUnknownDofGIDs[0]);
125 
126  // create array containing all DOF GIDs
127  size_t cntUnknownDofGIDs = 0;
128  for(int p = 0; p < comm->getSize(); p++) cntUnknownDofGIDs += numUnknownDofGIDs[p];
129  std::vector<GlobalOrdinal> lUnknownDofGIDs(cntUnknownDofGIDs,-1); // local version to be filled
130  std::vector<GlobalOrdinal> gUnknownDofGIDs(cntUnknownDofGIDs,-1); // global version after communication
131  // calculate the offset and fill chunk of memory with local data on each processor
132  size_t cntUnknownOffset = 0;
133  for(int p = 0; p < comm->getRank(); p++) cntUnknownOffset += numUnknownDofGIDs[p];
134  for(size_t k=0; k < Teuchos::as<size_t>(ovlUnknownStatusGids.size()); k++) {
135  lUnknownDofGIDs[k+cntUnknownOffset] = ovlUnknownStatusGids[k];
136  }
137  if(cntUnknownDofGIDs > 0)
138  Teuchos::reduceAll(*comm,Teuchos::REDUCE_MAX,Teuchos::as<int>(cntUnknownDofGIDs),&lUnknownDofGIDs[0],&gUnknownDofGIDs[0]);
139  std::sort(gUnknownDofGIDs.begin(), gUnknownDofGIDs.end());
140  gUnknownDofGIDs.erase(std::unique(gUnknownDofGIDs.begin(), gUnknownDofGIDs.end()), gUnknownDofGIDs.end());
141 
142  // loop through all GIDs with unknown status
143  for(size_t k=0; k < gUnknownDofGIDs.size(); k++) {
144  GlobalOrdinal curgid = gUnknownDofGIDs[k];
145  if(domainMap.isNodeGlobalElement(curgid)) {
146  ovlFoundStatusGids.push_back(curgid); // curgid is in special map (on this processor)
147  }
148  }
149 
150  std::vector<int> myFoundDofGIDs(comm->getSize(),0);
151  std::vector<int> numFoundDofGIDs(comm->getSize(),0);
152  myFoundDofGIDs[comm->getRank()] = ovlFoundStatusGids.size();
153  Teuchos::reduceAll(*comm,Teuchos::REDUCE_MAX,comm->getSize(),&myFoundDofGIDs[0],&numFoundDofGIDs[0]);
154 
155  // create array containing all DOF GIDs
156  size_t cntFoundDofGIDs = 0;
157  for(int p = 0; p < comm->getSize(); p++) cntFoundDofGIDs += numFoundDofGIDs[p];
158  std::vector<GlobalOrdinal> lFoundDofGIDs(cntFoundDofGIDs,-1); // local version to be filled
159  std::vector<GlobalOrdinal> gFoundDofGIDs(cntFoundDofGIDs,-1); // global version after communication
160  // calculate the offset and fill chunk of memory with local data on each processor
161  size_t cntFoundOffset = 0;
162  for(int p = 0; p < comm->getRank(); p++) cntFoundOffset += numFoundDofGIDs[p];
163  for(size_t k=0; k < Teuchos::as<size_t>(ovlFoundStatusGids.size()); k++) {
164  lFoundDofGIDs[k+cntFoundOffset] = ovlFoundStatusGids[k];
165  }
166  if(cntFoundDofGIDs > 0)
167  Teuchos::reduceAll(*comm,Teuchos::REDUCE_MAX,Teuchos::as<int>(cntFoundDofGIDs),&lFoundDofGIDs[0],&gFoundDofGIDs[0]);
168 
169  Teuchos::Array<GlobalOrdinal> ovlDomainMapArray;
170  for(size_t i = 0; i<input.getColMap()->getNodeNumElements(); i++) {
171  GlobalOrdinal gcid = input.getColMap()->getGlobalElement(i);
172  if(domainMap.isNodeGlobalElement(gcid) == true ||
173  std::find(gFoundDofGIDs.begin(), gFoundDofGIDs.end(), gcid) != gFoundDofGIDs.end()) {
174  ovlDomainMapArray.push_back(gcid);
175  }
176  }
177  RCP<Map> ovlDomainMap = MapFactory::Build (domainMap.lib(),Teuchos::OrdinalTraits<GlobalOrdinal>::invalid(),ovlDomainMapArray(),0,comm);
178  return ovlDomainMap;
179  }
180 
189  static Teuchos::RCP<Xpetra::BlockedCrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> > SplitMatrix(
191  Teuchos::RCP<const Xpetra::MapExtractor<Scalar, LocalOrdinal, GlobalOrdinal, Node> > rangeMapExtractor,
192  Teuchos::RCP<const Xpetra::MapExtractor<Scalar, LocalOrdinal, GlobalOrdinal, Node> > domainMapExtractor,
193  Teuchos::RCP<const Xpetra::MapExtractor<Scalar, LocalOrdinal, GlobalOrdinal, Node> > columnMapExtractor = Teuchos::null,
194  bool bThyraMode = false) {
199 
200  size_t numRows = rangeMapExtractor->NumMaps();
201  size_t numCols = domainMapExtractor->NumMaps();
202 
203  TEUCHOS_TEST_FOR_EXCEPTION(rangeMapExtractor->getThyraMode() == true, Xpetra::Exceptions::Incompatible, "Xpetra::MatrixUtils::Split: RangeMapExtractor must not use Thyra style numbering of GIDs. The MapExtractor must contain all GIDs of the full range map in order to define a proper splitting.")
204  TEUCHOS_TEST_FOR_EXCEPTION(domainMapExtractor->getThyraMode() == true, Xpetra::Exceptions::Incompatible, "Xpetra::MatrixUtils::Split: DomainMapExtractor must not use Thyra style numbering of GIDs. The MapExtractor must contain all GIDs of the full domain map in order to define a proper splitting.")
205 
206  RCP<const Map> fullRangeMap = rangeMapExtractor->getFullMap();
207  RCP<const Map> fullDomainMap = domainMapExtractor->getFullMap();
208 
209  TEUCHOS_TEST_FOR_EXCEPTION(fullRangeMap->getMaxAllGlobalIndex() != input.getRowMap()->getMaxAllGlobalIndex(), Xpetra::Exceptions::Incompatible, "Xpetra::MatrixUtils::Split: RangeMapExtractor incompatible to row map of input matrix.")
210  TEUCHOS_TEST_FOR_EXCEPTION(fullRangeMap->getGlobalNumElements() != input.getRowMap()->getGlobalNumElements(), Xpetra::Exceptions::Incompatible, "Xpetra::MatrixUtils::Split: RangeMapExtractor incompatible to row map of input matrix.")
211  TEUCHOS_TEST_FOR_EXCEPTION(fullRangeMap->getNodeNumElements() != input.getRowMap()->getNodeNumElements(), Xpetra::Exceptions::Incompatible, "Xpetra::MatrixUtils::Split: RangeMapExtractor incompatible to row map of input matrix.")
212  TEUCHOS_TEST_FOR_EXCEPTION(fullRangeMap->getMaxAllGlobalIndex() != input.getRangeMap()->getMaxAllGlobalIndex(), Xpetra::Exceptions::Incompatible, "Xpetra::MatrixUtils::Split: RangeMapExtractor incompatible to row map of input matrix.")
213  TEUCHOS_TEST_FOR_EXCEPTION(fullRangeMap->getGlobalNumElements() != input.getRangeMap()->getGlobalNumElements(), Xpetra::Exceptions::Incompatible, "Xpetra::MatrixUtils::Split: RangeMapExtractor incompatible to row map of input matrix.")
214  TEUCHOS_TEST_FOR_EXCEPTION(fullRangeMap->getNodeNumElements() != input.getRangeMap()->getNodeNumElements(), Xpetra::Exceptions::Incompatible, "Xpetra::MatrixUtils::Split: RangeMapExtractor incompatible to row map of input matrix.")
215 
216  TEUCHOS_TEST_FOR_EXCEPTION(fullDomainMap->getMaxAllGlobalIndex() != input.getColMap()->getMaxAllGlobalIndex(), Xpetra::Exceptions::Incompatible, "Xpetra::MatrixUtils::Split: DomainMapExtractor incompatible to domain map of input matrix.")
217  TEUCHOS_TEST_FOR_EXCEPTION(fullDomainMap->getMaxAllGlobalIndex() != input.getDomainMap()->getMaxAllGlobalIndex(), Xpetra::Exceptions::Incompatible, "Xpetra::MatrixUtils::Split: DomainMapExtractor incompatible to domain map of input matrix.")
218  TEUCHOS_TEST_FOR_EXCEPTION(fullDomainMap->getGlobalNumElements() != input.getDomainMap()->getGlobalNumElements(), Xpetra::Exceptions::Incompatible, "Xpetra::MatrixUtils::Split: DomainMapExtractor incompatible to domain map of input matrix.")
219  TEUCHOS_TEST_FOR_EXCEPTION(fullDomainMap->getNodeNumElements() != input.getDomainMap()->getNodeNumElements(), Xpetra::Exceptions::Incompatible, "Xpetra::MatrixUtils::Split: DomainMapExtractor incompatible to domain map of input matrix.")
220 
221  // check column map extractor
222  Teuchos::RCP<const MapExtractor> myColumnMapExtractor = Teuchos::null;
223  if(columnMapExtractor == Teuchos::null) {
224  TEUCHOS_TEST_FOR_EXCEPTION(domainMapExtractor->getThyraMode() == true, Xpetra::Exceptions::Incompatible, "Xpetra::MatrixUtils::Split: Auto generation of column map extractor not supported for Thyra style numbering.");
225  // This code is always executed, since we always provide map extractors in Xpetra numbering!
226  std::vector<Teuchos::RCP<const Map> > ovlxmaps(numCols, Teuchos::null);
227  for(size_t c = 0; c < numCols; c++) {
228  // TODO: is this routine working correctly?
229  Teuchos::RCP<const Map> colMap = MatrixUtils::findColumnSubMap(input, *(domainMapExtractor->getMap(c)));
230  ovlxmaps[c] = colMap;
231  }
232  RCP<const Map> fullColMap = MapUtils::concatenateMaps(ovlxmaps);
233  // This MapExtractor is always in Xpetra mode!
234  myColumnMapExtractor = MapExtractorFactory::Build(fullColMap,ovlxmaps);
235  } else
236  myColumnMapExtractor = columnMapExtractor; // use user-provided column map extractor.
237 
238  // all above MapExtractors are always in Xpetra mode
239  // build separate ones containing Thyra mode GIDs (if necessary)
240  Teuchos::RCP<const MapExtractor> thyRangeMapExtractor = Teuchos::null;
241  Teuchos::RCP<const MapExtractor> thyDomainMapExtractor = Teuchos::null;
242  Teuchos::RCP<const MapExtractor> thyColMapExtractor = Teuchos::null;
243  if(bThyraMode == true) {
244  // build Thyra-style map extractors
245  std::vector<Teuchos::RCP<const Map> > thyRgMapExtractorMaps(numRows, Teuchos::null);
246  for (size_t r = 0; r < numRows; r++) {
247  RCP<const Map> rMap = rangeMapExtractor->getMap(r);
248  RCP<const Map> shrinkedMap = MapUtils::shrinkMapGIDs(*rMap,*rMap);
249  RCP<const StridedMap> strRangeMap = Teuchos::rcp_dynamic_cast<const StridedMap>(rMap);
250  if(strRangeMap != Teuchos::null) {
251  std::vector<size_t> strInfo = strRangeMap->getStridingData();
252  GlobalOrdinal offset = strRangeMap->getOffset();
253  LocalOrdinal stridedBlockId = strRangeMap->getStridedBlockId();
254  RCP<const StridedMap> strShrinkedMap = Teuchos::rcp(new StridedMap(shrinkedMap, strInfo, shrinkedMap->getIndexBase(), stridedBlockId, offset));
255  thyRgMapExtractorMaps[r] = strShrinkedMap;
256  } else {
257  thyRgMapExtractorMaps[r] = shrinkedMap;
258  }
259  TEUCHOS_TEST_FOR_EXCEPTION(thyRgMapExtractorMaps[r]->getNodeNumElements() != rMap->getNodeNumElements(), Xpetra::Exceptions::Incompatible, "Xpetra::MatrixUtils::Split: Thyra-style range map extractor contains faulty data.")
260  }
261  RCP<const Map> fullThyRangeMap = MapUtils::concatenateMaps(thyRgMapExtractorMaps);
262  thyRangeMapExtractor = MapExtractorFactory::Build(fullThyRangeMap,thyRgMapExtractorMaps,true);
263  std::vector<Teuchos::RCP<const Map> > thyDoMapExtractorMaps (numCols, Teuchos::null);
264  std::vector<Teuchos::RCP<const Map> > thyColMapExtractorMaps(numCols, Teuchos::null);
265  for (size_t c = 0; c < numCols; c++) {
266  RCP<const Map> cMap = domainMapExtractor->getMap(c);
267 
268  RCP<const Map> shrinkedDomainMap = MapUtils::shrinkMapGIDs(*cMap,*cMap);
269  RCP<const StridedMap> strDomainMap = Teuchos::rcp_dynamic_cast<const StridedMap>(cMap);
270  if(strDomainMap != Teuchos::null) {
271  std::vector<size_t> strInfo = strDomainMap->getStridingData();
272  GlobalOrdinal offset = strDomainMap->getOffset();
273  LocalOrdinal stridedBlockId = strDomainMap->getStridedBlockId();
274  RCP<const StridedMap> strShrinkedDomainMap = Teuchos::rcp(new StridedMap(shrinkedDomainMap, strInfo, shrinkedDomainMap->getIndexBase(), stridedBlockId, offset));
275  thyDoMapExtractorMaps[c] = strShrinkedDomainMap;
276  } else {
277  thyDoMapExtractorMaps[c] = shrinkedDomainMap;
278  }
279  RCP<const Map> colMap = myColumnMapExtractor->getMap(c);
280  RCP<const Map> shrinkedColMap = MapUtils::shrinkMapGIDs(*colMap,*cMap);
281  RCP<const StridedMap> strColMap = Teuchos::rcp_dynamic_cast<const StridedMap>(colMap);
282  if(strColMap != Teuchos::null) {
283  std::vector<size_t> strInfo = strColMap->getStridingData();
284  GlobalOrdinal offset = strColMap->getOffset();
285  LocalOrdinal stridedBlockId = strColMap->getStridedBlockId();
286  RCP<const StridedMap> strShrinkedColMap = Teuchos::rcp(new StridedMap(shrinkedColMap, strInfo, shrinkedColMap->getIndexBase(), stridedBlockId, offset));
287  thyColMapExtractorMaps[c] = strShrinkedColMap;
288  } else {
289  thyColMapExtractorMaps[c] = shrinkedColMap;
290  }
291 
292  TEUCHOS_TEST_FOR_EXCEPTION(thyColMapExtractorMaps[c]->getNodeNumElements() != colMap->getNodeNumElements(), Xpetra::Exceptions::Incompatible, "Xpetra::MatrixUtils::Split: Thyra-style column map extractor contains faulty data.")
293  TEUCHOS_TEST_FOR_EXCEPTION(thyDoMapExtractorMaps[c]->getNodeNumElements() != cMap->getNodeNumElements(), Xpetra::Exceptions::Incompatible, "Xpetra::MatrixUtils::Split: Thyra-style domain map extractor contains faulty data.")
294  }
295  RCP<const Map> fullThyDomainMap = MapUtils::concatenateMaps(thyDoMapExtractorMaps);
296  RCP<const Map> fullThyColumnMap = MapUtils::concatenateMaps(thyColMapExtractorMaps);
297  thyDomainMapExtractor = MapExtractorFactory::Build(fullThyDomainMap,thyDoMapExtractorMaps,true);
298  thyColMapExtractor = MapExtractorFactory::Build(fullThyColumnMap,thyColMapExtractorMaps,true);
299  }
300  // create submatrices
301  std::vector<Teuchos::RCP<Matrix> > subMatrices(numRows*numCols, Teuchos::null);
302  for (size_t r = 0; r < numRows; r++) {
303  for (size_t c = 0; c < numCols; c++) {
304  // create empty CrsMatrix objects
305  // make sure that the submatrices are defined using the right row maps (either Thyra or xpetra style)
306  // Note: we're reserving a little bit too much memory for the submatrices, but should be still reasonable
307  if(bThyraMode == true)
308  subMatrices[r*numCols+c] = MatrixFactory::Build (thyRangeMapExtractor->getMap(r,true),input.getNodeMaxNumRowEntries());
309  else
310  subMatrices[r*numCols+c] = MatrixFactory::Build (rangeMapExtractor->getMap(r),input.getNodeMaxNumRowEntries());
311  }
312  }
313 
314  // We need a vector which lives on the column map of input and stores the block id that the column belongs to.
315  // create a vector on the domain map. Loop over it and fill in the corresponding block id numbers
316  // create a vector on the column map and import the data
317  // Importer: source map is non-overlapping. Target map is overlapping
318  // call colMap.Import(domMap,Importer,Insert)
319  // do the same with "Add" to make sure only one processor is responsible for the different GIDs!
320 #if 0 // TAW needs to be fixed (does not compile for Scalar=complex)
325 
326  RCP<Vector> doCheck = VectorFactory::Build(input.getDomainMap(), true);
327  doCheck->putScalar(1.0);
328  RCP<Vector> coCheck = VectorFactory::Build(input.getColMap(), true);
329  RCP<Import> imp = ImportFactory::Build(input.getDomainMap(), input.getColMap());
330  coCheck->doImport(*doCheck, *imp, Xpetra::ADD);
331  TEUCHOS_TEST_FOR_EXCEPTION(coCheck->normInf() != Teuchos::ScalarTraits< Scalar >::magnitude(1.0), Xpetra::Exceptions::RuntimeError, "Xpetra::MatrixUtils::SplitMatrix: error when distributing data.");
332 
333  doCheck->putScalar(-1.0);
334  coCheck->putScalar(-1.0);
335 
336  Teuchos::ArrayRCP< Scalar > doCheckData = doCheck->getDataNonConst(0);
337  for (size_t rrr = 0; rrr < input.getDomainMap()->getNodeNumElements(); rrr++) {
338  // global row id to extract data from global monolithic matrix
339  GlobalOrdinal id = input.getDomainMap()->getGlobalElement(rrr); // LID -> GID (column)
340 
341  // Find the block id in range map extractor that belongs to same global id.
342  size_t BlockId = domainMapExtractor->getMapIndexForGID(id);
343 
344  doCheckData[rrr] = Teuchos::as<Scalar>(BlockId);
345  }
346 
347  coCheck->doImport(*doCheck, *imp, Xpetra::INSERT);
348 
349  Teuchos::ArrayRCP< Scalar > coCheckData = coCheck->getDataNonConst(0);
350 #endif
351  // loop over all rows of input matrix
352  for (size_t rr = 0; rr < input.getRowMap()->getNodeNumElements(); rr++) {
353 
354  // global row id to extract data from global monolithic matrix
355  GlobalOrdinal growid = input.getRowMap()->getGlobalElement(rr); // LID -> GID (column)
356 
357  // Find the block id in range map extractor that belongs to same global row id.
358  // We assume the global ids to be unique in the map extractor
359  // The MapExtractor objects may be constructed using the thyra mode. However, we use
360  // global GID ids (as we have a single input matrix). The only problematic thing could
361  // be that the GIDs of the input matrix are inconsistent with the GIDs in the map extractor.
362  // Note, that for the Thyra mode the GIDs in the map extractor are generated using the ordering
363  // of the blocks.
364  size_t rowBlockId = rangeMapExtractor->getMapIndexForGID(growid);
365 
366  // global row id used for subblocks to insert information
367  GlobalOrdinal subblock_growid = growid; // for Xpetra-style numbering the global row ids are not changing
368  if(bThyraMode == true) {
369  // find local row id associated with growid in the corresponding subblock
370  LocalOrdinal lrowid = rangeMapExtractor->getMap(rowBlockId)->getLocalElement(growid);
371  // translate back local row id to global row id for the subblock
372  subblock_growid = thyRangeMapExtractor->getMap(rowBlockId,true)->getGlobalElement(lrowid);
373  }
374 
375  // extract matrix entries from input matrix
376  // we use global ids since we have to determine the corresponding
377  // block column ids using the global ids anyway
378  Teuchos::ArrayView<const LocalOrdinal> indices;
379  Teuchos::ArrayView<const Scalar> vals;
380  input.getLocalRowView(rr, indices, vals);
381 
382  std::vector<Teuchos::Array<GlobalOrdinal> > blockColIdx (numCols, Teuchos::Array<GlobalOrdinal>());
383  std::vector<Teuchos::Array<Scalar> > blockColVals(numCols, Teuchos::Array<Scalar>());
384 
385  for(size_t i=0; i<(size_t)indices.size(); i++) {
386  // gobal column id to extract data from full monolithic matrix
387  GlobalOrdinal gcolid = input.getColMap()->getGlobalElement(indices[i]);
388 
389  size_t colBlockId = myColumnMapExtractor->getMapIndexForGID(gcolid); // old buggy thing
390  //size_t colBlockId = Teuchos::as<size_t>(coCheckData[indices[i]]);
391 
392  // global column id used for subblocks to insert information
393  GlobalOrdinal subblock_gcolid = gcolid; // for Xpetra-style numbering the global col ids are not changing
394  if(bThyraMode == true) {
395  // find local col id associated with gcolid in the corresponding subblock
396  LocalOrdinal lcolid = myColumnMapExtractor->getMap(colBlockId)->getLocalElement(gcolid);
397  // translate back local col id to global col id for the subblock
398  subblock_gcolid = thyColMapExtractor->getMap(colBlockId,true)->getGlobalElement(lcolid);
399  }
400  blockColIdx [colBlockId].push_back(subblock_gcolid);
401  blockColVals[colBlockId].push_back(vals[i]);
402  }
403 
404  for (size_t c = 0; c < numCols; c++) {
405  subMatrices[rowBlockId*numCols+c]->insertGlobalValues(subblock_growid,blockColIdx[c].view(0,blockColIdx[c].size()),blockColVals[c].view(0,blockColVals[c].size()));
406  }
407 
408  }
409 
410  // call fill complete on subblocks and create BlockedCrsOperator
411  RCP<BlockedCrsMatrix> bA = Teuchos::null;
412  if(bThyraMode == true) {
413  for (size_t r = 0; r < numRows; r++) {
414  for (size_t c = 0; c < numCols; c++) {
415  subMatrices[r*numCols+c]->fillComplete(thyDomainMapExtractor->getMap(c,true), thyRangeMapExtractor->getMap(r,true));
416  }
417  }
418  bA = Teuchos::rcp(new BlockedCrsMatrix(thyRangeMapExtractor, thyDomainMapExtractor, 10 /*input.getRowMap()->getNodeMaxNumRowEntries()*/));
419  } else {
420  for (size_t r = 0; r < numRows; r++) {
421  for (size_t c = 0; c < numCols; c++) {
422  subMatrices[r*numCols+c]->fillComplete(domainMapExtractor->getMap(c), rangeMapExtractor->getMap(r));
423  }
424  }
425  bA = Teuchos::rcp(new BlockedCrsMatrix(rangeMapExtractor, domainMapExtractor, 10 /*input.getRowMap()->getNodeMaxNumRowEntries()*/));
426  }
427 
428  for (size_t r = 0; r < numRows; r++) {
429  for (size_t c = 0; c < numCols; c++) {
430  bA->setMatrix(r,c,subMatrices[r*numCols+c]);
431  }
432  }
433  return bA;
434  }
435 };
436 
437 } // end namespace Xpetra
438 
439 #define XPETRA_MATRIXUTILS_SHORT
440 
441 #endif // PACKAGES_XPETRA_SUP_MATRIX_UTILS_HPP_
virtual size_t getNodeMaxNumRowEntries() const =0
Returns the maximum number of entries across all rows/columns on this node.
static Teuchos::RCP< Xpetra::Map< LocalOrdinal, GlobalOrdinal, Node > > findColumnSubMap(const Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > &input, const Xpetra::Map< LocalOrdinal, GlobalOrdinal, Node > &domainMap)
static RCP< Import< LocalOrdinal, GlobalOrdinal, Node > > Build(const RCP< const Map< LocalOrdinal, GlobalOrdinal, Node > > &source, const RCP< const Map< LocalOrdinal, GlobalOrdinal, Node > > &target)
Constructor specifying the number of non-zeros for all rows.
std::vector< size_t > getStridingData() const
static Teuchos::RCP< Map< LocalOrdinal, GlobalOrdinal, Node > > Build(UnderlyingLib lib, global_size_t numGlobalElements, GlobalOrdinal indexBase, const Teuchos::RCP< const Teuchos::Comm< int > > &comm, LocalGlobal lg=Xpetra::GloballyDistributed, const Teuchos::RCP< Node > &node=defaultArgNode())
Map constructor with Xpetra-defined contiguous uniform distribution.
virtual const RCP< const Map > & getRowMap() const
Returns the Map that describes the row distribution in this matrix.
Xpetra namespace
Exception throws to report errors in the internal logical of the program.
Xpetra utility class for common map-related routines.
static Teuchos::RCP< const Xpetra::Map< LocalOrdinal, GlobalOrdinal, Node > > concatenateMaps(const std::vector< Teuchos::RCP< const Xpetra::Map< LocalOrdinal, GlobalOrdinal, Node > > > &subMaps)
Helper function to concatenate several maps.
static Teuchos::RCP< Xpetra::BlockedCrsMatrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > > SplitMatrix(const Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > &input, Teuchos::RCP< const Xpetra::MapExtractor< Scalar, LocalOrdinal, GlobalOrdinal, Node > > rangeMapExtractor, Teuchos::RCP< const Xpetra::MapExtractor< Scalar, LocalOrdinal, GlobalOrdinal, Node > > domainMapExtractor, Teuchos::RCP< const Xpetra::MapExtractor< Scalar, LocalOrdinal, GlobalOrdinal, Node > > columnMapExtractor=Teuchos::null, bool bThyraMode=false)
virtual const RCP< const Map > & getColMap() const
Returns the Map that describes the column distribution in this matrix. This might be null until fillC...
virtual bool isNodeGlobalElement(GlobalOrdinal globalIndex) const =0
Whether the given global index is valid for this Map on this process.
virtual void getLocalRowView(LocalOrdinal LocalRow, ArrayView< const LocalOrdinal > &indices, ArrayView< const Scalar > &values) const =0
Extract a const, non-persisting view of local indices in a specified row of the matrix.
virtual Teuchos::RCP< const Map< LocalOrdinal, GlobalOrdinal, Node > > getMap() const=0
The Map describing the parallel distribution of this object.
virtual size_t getLocalLength() const =0
Local number of rows on the calling process.
static Teuchos::RCP< Xpetra::Map< LocalOrdinal, GlobalOrdinal, Node > > shrinkMapGIDs(const Xpetra::Map< LocalOrdinal, GlobalOrdinal, Node > &input, const Xpetra::Map< LocalOrdinal, GlobalOrdinal, Node > &nonOvlInput)
Helper function to shrink the GIDs and generate a standard map whith GIDs starting at 0...
static RCP< Vector > Build(const Teuchos::RCP< const Map > &map, bool zeroOut=true)
Constructor specifying the number of non-zeros for all rows.
static Teuchos::RCP< Xpetra::MultiVector< Scalar, LocalOrdinal, GlobalOrdinal, Node > > xpetraGidNumbering2ThyraGidNumbering(const Xpetra::MultiVector< Scalar, LocalOrdinal, GlobalOrdinal, Node > &input)
static Teuchos::RCP< Xpetra::MapExtractor< Scalar, LocalOrdinal, GlobalOrdinal, Node > > Build(const Teuchos::RCP< const Map > &fullmap, const std::vector< Teuchos::RCP< const Map > > &maps, bool bThyraMode=false)
Constructor specifying the Maps.
Exception throws to report incompatible objects (like maps).
virtual Teuchos::RCP< const Map > getRangeMap() const =0
The Map associated with the range of this operator, which must be compatible with Y...
static Teuchos::RCP< MultiVector< Scalar, LocalOrdinal, GlobalOrdinal, Node > > Build(const Teuchos::RCP< const Map< LocalOrdinal, GlobalOrdinal, Node > > &map, size_t NumVectors, bool zeroOut=true)
Constructor specifying the number of non-zeros for all rows.
static RCP< Matrix > Build(const RCP< const Map > &rowMap, size_t maxNumEntriesPerRow, Xpetra::ProfileType pftype=Xpetra::DynamicProfile)
Constructor specifying the number of non-zeros for all rows.
virtual size_t getNumVectors() const =0
Number of columns in the multivector.
virtual UnderlyingLib lib() const =0
Get the library used by this object (Tpetra or Epetra?)
virtual Teuchos::ArrayRCP< const Scalar > getData(size_t j) const =0
Const view of the local values in a particular vector of this multivector.
Xpetra utility class for common matrix-related routines.
virtual Teuchos::RCP< const Map > getDomainMap() const =0
The Map associated with the domain of this operator, which must be compatible with X...
Xpetra-specific matrix class.
Class that stores a strided map.