EpetraExt Package Browser (Single Doxygen Collection) Development
EpetraExt_SolverMap_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_SolverMap_CrsMatrix.h>
00043 
00044 #include <Epetra_CrsGraph.h>
00045 #include <Epetra_CrsMatrix.h>
00046 #include <Epetra_Map.h>
00047 #include <Epetra_Comm.h>
00048 
00049 #include <vector>
00050 
00051 namespace EpetraExt {
00052 
00053 CrsMatrix_SolverMap::
00054 ~CrsMatrix_SolverMap()
00055 {
00056   if( newObj_ && newObj_ != origObj_ ) delete newObj_;
00057   if( NewGraph_ ) delete NewGraph_;
00058   if( NewColMap_ ) delete NewColMap_;
00059 }
00060 
00061   template<typename int_type>
00062 CrsMatrix_SolverMap::NewTypeRef
00063 CrsMatrix_SolverMap::
00064 construct( OriginalTypeRef orig )
00065 {
00066   origObj_ = &orig;
00067 
00068   assert( !orig.IndicesAreGlobal() );
00069 
00070   //test if matrix has missing local columns in its col std::map
00071   const Epetra_Map & RowMap = orig.RowMap();
00072   const Epetra_Map & DomainMap = orig.DomainMap();
00073   const Epetra_Map & ColMap = orig.ColMap();
00074   const Epetra_Comm & Comm = RowMap.Comm();
00075   int NumMyRows = RowMap.NumMyElements();
00076   int NumCols = DomainMap.NumMyElements();
00077   int Match = 0;
00078   for( int i = 0; i < NumCols; ++i )
00079     if( DomainMap.GID64(i) != ColMap.GID64(i) )
00080     {
00081       Match = 1;
00082       break;
00083     }
00084 
00085   int MatchAll = 0;
00086   Comm.SumAll( &Match, &MatchAll, 1 );
00087 
00088   if( !MatchAll )
00089   {
00090     newObj_ = origObj_;
00091   }
00092   else
00093   {
00094     //create ColMap with all local rows included
00095     std::vector<int_type> Cols(NumCols);
00096     //fill Cols list with GIDs of all local columns 
00097     for( int i = 0; i < NumCols; ++i )
00098       Cols[i] = (int_type) DomainMap.GID64(i);
00099 
00100     //now append to Cols any ghost column entries
00101     int NumMyCols = ColMap.NumMyElements();
00102     for( int i = 0; i < NumMyCols; ++i )
00103       if( !DomainMap.MyGID( ColMap.GID64(i) ) ) Cols.push_back( (int_type) ColMap.GID64(i) );
00104     
00105     int NewNumMyCols = Cols.size();
00106     int NewNumGlobalCols;
00107     Comm.SumAll( &NewNumMyCols, &NewNumGlobalCols, 1 );
00108     //create new column std::map
00109     NewColMap_ = new Epetra_Map( NewNumGlobalCols, NewNumMyCols,&Cols[0], DomainMap.IndexBase64(), Comm );
00110 
00111     //New Graph
00112     std::vector<int> NumIndicesPerRow( NumMyRows );
00113     for( int i = 0; i < NumMyRows; ++i )
00114       NumIndicesPerRow[i] = orig.NumMyEntries(i);
00115     NewGraph_ = new Epetra_CrsGraph( Copy, RowMap, *NewColMap_, &NumIndicesPerRow[0] );
00116 
00117     int MaxNumEntries = orig.MaxNumEntries();
00118     int NumEntries;
00119     std::vector<int_type> Indices( MaxNumEntries );
00120     for( int i = 0; i < NumMyRows; ++i )
00121     {
00122       int_type RowGID = (int_type) RowMap.GID64(i);
00123       orig.Graph().ExtractGlobalRowCopy( RowGID, MaxNumEntries, NumEntries, &Indices[0] );
00124       NewGraph_->InsertGlobalIndices( RowGID, NumEntries, &Indices[0] );
00125     }
00126     const Epetra_Map & RangeMap = orig.RangeMap();
00127     NewGraph_->FillComplete(DomainMap,RangeMap);
00128 
00129     //intial construction of matrix 
00130     Epetra_CrsMatrix * NewMatrix = new Epetra_CrsMatrix( View, *NewGraph_ );
00131 
00132     //insert views of row values
00133     int * myIndices;
00134     double * myValues;
00135     int indicesCnt;
00136     int numMyRows = NewMatrix->NumMyRows();
00137     for( int i = 0; i < numMyRows; ++i )
00138     {
00139       orig.ExtractMyRowView( i, indicesCnt, myValues, myIndices );
00140       NewGraph_->ExtractMyRowView( i, indicesCnt, myIndices );
00141 
00142       NewMatrix->InsertMyValues( i, indicesCnt, myValues, myIndices );
00143     }
00144 
00145     NewMatrix->FillComplete(DomainMap,RangeMap);
00146 
00147     newObj_ = NewMatrix;
00148   }
00149 
00150   return *newObj_;
00151 }
00152 
00153 CrsMatrix_SolverMap::NewTypeRef
00154 CrsMatrix_SolverMap::
00155 operator()( OriginalTypeRef orig )
00156 {
00157 #ifndef EPETRA_NO_32BIT_GLOBAL_INDICES
00158   if(orig.RowMap().GlobalIndicesInt()) {
00159     return construct<int>(orig);
00160   }
00161   else
00162 #endif
00163 #ifndef EPETRA_NO_64BIT_GLOBAL_INDICES
00164   if(orig.RowMap().GlobalIndicesLongLong()) {
00165     return construct<long long>(orig);
00166   }
00167   else
00168 #endif
00169     throw "CrsMatrix_SolverMap::operator(): GlobalIndices type unknown";
00170 }
00171 } // namespace EpetraExt
00172 
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines