Tpetra parallel linear algebra  Version of the Day
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
Tpetra_applyDirichletBoundaryCondition.hpp
Go to the documentation of this file.
1 /*
2 // @HEADER
3 // ***********************************************************************
4 //
5 // Tpetra: Templated Linear Algebra Services Package
6 // Copyright (2008) 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 #ifndef TPETRA_APPLYDIRICHLETBOUNDARYCONDITION_HPP
45 #define TPETRA_APPLYDIRICHLETBOUNDARYCONDITION_HPP
46 
50 
51 #include "Tpetra_CrsMatrix.hpp"
52 #include "Tpetra_Vector.hpp"
53 #include "Tpetra_Map.hpp"
54 #include "KokkosSparse_CrsMatrix.hpp"
55 #include "Kokkos_ArithTraits.hpp"
56 
57 namespace Tpetra {
58 
68 template<class CrsMatrixType>
69 void
70 applyDirichletBoundaryConditionToLocalMatrixRows
71 (const typename CrsMatrixType::execution_space& execSpace,
72  CrsMatrixType& A,
73  const Kokkos::View<
74  typename CrsMatrixType::local_ordinal_type*,
75  typename CrsMatrixType::device_type> & lclRowInds);
76 
86 template<class CrsMatrixType>
87 void
88 applyDirichletBoundaryConditionToLocalMatrixRows
89 (CrsMatrixType& A,
90  const Kokkos::View<
91  typename CrsMatrixType::local_ordinal_type*,
92  typename CrsMatrixType::device_type>& lclRowInds);
93 
103 template<class CrsMatrixType>
104 void
105 applyDirichletBoundaryConditionToLocalMatrixRows
106 (CrsMatrixType& A,
107  const Kokkos::View<
108  typename CrsMatrixType::local_ordinal_type*,
109  Kokkos::HostSpace> & lclRowInds);
110 
111 
121 template<class CrsMatrixType>
122 void
123 applyDirichletBoundaryConditionToLocalMatrixRowsAndColumns
124 (const typename CrsMatrixType::execution_space& execSpace,
125  CrsMatrixType& A,
126  const Kokkos::View<
127  typename CrsMatrixType::local_ordinal_type*,
128  typename CrsMatrixType::device_type> & lclRowInds);
129 
130 
131 
132 
141 
143 template<class CrsMatrixType>
144 void
145 applyDirichletBoundaryConditionToLocalMatrixRowsAndColumns
146 (CrsMatrixType& A,
147  const Kokkos::View<
148  typename CrsMatrixType::local_ordinal_type*,
149  Kokkos::HostSpace> & lclRowInds);
150 
160 template<class CrsMatrixType>
161 void
162 applyDirichletBoundaryConditionToLocalMatrixRowsAndColumns
163 (CrsMatrixType& A,
164  const Kokkos::View<
165  typename CrsMatrixType::local_ordinal_type*,
166  typename CrsMatrixType::device_type> & lclRowInds);
167 
168 
169 namespace Details {
170 
171 template<class SC, class LO, class GO, class NT>
172 struct ApplyDirichletBoundaryConditionToLocalMatrixRows {
173  using crs_matrix_type = ::Tpetra::CrsMatrix<SC, LO, GO, NT>;
174  using execution_space = typename crs_matrix_type::execution_space;
175  using local_row_indices_type =
176  Kokkos::View<const LO*, Kokkos::AnonymousSpace>;
177 
178  static void
179  run (const execution_space& execSpace,
180  crs_matrix_type& A,
181  const local_row_indices_type& lclRowInds,
182  const bool runOnHost)
183  {
184  // Notes for future refactoring: This routine seems to have one more layer
185  // of options than it probably needs. For instance, if you passed a Kokkos::Serial
186  // execution_space instance as the first argument you probably wound't need the runOnHost
187  // option and then the code below could be collapsed out removing one of the parallel_for's
188 
189  using IST = typename crs_matrix_type::impl_scalar_type;
190  using KAT = Kokkos::ArithTraits<IST>;
191 
192  const auto rowMap = A.getRowMap ();
193  TEUCHOS_TEST_FOR_EXCEPTION
194  (rowMap.get () == nullptr, std::invalid_argument,
195  "The matrix must have a row Map.");
196  const auto colMap = A.getColMap ();
197  TEUCHOS_TEST_FOR_EXCEPTION
198  (colMap.get () == nullptr, std::invalid_argument,
199  "The matrix must have a column Map.");
200  auto A_lcl = A.getLocalMatrixDevice ();
201 
202  const LO lclNumRows = static_cast<LO> (rowMap->getLocalNumElements ());
203  TEUCHOS_TEST_FOR_EXCEPTION
204  (lclNumRows != 0 && static_cast<LO>(A_lcl.graph.numRows ()) != lclNumRows,
205  std::invalid_argument, "The matrix must have been either created "
206  "with a KokkosSparse::CrsMatrix, or must have been fill-completed "
207  "at least once.");
208 
209  auto lclRowMap = A.getRowMap ()->getLocalMap ();
210  auto lclColMap = A.getColMap ()->getLocalMap ();
211  auto rowptr = A_lcl.graph.row_map;
212  auto colind = A_lcl.graph.entries;
213  auto values = A_lcl.values;
214 
215  const bool wasFillComplete = A.isFillComplete ();
216  if (wasFillComplete) {
217  A.resumeFill ();
218  }
219 
220  const LO numInputRows = lclRowInds.extent (0);
221  if (! runOnHost) {
222  using range_type = Kokkos::RangePolicy<execution_space, LO>;
223  Kokkos::parallel_for
224  ("Tpetra::CrsMatrix apply Dirichlet: Device",
225  range_type (execSpace, 0, numInputRows),
226  KOKKOS_LAMBDA (const LO i) {
227  LO row = lclRowInds(i);
228  const GO row_gid = lclRowMap.getGlobalElement(row);
229  for (auto j = rowptr(row); j < rowptr(row+1); ++j) {
230  const bool diagEnt =
231  lclColMap.getGlobalElement (colind(j)) == row_gid;
232  values(j) = diagEnt ? KAT::one () : KAT::zero ();
233  }
234  });
235  }
236  else {
237  using range_type =
238  Kokkos::RangePolicy<Kokkos::DefaultHostExecutionSpace, LO>;
239  Kokkos::parallel_for
240  ("Tpetra::CrsMatrix apply Dirichlet: Host",
241  range_type (0, numInputRows),
242  [&] (const LO i) {
243  LO row = lclRowInds(i);
244  const GO row_gid = lclRowMap.getGlobalElement(row);
245  for (auto j = rowptr(row); j < rowptr(row+1); ++j) {
246  const bool diagEnt =
247  lclColMap.getGlobalElement (colind(j)) == row_gid;
248  values(j) = diagEnt ? KAT::one () : KAT::zero ();
249  }
250  });
251  }
252  if (wasFillComplete) {
253  A.fillComplete (A.getDomainMap (), A.getRangeMap ());
254  }
255  }
256 };
257 
258 
259 
260 template<class SC, class LO, class GO, class NT>
261 struct ApplyDirichletBoundaryConditionToLocalMatrixColumns {
262  using crs_matrix_type = ::Tpetra::CrsMatrix<SC, LO, GO, NT>;
263  using execution_space = typename crs_matrix_type::execution_space;
264  using local_col_flag_type =
265  Kokkos::View<bool*, Kokkos::AnonymousSpace>;
266 
267 
268  static void
269  run (const execution_space& execSpace,
270  crs_matrix_type& A,
271  const local_col_flag_type& lclColFlags,
272  const bool runOnHost)
273  {
274  // Notes for future refactoring: This routine seems to have one more layer
275  // of options than it probably needs. For instance, if you passed a Kokkos::Serial
276  // execution_space instance as the first argument you probably wound't need the runOnHost
277  // option and then the code below could be collapsed out removing one of the parallel_for's
278 
279  using IST = typename crs_matrix_type::impl_scalar_type;
280  using KAT = Kokkos::ArithTraits<IST>;
281 
282  const auto rowMap = A.getRowMap ();
283  TEUCHOS_TEST_FOR_EXCEPTION
284  (rowMap.get () == nullptr, std::invalid_argument,
285  "The matrix must have a row Map.");
286  const auto colMap = A.getColMap ();
287  TEUCHOS_TEST_FOR_EXCEPTION
288  (colMap.get () == nullptr, std::invalid_argument,
289  "The matrix must have a column Map.");
290  auto A_lcl = A.getLocalMatrixDevice ();
291 
292  const LO lclNumRows = static_cast<LO> (rowMap->getLocalNumElements ());
293  TEUCHOS_TEST_FOR_EXCEPTION
294  (lclNumRows != 0 && static_cast<LO>(A_lcl.graph.numRows ()) != lclNumRows,
295  std::invalid_argument, "The matrix must have been either created "
296  "with a KokkosSparse::CrsMatrix, or must have been fill-completed "
297  "at least once.");
298 
299  auto lclRowMap = A.getRowMap()->getLocalMap ();
300  auto lclColMap = A.getColMap()->getLocalMap ();
301  auto rowptr = A_lcl.graph.row_map;
302  auto colind = A_lcl.graph.entries;
303  auto values = A_lcl.values;
304 
305  const bool wasFillComplete = A.isFillComplete ();
306  if (wasFillComplete) {
307  A.resumeFill ();
308  }
309 
310  const LO numRows = (LO) A.getRowMap()->getLocalNumElements();
311  if (! runOnHost) {
312  using range_type = Kokkos::RangePolicy<execution_space, LO>;
313  Kokkos::parallel_for
314  ("Tpetra::CrsMatrix apply Dirichlet cols: Device",
315  range_type (execSpace, 0, numRows),
316  KOKKOS_LAMBDA (const LO i) {
317  for (auto j = rowptr(i); j < rowptr(i+1); ++j) {
318  if(lclColFlags[colind[j]])
319  values(j) = KAT::zero();
320  }
321  });
322  }
323  else {
324  using range_type =
325  Kokkos::RangePolicy<Kokkos::DefaultHostExecutionSpace, LO>;
326  Kokkos::parallel_for
327  ("Tpetra::CrsMatrix apply Dirichlet cols: Host",
328  range_type (0, numRows),
329  KOKKOS_LAMBDA (const LO i) {
330  for (auto j = rowptr(i); j < rowptr(i+1); ++j) {
331  if(lclColFlags[colind[j]])
332  values(j) = KAT::zero();
333  }
334  });
335  }
336  if (wasFillComplete) {
337  A.fillComplete (A.getDomainMap (), A.getRangeMap ());
338  }
339  }
340 };
341 
342 
343 
344 template<class SC, class LO, class GO, class NT>
345 void localRowsToColumns(const typename ::Tpetra::CrsMatrix<SC, LO, GO, NT>::execution_space& execSpace, const ::Tpetra::CrsMatrix<SC, LO, GO, NT>& A,const Kokkos::View<const LO*,Kokkos::AnonymousSpace> & dirichletRowIds, Kokkos::View<bool*,Kokkos::AnonymousSpace> & dirichletColFlags) {
346  using crs_matrix_type = ::Tpetra::CrsMatrix<SC, LO, GO, NT>;
347  using execution_space = typename crs_matrix_type::execution_space;
348  using memory_space = typename crs_matrix_type::device_type::memory_space;
349 
350  // Need a colMap
351  TEUCHOS_TEST_FOR_EXCEPTION(A.getColMap().get () == nullptr, std::invalid_argument,"The matrix must have a column Map.");
352 
353  // NOTE: We assume that the RowMap and DomainMap of the matrix match.
354  // This could get relaxed at a later date, if we need that functionality
355  TEUCHOS_TEST_FOR_EXCEPTION(!A.getRowMap()->isSameAs(*A.getDomainMap()),std::invalid_argument, "localRowsToColumns: Row/Domain maps do not match");
356 
357  // Assume that the dirichletColFlags array is the correct size
358  TEUCHOS_TEST_FOR_EXCEPTION(A.getColMap()->getLocalNumElements() != dirichletColFlags.size(), std::invalid_argument,"localRowsToColumns: dirichletColFlags must be the correct size");
359 
360  LO numDirichletRows = (LO) dirichletRowIds.size();
361  LO LO_INVALID = Teuchos::OrdinalTraits<LO>::invalid();
362 
363  if(A.getCrsGraph()->getImporter().is_null()) {
364  // Serial case: If A doesn't have an importer, just set the flags from the dirichletRowIds
365  using range_type = Kokkos::RangePolicy<execution_space, LO>;
366  auto lclRowMap = A.getRowMap()->getLocalMap();
367  auto lclColMap = A.getColMap()->getLocalMap();
368 
369  Kokkos::deep_copy(execSpace,dirichletColFlags,false);
370  using range_type = Kokkos::RangePolicy<execution_space, LO>;
371  Kokkos::parallel_for
372  ("Tpetra::CrsMatrix flag Dirichlet cols",
373  range_type (execSpace, 0, numDirichletRows),
374  KOKKOS_LAMBDA (const LO i) {
375  GO row_gid = lclRowMap.getGlobalElement(dirichletRowIds[i]);
376  LO col_lid = lclColMap.getLocalElement(row_gid);
377  if(col_lid != LO_INVALID)
378  dirichletColFlags[col_lid] = true;
379  });
380  }
381  else {
382  // Parallel case: Use A's importer to halo-exchange Dirichlet information
383  auto Importer = A.getCrsGraph()->getImporter();
384  auto lclRowMap = A.getRowMap()->getLocalMap();
385  auto lclColMap = A.getColMap()->getLocalMap();
386  ::Tpetra::Vector<LO,LO,GO,NT> domainDirichlet(A.getDomainMap());
387  ::Tpetra::Vector<LO,LO,GO,NT> colDirichlet(A.getColMap());
388  const LO one = Teuchos::OrdinalTraits<LO>::one();
389  using range_type = Kokkos::RangePolicy<execution_space, LO>;
390  {
391  auto domain_data = domainDirichlet.template getLocalView<memory_space>(Access::ReadWrite);
392  Kokkos::parallel_for
393  ("Tpetra::CrsMatrix flag Dirichlet domains",
394  range_type (execSpace, 0, numDirichletRows),
395  KOKKOS_LAMBDA (const LO i) {
396  GO row_gid = lclRowMap.getGlobalElement(dirichletRowIds[i]);
397  LO col_lid = lclColMap.getLocalElement(row_gid);
398  if(col_lid != LO_INVALID)
399  domain_data(col_lid,0) = one;
400  });
401  }
402  colDirichlet.doImport(domainDirichlet,*Importer,::Tpetra::INSERT);
403  LO numCols = (LO) A.getColMap()->getLocalNumElements();
404  {
405  auto col_data = colDirichlet.template getLocalView<memory_space>(Access::ReadOnly);
406  Kokkos::parallel_for
407  ("Tpetra::CrsMatrix flag Dirichlet cols",
408  range_type (execSpace, 0, numCols),
409  KOKKOS_LAMBDA (const LO i) {
410  dirichletColFlags[i] = (col_data(i,0) == one) ? true : false;
411  });
412  }
413  }
414 }
415 
416 } // namespace Details
417 
418 template<class CrsMatrixType>
419 void
420 applyDirichletBoundaryConditionToLocalMatrixRows
421 (const typename CrsMatrixType::execution_space& execSpace,
422  CrsMatrixType& A,
423  const Kokkos::View<
424  typename CrsMatrixType::local_ordinal_type*,
425  typename CrsMatrixType::device_type> & lclRowInds)
426 {
427  using SC = typename CrsMatrixType::scalar_type;
428  using LO = typename CrsMatrixType::local_ordinal_type;
429  using GO = typename CrsMatrixType::global_ordinal_type;
430  using NT = typename CrsMatrixType::node_type;
431 
432  using local_row_indices_type =
433  Kokkos::View<const LO*, Kokkos::AnonymousSpace>;
434  const local_row_indices_type lclRowInds_a (lclRowInds);
435 
436  using Details::ApplyDirichletBoundaryConditionToLocalMatrixRows;
437  using impl_type =
438  ApplyDirichletBoundaryConditionToLocalMatrixRows<SC, LO, GO, NT>;
439  const bool runOnHost = false;
440  impl_type::run (execSpace, A, lclRowInds_a, runOnHost);
441 }
442 
443 template<class CrsMatrixType>
444 void
445 applyDirichletBoundaryConditionToLocalMatrixRows
446 (CrsMatrixType& A,
447  const Kokkos::View<
448  typename CrsMatrixType::local_ordinal_type*,
449  typename CrsMatrixType::device_type> & lclRowInds)
450 {
451  using execution_space = typename CrsMatrixType::execution_space;
452  applyDirichletBoundaryConditionToLocalMatrixRows (execution_space (), A, lclRowInds);
453 }
454 
455 template<class CrsMatrixType>
456 void
457 applyDirichletBoundaryConditionToLocalMatrixRows
458 (CrsMatrixType& A,
459  const Kokkos::View<
460  typename CrsMatrixType::local_ordinal_type*,
461  Kokkos::HostSpace> & lclRowInds)
462 {
463  using SC = typename CrsMatrixType::scalar_type;
464  using LO = typename CrsMatrixType::local_ordinal_type;
465  using GO = typename CrsMatrixType::global_ordinal_type;
466  using NT = typename CrsMatrixType::node_type;
467  using crs_matrix_type = ::Tpetra::CrsMatrix<SC, LO, GO, NT>;
468  using execution_space = typename crs_matrix_type::execution_space;
469  using memory_space = typename crs_matrix_type::device_type::memory_space;
470 
471  using Details::ApplyDirichletBoundaryConditionToLocalMatrixRows;
472  using impl_type =
473  ApplyDirichletBoundaryConditionToLocalMatrixRows<SC, LO, GO, NT>;
474 
475  // Only run on host if we can access the data
476  const bool runOnHost = Kokkos::SpaceAccessibility<Kokkos::Serial,memory_space>::accessible;
477  if(runOnHost) {
478  using local_row_indices_type = Kokkos::View<const LO*, Kokkos::AnonymousSpace>;
479  const local_row_indices_type lclRowInds_a (lclRowInds);
480  impl_type::run (execution_space (), A, lclRowInds_a, true);
481  }
482  else {
483  auto lclRowInds_a = Kokkos::create_mirror_view_and_copy(execution_space(),lclRowInds);
484  impl_type::run (execution_space (), A, lclRowInds_a, false);
485  }
486 }
487 
488 
489 template<class CrsMatrixType>
490 void
491 applyDirichletBoundaryConditionToLocalMatrixRowsAndColumns
492 (CrsMatrixType& A,
493  const Kokkos::View<
494  typename CrsMatrixType::local_ordinal_type*,
495  Kokkos::HostSpace> & lclRowInds) {
496  using SC = typename CrsMatrixType::scalar_type;
497  using LO = typename CrsMatrixType::local_ordinal_type;
498  using GO = typename CrsMatrixType::global_ordinal_type;
499  using NT = typename CrsMatrixType::node_type;
500  using crs_matrix_type = ::Tpetra::CrsMatrix<SC, LO, GO, NT>;
501  using execution_space = typename crs_matrix_type::execution_space;
502  using memory_space = typename crs_matrix_type::device_type::memory_space;
503 
504  TEUCHOS_TEST_FOR_EXCEPTION(A.getColMap().get () == nullptr, std::invalid_argument,"The matrix must have a column Map.");
505 
506  // Copy the Host array to device
507  auto lclRowInds_d = Kokkos::create_mirror_view_and_copy(execution_space(),lclRowInds);
508 
509  Kokkos::View<bool*,memory_space> dirichletColFlags("dirichletColFlags",A.getColMap()->getLocalNumElements());
510  Kokkos::View<bool*, Kokkos::AnonymousSpace> dirichletColFlags_a(dirichletColFlags);
511  Details::localRowsToColumns<SC,LO,GO,NT>(execution_space(),A,lclRowInds_d,dirichletColFlags_a);
512 
513  Details::ApplyDirichletBoundaryConditionToLocalMatrixColumns<SC, LO, GO, NT>::run(execution_space(),A,dirichletColFlags,false);
514  Details::ApplyDirichletBoundaryConditionToLocalMatrixRows<SC, LO, GO, NT>::run(execution_space(),A,lclRowInds_d,false);
515 }
516 
517 
518 template<class CrsMatrixType>
519 void
520 applyDirichletBoundaryConditionToLocalMatrixRowsAndColumns
521 (CrsMatrixType& A,
522  const Kokkos::View<
523  typename CrsMatrixType::local_ordinal_type*,
524  typename CrsMatrixType::device_type> & lclRowInds_d) {
525  using SC = typename CrsMatrixType::scalar_type;
526  using LO = typename CrsMatrixType::local_ordinal_type;
527  using GO = typename CrsMatrixType::global_ordinal_type;
528  using NT = typename CrsMatrixType::node_type;
529  using crs_matrix_type = ::Tpetra::CrsMatrix<SC, LO, GO, NT>;
530  using execution_space = typename crs_matrix_type::execution_space;
531  using memory_space = typename crs_matrix_type::device_type::memory_space;
532 
533  TEUCHOS_TEST_FOR_EXCEPTION(A.getColMap().get () == nullptr, std::invalid_argument,"The matrix must have a column Map.");
534 
535  Kokkos::View<bool*,memory_space> dirichletColFlags("dirichletColFlags",A.getColMap()->getLocalNumElements());
536  Kokkos::View<bool*, Kokkos::AnonymousSpace> dirichletColFlags_a(dirichletColFlags);
537  Details::localRowsToColumns<SC,LO,GO,NT>(execution_space(),A,lclRowInds_d,dirichletColFlags_a);
538 
539  Details::ApplyDirichletBoundaryConditionToLocalMatrixColumns<SC, LO, GO, NT>::run(execution_space(),A,dirichletColFlags,false);
540  Details::ApplyDirichletBoundaryConditionToLocalMatrixRows<SC, LO, GO, NT>::run(execution_space(),A,lclRowInds_d,false);
541 
542 }
543 
544 
545 template<class CrsMatrixType>
546 void
547 applyDirichletBoundaryConditionToLocalMatrixRowsAndColumns
548 (const typename CrsMatrixType::execution_space& execSpace,
549  CrsMatrixType& A,
550  const Kokkos::View<
551  typename CrsMatrixType::local_ordinal_type*,
552  typename CrsMatrixType::device_type> & lclRowInds_d) {
553  using SC = typename CrsMatrixType::scalar_type;
554  using LO = typename CrsMatrixType::local_ordinal_type;
555  using GO = typename CrsMatrixType::global_ordinal_type;
556  using NT = typename CrsMatrixType::node_type;
557  using crs_matrix_type = ::Tpetra::CrsMatrix<SC, LO, GO, NT>;
558  // using execution_space = typename crs_matrix_type::execution_space;
559  using memory_space = typename crs_matrix_type::device_type::memory_space;
560 
561  TEUCHOS_TEST_FOR_EXCEPTION(A.getColMap().get () == nullptr, std::invalid_argument,"The matrix must have a column Map.");
562 
563  Kokkos::View<bool*,memory_space> dirichletColFlags("dirichletColFlags",A.getColMap()->getLocalNumElements());
564  Kokkos::View<bool*, Kokkos::AnonymousSpace> dirichletColFlags_a(dirichletColFlags);
565  Details::localRowsToColumns<SC,LO,GO,NT>(execSpace,A,lclRowInds_d,dirichletColFlags_a);
566 
567  Details::ApplyDirichletBoundaryConditionToLocalMatrixColumns<SC, LO, GO, NT>::run(execSpace,A,dirichletColFlags,false);
568  Details::ApplyDirichletBoundaryConditionToLocalMatrixRows<SC, LO, GO, NT>::run(execSpace,A,lclRowInds_d,false);
569 
570 }
571 
572 
573 } // namespace Tpetra
574 
575 #endif // TPETRA_APPLYDIRICHLETBOUNDARYCONDITION_HPP
Tpetra::CrsMatrix
Sparse matrix that presents a row-oriented interface that lets users read or modify entries...
Definition: Tpetra_CrsMatrix_decl.hpp:430
Tpetra::CrsMatrix::execution_space
typename device_type::execution_space execution_space
The Kokkos execution space.
Definition: Tpetra_CrsMatrix_decl.hpp:456
Tpetra::applyDirichletBoundaryConditionToLocalMatrixRowsAndColumns
void applyDirichletBoundaryConditionToLocalMatrixRowsAndColumns(const typename CrsMatrixType::execution_space &execSpace, CrsMatrixType &A, const Kokkos::View< typename CrsMatrixType::local_ordinal_type *, typename CrsMatrixType::device_type > &lclRowInds)
For all k in [0, lclRowInds.extent(0)), set local row and column lclRowInds[k] of A to have 1 on the ...
Definition: Tpetra_applyDirichletBoundaryCondition.hpp:548
Tpetra::deep_copy
void deep_copy(MultiVector< DS, DL, DG, DN > &dst, const MultiVector< SS, SL, SG, SN > &src)
Copy the contents of the MultiVector src into dst.
Definition: Tpetra_MultiVector_decl.hpp:2488
Tpetra::INSERT
Insert new values that don&#39;t currently exist.
Definition: Tpetra_CombineMode.hpp:99
Tpetra::CrsMatrix::impl_scalar_type
typename row_matrix_type::impl_scalar_type impl_scalar_type
The type used internally in place of Scalar.
Definition: Tpetra_CrsMatrix_decl.hpp:487
Tpetra::Vector
A distributed dense vector.
Definition: Tpetra_Vector_decl.hpp:76
Tpetra::applyDirichletBoundaryConditionToLocalMatrixRows
void applyDirichletBoundaryConditionToLocalMatrixRows(const typename CrsMatrixType::execution_space &execSpace, CrsMatrixType &A, const Kokkos::View< typename CrsMatrixType::local_ordinal_type *, typename CrsMatrixType::device_type > &lclRowInds)
For all k in [0, lclRowInds.extent(0)), set local row lclRowInds[k] of A to have 1 on the diagonal an...
Definition: Tpetra_applyDirichletBoundaryCondition.hpp:421
Generated on Fri Apr 26 2024 09:19:49 for Tpetra parallel linear algebra by   doxygen 1.8.5