EpetraExt Package Browser (Single Doxygen Collection) Development
Loading...
Searching...
No Matches
EpetraExt_RestrictedCrsMatrixWrapper.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#include "Epetra_ConfigDefs.h"
46
47
48#ifdef HAVE_MPI
49
50
52#include "Epetra_MpiComm.h"
53#include "Epetra_Map.h"
54#include "Epetra_CrsMatrix.h"
55
56
57namespace EpetraExt{
58
59RestrictedCrsMatrixWrapper::RestrictedCrsMatrixWrapper()
60 : proc_is_active(true),
61 subcomm_is_set(false),
62 MPI_SubComm_(MPI_COMM_NULL),
63 RestrictedComm_(0),
64 ResRowMap_(0),
65 ResColMap_(0),
66 input_matrix_(),
67 restricted_matrix_()
68{
69
70}
71
72RestrictedCrsMatrixWrapper::~RestrictedCrsMatrixWrapper() {
73 delete ResRowMap_;
74 delete ResColMap_;
75 delete RestrictedComm_;
76}
77
78
79
80int RestrictedCrsMatrixWrapper::SetMPISubComm(MPI_Comm MPI_SubComm){
81 if(!subcomm_is_set){
82 MPI_SubComm_=MPI_SubComm; delete RestrictedComm_; subcomm_is_set=true;
83 return 0;
84 }
85 else return -1;
86}
87
88
89template<typename int_type>
90int RestrictedCrsMatrixWrapper::Trestrict_comm(Teuchos::RCP<Epetra_CrsMatrix> input_matrix){
91 /* Pull the Matrix Info */
92 input_matrix_=input_matrix;
93
94 const Epetra_MpiComm *InComm = dynamic_cast<const Epetra_MpiComm*>(& input_matrix_->Comm());
95 const Epetra_Map *InRowMap= dynamic_cast<const Epetra_Map* >(& input_matrix_->RowMap());
96 const Epetra_Map *InColMap= dynamic_cast<const Epetra_Map* >(& input_matrix_->ColMap());
97
98 if(!InComm || !InRowMap || !InColMap) return (-1);
99
100 int_type Nrows = (int_type) InRowMap->NumGlobalElements64();
101 int_type Ncols = (int_type) InColMap->NumGlobalElements64();
102
103 if(!subcomm_is_set){
104 /* Build the Split Communicators, If Needed */
105 int color;
106 if(InRowMap->NumMyElements()) color=1;
107 else color=MPI_UNDEFINED;
108 MPI_Comm_split(InComm->Comm(),color,InComm->MyPID(),&MPI_SubComm_);
109 }
110 else{
111 /* Sanity check user-provided subcomm - drop an error if the MPISubComm
112 does not include a processor with data. */
113 if (input_matrix->NumMyRows() && MPI_SubComm_ == MPI_COMM_NULL)
114 return(-2);
115 }
116
117 /* Mark active processors */
118 if(MPI_SubComm_ == MPI_COMM_NULL) proc_is_active=false;
119 else proc_is_active=true;
120
121
122 if(proc_is_active){
123 RestrictedComm_=new Epetra_MpiComm(MPI_SubComm_);
124
125 int_type* RowMapGlobalElements = 0;
126 InRowMap->MyGlobalElementsPtr(RowMapGlobalElements);
127 int_type* ColMapGlobalElements = 0;
128 InColMap->MyGlobalElementsPtr(ColMapGlobalElements);
129
130 /* Build the Restricted Maps */
131 ResRowMap_ = new Epetra_Map(Nrows,InRowMap->NumMyElements(),RowMapGlobalElements,
132 (int_type) InRowMap->IndexBase64(),*RestrictedComm_);
133 ResColMap_ = new Epetra_Map(Ncols,InColMap->NumMyElements(),ColMapGlobalElements,
134 (int_type) InColMap->IndexBase64(),*RestrictedComm_);
135
136 int *colind,Nr;
137 double *values;
138
139 /* Allocate the Restricted Matrix */
140 restricted_matrix_= Teuchos::rcp(new Epetra_CrsMatrix(View,*ResRowMap_,*ResColMap_,0));
141 for(int i=0;i<input_matrix_->NumMyRows();i++) {
142 input_matrix_->ExtractMyRowView(i,Nr,values,colind);
143 restricted_matrix_->InsertMyValues(i,Nr,values,colind);
144 }
145 restricted_matrix_->FillComplete();
146 }
147
148 return 0;
149}/*end restrict_comm*/
150
151int RestrictedCrsMatrixWrapper::restrict_comm(Teuchos::RCP<Epetra_CrsMatrix> input_matrix)
152{
153#ifndef EPETRA_NO_32BIT_GLOBAL_INDICES
154 if(input_matrix->RowMap().GlobalIndicesInt()) {
155 return Trestrict_comm<int>(input_matrix);
156 }
157 else
158#endif
159#ifndef EPETRA_NO_64BIT_GLOBAL_INDICES
160 if(input_matrix->RowMap().GlobalIndicesLongLong()) {
161 return Trestrict_comm<long long>(input_matrix);
162 }
163 else
164#endif
165 throw "EpetraExt::Trestrict_comm: ERROR, GlobalIndices type unknown.";
166}
167
168
169}
170#endif
long long NumGlobalElements64() const
EpetraExt::BlockCrsMatrix: A class for constructing a distributed block matrix.