EpetraExt Development
Loading...
Searching...
No Matches
EpetraExt_Reindex_CrsMatrix.cpp
Go to the documentation of this file.
1//@HEADER
2// ***********************************************************************
3//
4// EpetraExt: Epetra Extended - Linear Algebra Services Package
5// Copyright (2011) Sandia Corporation
6//
7// Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
8// the U.S. Government retains certain rights in this software.
9//
10// Redistribution and use in source and binary forms, with or without
11// modification, are permitted provided that the following conditions are
12// met:
13//
14// 1. Redistributions of source code must retain the above copyright
15// notice, this list of conditions and the following disclaimer.
16//
17// 2. Redistributions in binary form must reproduce the above copyright
18// notice, this list of conditions and the following disclaimer in the
19// documentation and/or other materials provided with the distribution.
20//
21// 3. Neither the name of the Corporation nor the names of the
22// contributors may be used to endorse or promote products derived from
23// this software without specific prior written permission.
24//
25// THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
26// EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
27// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
28// PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
29// CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
30// EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
31// PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
32// PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
33// LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
34// NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
35// SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
36//
37// Questions? Contact Michael A. Heroux (maherou@sandia.gov)
38//
39// ***********************************************************************
40//@HEADER
41
42#include <Epetra_ConfigDefs.h>
44
45#include <vector>
46
47#include <Epetra_CrsMatrix.h>
48#include <Epetra_Map.h>
49#include <Epetra_GIDTypeVector.h>
50
51#ifndef EPETRA_NO_64BIT_GLOBAL_INDICES
52#include <Epetra_LongLongVector.h>
53#endif
54
55#include <Epetra_Export.h>
56#include <Epetra_Import.h>
57
58namespace EpetraExt {
59
62{
63 if( newObj_ ) delete newObj_;
64 if( NewColMap_ ) delete NewColMap_;
65}
66
67template<typename int_type>
69CrsMatrix_Reindex::
70Toperator( OriginalTypeRef orig )
71{
72 origObj_ = &orig;
73
74 //test std::map, must have same number of local and global elements as original row std::map
75 //Epetra_Map & OldRowMap = const_cast<Epetra_Map&>(orig.RowMap());
76 Epetra_Map & OldDomainMap = const_cast<Epetra_Map&>(orig.OperatorDomainMap());
77 Epetra_Map & OldColMap = const_cast<Epetra_Map&>(orig.ColMap());
78 int NumMyElements = OldDomainMap.NumMyElements();
79 int_type NumGlobalElements = (int_type) OldDomainMap.NumGlobalElements64();
80 assert( orig.RowMap().NumMyElements() == NewRowMap_.NumMyElements() );
81
82 if (NumGlobalElements == 0 && orig.RowMap().NumGlobalElements64() == 0 )
83 {
84 //construct a zero matrix as a placeholder, don't do reindexing analysis.
85 Epetra_CrsMatrix * NewMatrix = new Epetra_CrsMatrix( View, orig.RowMap(), orig.ColMap(), 0 );
86 newObj_ = NewMatrix;
87 }
88 else {
89
90 //Construct new Column Map
91 typename Epetra_GIDTypeVector<int_type>::impl Cols( OldDomainMap );
92 typename Epetra_GIDTypeVector<int_type>::impl NewCols( OldColMap );
93 Epetra_Import Importer( OldColMap, OldDomainMap );
94
95 Epetra_Map tmpColMap( NumGlobalElements, NumMyElements, 0, OldDomainMap.Comm() );
96
97 for( int i = 0; i < NumMyElements; ++i )
98 Cols[i] = (int_type) tmpColMap.GID64(i);
99
100 NewCols.Import( Cols, Importer, Insert );
101
102 std::vector<int_type*> NewColIndices(1);
103 NewCols.ExtractView( &NewColIndices[0] );
104
105 int NumMyColElements = OldColMap.NumMyElements();
106 int_type NumGlobalColElements = (int_type) OldColMap.NumGlobalElements64();
107
108 NewColMap_ = new Epetra_Map( NumGlobalColElements, NumMyColElements, NewColIndices[0], (int_type) OldColMap.IndexBase64(), OldColMap.Comm() );
109
110 //intial construction of matrix
111 Epetra_CrsMatrix * NewMatrix = new Epetra_CrsMatrix( View, NewRowMap_, *NewColMap_, 0 );
112
113 //insert views of row values
114 int * myIndices;
115 double * myValues;
116 int indicesCnt;
117 int numMyRows = NewMatrix->NumMyRows();
118 for( int i = 0; i < numMyRows; ++i )
119 {
120 orig.ExtractMyRowView( i, indicesCnt, myValues, myIndices );
121 NewMatrix->InsertMyValues( i, indicesCnt, myValues, myIndices );
122 }
123
124 NewMatrix->FillComplete();
125
126 newObj_ = NewMatrix;
127
128 }
129
130 return *newObj_;
131}
132
135operator()( OriginalTypeRef orig )
136{
137#ifndef EPETRA_NO_32BIT_GLOBAL_INDICES
138 if(orig.RowMatrixRowMap().GlobalIndicesInt())
139 return Toperator<int>(orig);
140 else
141#endif
142#ifndef EPETRA_NO_64BIT_GLOBAL_INDICES
143 if(orig.RowMatrixRowMap().GlobalIndicesLongLong())
144 return Toperator<long long>(orig);
145 else
146#endif
147 throw "EpetraExt::CrsMatrix_Reindex::operator(): Global indices unknown.";
148}
149
150} // namespace EpetraExt
151
NewTypeRef operator()(OriginalTypeRef orig)
Constructs "reindexed" Matrix.
const Epetra_Comm & Comm() const
int NumMyElements() const
int FillComplete(bool OptimizeDataStorage=true)
int NumMyRows() const
int InsertMyValues(int MyRow, int NumEntries, const double *Values, const int *Indices)
EpetraExt::BlockCrsMatrix: A class for constructing a distributed block matrix.