EpetraExt Development
EpetraExt_Reindex_CrsMatrix.cpp
Go to the documentation of this file.
00001 //@HEADER
00002 // ***********************************************************************
00003 //
00004 //     EpetraExt: Epetra Extended - Linear Algebra Services Package
00005 //                 Copyright (2011) Sandia Corporation
00006 //
00007 // Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
00008 // the U.S. Government retains certain rights in this software.
00009 //
00010 // Redistribution and use in source and binary forms, with or without
00011 // modification, are permitted provided that the following conditions are
00012 // met:
00013 //
00014 // 1. Redistributions of source code must retain the above copyright
00015 // notice, this list of conditions and the following disclaimer.
00016 //
00017 // 2. Redistributions in binary form must reproduce the above copyright
00018 // notice, this list of conditions and the following disclaimer in the
00019 // documentation and/or other materials provided with the distribution.
00020 //
00021 // 3. Neither the name of the Corporation nor the names of the
00022 // contributors may be used to endorse or promote products derived from
00023 // this software without specific prior written permission.
00024 //
00025 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
00026 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
00027 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
00028 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
00029 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
00030 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
00031 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
00032 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
00033 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
00034 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
00035 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
00036 //
00037 // Questions? Contact Michael A. Heroux (maherou@sandia.gov)
00038 //
00039 // ***********************************************************************
00040 //@HEADER
00041 
00042 #include <Epetra_ConfigDefs.h>
00043 #include <EpetraExt_Reindex_CrsMatrix.h>
00044 
00045 #include <vector>
00046 
00047 #include <Epetra_CrsMatrix.h>
00048 #include <Epetra_Map.h>
00049 #include <Epetra_GIDTypeVector.h>
00050 
00051 #ifndef EPETRA_NO_64BIT_GLOBAL_INDICES
00052 #include <Epetra_LongLongVector.h>
00053 #endif
00054 
00055 #include <Epetra_Export.h>
00056 #include <Epetra_Import.h>
00057 
00058 namespace EpetraExt {
00059 
00060 CrsMatrix_Reindex::
00061 ~CrsMatrix_Reindex()
00062 {
00063   if( newObj_ ) delete newObj_;
00064   if( NewColMap_ ) delete NewColMap_;
00065 }
00066 
00067 template<typename int_type>
00068 CrsMatrix_Reindex::NewTypeRef
00069 CrsMatrix_Reindex::
00070 Toperator( OriginalTypeRef orig )
00071 {
00072   origObj_ = &orig;
00073 
00074   //test std::map, must have same number of local and global elements as original row std::map
00075   //Epetra_Map & OldRowMap = const_cast<Epetra_Map&>(orig.RowMap());
00076   Epetra_Map & OldDomainMap = const_cast<Epetra_Map&>(orig.OperatorDomainMap());
00077   Epetra_Map & OldColMap = const_cast<Epetra_Map&>(orig.ColMap());
00078   int NumMyElements = OldDomainMap.NumMyElements();
00079   int_type NumGlobalElements = (int_type) OldDomainMap.NumGlobalElements64();
00080   assert( orig.RowMap().NumMyElements() == NewRowMap_.NumMyElements() );
00081 
00082   if (NumGlobalElements == 0 && orig.RowMap().NumGlobalElements64() == 0 )
00083   {
00084     //construct a zero matrix as a placeholder, don't do reindexing analysis.
00085     Epetra_CrsMatrix * NewMatrix = new Epetra_CrsMatrix( View, orig.RowMap(), orig.ColMap(), 0 );
00086     newObj_ = NewMatrix;
00087   }
00088   else {
00089 
00090     //Construct new Column Map
00091     typename Epetra_GIDTypeVector<int_type>::impl Cols( OldDomainMap );
00092     typename Epetra_GIDTypeVector<int_type>::impl NewCols( OldColMap );
00093     Epetra_Import Importer( OldColMap, OldDomainMap );
00094  
00095     Epetra_Map tmpColMap( NumGlobalElements, NumMyElements, 0, OldDomainMap.Comm() );
00096  
00097     for( int i = 0; i < NumMyElements; ++i )
00098       Cols[i] = (int_type) tmpColMap.GID64(i);
00099 
00100     NewCols.Import( Cols, Importer, Insert );
00101 
00102     std::vector<int_type*> NewColIndices(1);
00103     NewCols.ExtractView( &NewColIndices[0] );
00104 
00105     int NumMyColElements = OldColMap.NumMyElements();
00106     int_type NumGlobalColElements = (int_type) OldColMap.NumGlobalElements64();
00107 
00108     NewColMap_ = new Epetra_Map( NumGlobalColElements, NumMyColElements, NewColIndices[0], (int_type) OldColMap.IndexBase64(), OldColMap.Comm() );
00109 
00110     //intial construction of matrix 
00111     Epetra_CrsMatrix * NewMatrix = new Epetra_CrsMatrix( View, NewRowMap_, *NewColMap_, 0 );
00112 
00113     //insert views of row values
00114     int * myIndices;
00115     double * myValues;
00116     int indicesCnt;
00117     int numMyRows = NewMatrix->NumMyRows();
00118     for( int i = 0; i < numMyRows; ++i )
00119     {
00120       orig.ExtractMyRowView( i, indicesCnt, myValues, myIndices );
00121       NewMatrix->InsertMyValues( i, indicesCnt, myValues, myIndices );
00122     }
00123 
00124     NewMatrix->FillComplete();
00125 
00126     newObj_ = NewMatrix;
00127 
00128   }
00129 
00130   return *newObj_;
00131 }
00132 
00133 CrsMatrix_Reindex::NewTypeRef
00134 CrsMatrix_Reindex::
00135 operator()( OriginalTypeRef orig )
00136 {
00137 #ifndef EPETRA_NO_32BIT_GLOBAL_INDICES
00138   if(orig.RowMatrixRowMap().GlobalIndicesInt())
00139     return Toperator<int>(orig);
00140   else
00141 #endif
00142 #ifndef EPETRA_NO_64BIT_GLOBAL_INDICES
00143   if(orig.RowMatrixRowMap().GlobalIndicesLongLong())
00144     return Toperator<long long>(orig);
00145   else
00146 #endif
00147     throw "EpetraExt::CrsMatrix_Reindex::operator(): Global indices unknown.";
00148 }
00149 
00150 } // namespace EpetraExt
00151 
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines