EpetraExt Development
EpetraExt_SubCopy_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 "EpetraExt_SubCopy_CrsMatrix.h"
00043 
00044 #include <Epetra_CrsMatrix.h>
00045 #include <Epetra_Map.h>
00046 #include <Epetra_Import.h>
00047 #include <Epetra_IntSerialDenseVector.h>
00048 #include <Epetra_LongLongSerialDenseVector.h>
00049 #include <Epetra_GIDTypeSerialDenseVector.h>
00050 #include <vector>
00051 
00052 namespace EpetraExt {
00053 
00054 CrsMatrix_SubCopy::
00055 ~CrsMatrix_SubCopy()
00056 {
00057   if( newObj_ ) delete newObj_;
00058 }
00059 
00060 template<typename int_type>
00061 CrsMatrix_SubCopy::NewTypeRef
00062 CrsMatrix_SubCopy::
00063 transform( OriginalTypeRef orig )
00064 {
00065   origObj_ = &orig;
00066 
00067   //Error, must be local indices
00068   assert( orig.Filled() );
00069 
00070   //test maps, new map must be subset of old
00071   const Epetra_Map & oRowMap = orig.RowMap();
00072   const Epetra_Map & oColMap = orig.ColMap();
00073 
00074   int oNumRows = oRowMap.NumMyElements();
00075   (void) oNumRows; // Silence "unused variable" compiler warning.
00076   int oNumCols = oColMap.NumMyElements();
00077   int nNumRows = newRowMap_.NumMyElements();
00078   int nNumDomain = newDomainMap_.NumMyElements();
00079 
00080   bool matched = true;
00081 
00082   // Make sure all rows in newRowMap are already on this processor
00083   for( int i = 0; i < nNumRows; ++i )
00084     matched = matched && ( oRowMap.MyGID(newRowMap_.GID64(i)) );
00085   if( !matched ) std::cerr << "EDT_CrsMatrix_SubCopy: Bad new_row_Map.  GIDs of new row map must be GIDs of the original row map on the same processor.\n";
00086 
00087   // Make sure all GIDs in the new domain map are GIDs in the old domain map
00088   if( !newRangeMap_.SameAs(newDomainMap_) ) {
00089     Epetra_IntSerialDenseVector pidList(nNumDomain);
00090   int_type* newDomainMap_MyGlob = 0;
00091   newDomainMap_.MyGlobalElementsPtr(newDomainMap_MyGlob);
00092     oColMap.RemoteIDList(newDomainMap_.NumMyElements(), newDomainMap_MyGlob, pidList.Values(), 0);
00093     for( int i = 0; i < nNumDomain; ++i )
00094       matched = matched && ( pidList[i]>=0 );
00095   }
00096 
00097   if( !matched ) std::cout << "EDT_CrsMatrix_SubCopy: Bad newDomainMap.  One or more GIDs in new domain map are not part of original domain map.\n";
00098   assert( matched );
00099 
00100 
00101   // Next build new column map
00102   Epetra_IntSerialDenseVector pidList(oNumCols);
00103   Epetra_IntSerialDenseVector lidList(oNumCols);
00104   Epetra_IntSerialDenseVector sizeList(oNumCols);
00105   int_type* oColMap_MyGlob = 0;
00106   oColMap.MyGlobalElementsPtr(oColMap_MyGlob);
00107   newDomainMap_.RemoteIDList(oColMap.NumMyElements(), oColMap_MyGlob, pidList.Values(), 0);
00108   int numNewCols = 0;
00109   typename Epetra_GIDTypeSerialDenseVector<int_type>::impl newColMapGidList(oNumCols);
00110   int_type * origColGidList = 0;
00111   oColMap.MyGlobalElementsPtr(origColGidList);
00112   for( int i = 0; i < oNumCols; ++i )
00113     if (pidList[i] >=0)
00114       newColMapGidList[numNewCols++]= origColGidList[i];
00115   newColMap_ = Epetra_Map(-1, numNewCols, newColMapGidList.Values(), 0, oColMap.Comm());
00116 
00117   importer_ = new Epetra_Import(newRowMap_, oRowMap);
00118 
00119   Epetra_CrsMatrix * newMatrix = new Epetra_CrsMatrix(Copy, newRowMap_, newColMap_, 0);
00120 
00121   newObj_ = newMatrix;
00122 
00123   newObj_->Import(*origObj_, *importer_, Add);
00124 
00125   newObj_->FillComplete();
00126 
00127   return *newObj_;
00128 }
00129 
00130 //==============================================================================
00131 
00132 CrsMatrix_SubCopy::NewTypeRef
00133 CrsMatrix_SubCopy::
00134 operator()( OriginalTypeRef orig )
00135 {
00136   const Epetra_Map & oRowMap = orig.RowMap();
00137   const Epetra_Map & oColMap = orig.ColMap();
00138 
00139 #ifndef EPETRA_NO_32BIT_GLOBAL_INDICES
00140   if(oRowMap.GlobalIndicesInt() && oColMap.GlobalIndicesInt()) {
00141     return transform<int>(orig);
00142   }
00143   else
00144 #endif
00145 #ifndef EPETRA_NO_64BIT_GLOBAL_INDICES
00146   if(oRowMap.GlobalIndicesLongLong() && oColMap.GlobalIndicesLongLong()) {
00147     return transform<long long>(orig);
00148   }
00149   else
00150 #endif
00151     throw "CrsMatrix_SubCopy::operator(): GlobalIndices type unknown";
00152 }
00153 
00154 //==============================================================================
00155 bool CrsMatrix_SubCopy::fwd()
00156 {
00157 
00158   if (newObj_->Filled()) newObj_->PutScalar(0.0); // zero contents
00159 
00160   newObj_->Import(*origObj_, *importer_, Add);
00161 
00162   newObj_->FillComplete();
00163 
00164 
00165   return (true);
00166 }
00167 
00168 //==============================================================================
00169 bool CrsMatrix_SubCopy::rvs()
00170 {
00171   if (!newObj_->Filled()) return(false); // Must have fillCompleted
00172 
00173   origObj_->Export(*newObj_, *importer_, Add);
00174 
00175   origObj_->FillComplete();
00176 
00177   return (true);
00178 }
00179 
00180 } // namespace EpetraExt
00181 
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines