EpetraExt Development
Loading...
Searching...
No Matches
EpetraExt_RestrictedMultiVectorWrapper.cpp
Go to the documentation of this file.
1/*
2//@HEADER
3// ***********************************************************************
4//
5// EpetraExt: Epetra Extended - Linear Algebra Services 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//@HEADER
42*/
43
45
46
47#ifdef HAVE_MPI
48
49
51#include "Epetra_MpiComm.h"
52#include "Epetra_BlockMap.h"
53#include "Epetra_MultiVector.h"
54
55
56namespace EpetraExt {
57
58
59RestrictedMultiVectorWrapper::RestrictedMultiVectorWrapper()
60 : proc_is_active(true),
61 subcomm_is_set(false),
62 MPI_SubComm_(MPI_COMM_NULL),
63 RestrictedComm_(0),
64 ResMap_(0),
65 input_mv_(),
66 restricted_mv_()
67 {}
68
69
70RestrictedMultiVectorWrapper::~RestrictedMultiVectorWrapper() {
71 delete ResMap_;
72 delete RestrictedComm_;
73}
74
75
76int RestrictedMultiVectorWrapper::SetMPISubComm(MPI_Comm MPI_SubComm){
77 if(!subcomm_is_set){
78 MPI_SubComm_=MPI_SubComm; delete RestrictedComm_; subcomm_is_set=true;
79 return 0;
80 }
81 else return -1;
82}
83
84
85int
86RestrictedMultiVectorWrapper::
87restrict_comm (Teuchos::RCP<Epetra_MultiVector> input_mv)
88{
89 using Teuchos::rcp;
90 input_mv_ = input_mv;
91
92 // Extract the input MV's communicator and Map.
93 const Epetra_MpiComm *InComm =
94 dynamic_cast<const Epetra_MpiComm*> (& (input_mv_->Comm ()));
95 const Epetra_BlockMap *InMap =
96 dynamic_cast<const Epetra_BlockMap*> (& (input_mv_->Map ()));
97
98 if (! InComm || ! InMap) {
99 return -1; // At least one dynamic cast failed.
100 }
101
102 if (! subcomm_is_set) {
103 /* Build the Split Communicators, If Needed */
104 int color;
105 if (InMap->NumMyElements()) {
106 color = 1;
107 }
108 else {
109 color = MPI_UNDEFINED;
110 }
111 const int err =
112 MPI_Comm_split (InComm->Comm(), color, InComm->MyPID(), &MPI_SubComm_);
113 if (err != MPI_SUCCESS) {
114 return -2;
115 }
116 }
117 else {
118 /* Sanity check user-provided subcomm - drop an error if the MPISubComm
119 does not include a processor with data. */
120 if (input_mv->MyLength() && MPI_SubComm_ == MPI_COMM_NULL) {
121 return -2;
122 }
123 }
124
125 /* Mark active processors */
126 if (MPI_SubComm_ == MPI_COMM_NULL) {
127 proc_is_active = false;
128 }
129 else {
130 proc_is_active = true;
131 }
132
133 if (proc_is_active) {
134
135#ifndef EPETRA_NO_32BIT_GLOBAL_INDICES
136 if(InMap->GlobalIndicesInt()) {
137 int Nrows = InMap->NumGlobalElements ();
138 RestrictedComm_ = new Epetra_MpiComm (MPI_SubComm_);
139
140 // Build the restricted Maps
141 ResMap_ = new Epetra_BlockMap (Nrows, InMap->NumMyElements(),
142 InMap->MyGlobalElements(),
143 InMap->ElementSizeList(),
144 InMap->IndexBase(), *RestrictedComm_);
145 }
146 else
147#endif
148#ifndef EPETRA_NO_64BIT_GLOBAL_INDICES
149 if(InMap->GlobalIndicesLongLong()) {
150 long long Nrows = InMap->NumGlobalElements64 ();
151 RestrictedComm_ = new Epetra_MpiComm (MPI_SubComm_);
152
153 // Build the restricted Maps
154 ResMap_ = new Epetra_BlockMap (Nrows, InMap->NumMyElements(),
155 InMap->MyGlobalElements64(),
156 InMap->ElementSizeList(),
157 InMap->IndexBase64(), *RestrictedComm_);
158 }
159 else
160#endif
161 throw "EpetraExt::RestrictedMultiVectorWrapper::restrict_comm ERROR: GlobalIndices type unknown";
162
163 // Allocate the restricted matrix
164 double *A;
165 int LDA;
166 input_mv_->ExtractView (&A,&LDA);
167 restricted_mv_ = rcp (new Epetra_MultiVector (View, *ResMap_, A, LDA,
168 input_mv_->NumVectors ()));
169 }
170 return 0; // Success!
171}/*end restrict_comm*/
172
173}
174
175#endif
int MyGlobalElements(int *MyGlobalElementList) const
int * ElementSizeList() const
int IndexBase() const
bool GlobalIndicesInt() const
int NumGlobalElements() const
int NumMyElements() const
bool GlobalIndicesLongLong() const
MPI_Comm Comm() const
int MyPID() const
EpetraExt::BlockCrsMatrix: A class for constructing a distributed block matrix.