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 CrsMatrix_SolverMap::NewTypeRef
00062 CrsMatrix_SolverMap::
00063 operator()( OriginalTypeRef orig )
00064 {
00065   origObj_ = &orig;
00066 
00067   assert( !orig.IndicesAreGlobal() );
00068 
00069   //test if matrix has missing local columns in its col std::map
00070   const Epetra_Map & RowMap = orig.RowMap();
00071   const Epetra_Map & DomainMap = orig.DomainMap();
00072   const Epetra_Map & ColMap = orig.ColMap();
00073   const Epetra_Comm & Comm = RowMap.Comm();
00074   int NumMyRows = RowMap.NumMyElements();
00075   int NumCols = DomainMap.NumMyElements();
00076   int Match = 0;
00077   for( int i = 0; i < NumCols; ++i )
00078     if( DomainMap.GID(i) != ColMap.GID(i) )
00079     {
00080       Match = 1;
00081       break;
00082     }
00083 
00084   int MatchAll = 0;
00085   Comm.SumAll( &Match, &MatchAll, 1 );
00086 
00087   if( !MatchAll )
00088   {
00089     newObj_ = origObj_;
00090   }
00091   else
00092   {
00093     //create ColMap with all local rows included
00094     std::vector<int> Cols(NumCols);
00095     //fill Cols list with GIDs of all local columns 
00096     for( int i = 0; i < NumCols; ++i )
00097       Cols[i] = DomainMap.GID(i);
00098 
00099     //now append to Cols any ghost column entries
00100     int NumMyCols = ColMap.NumMyElements();
00101     for( int i = 0; i < NumMyCols; ++i )
00102       if( !DomainMap.MyGID( ColMap.GID(i) ) ) Cols.push_back( ColMap.GID(i) );
00103     
00104     int NewNumMyCols = Cols.size();
00105     int NewNumGlobalCols;
00106     Comm.SumAll( &NewNumMyCols, &NewNumGlobalCols, 1 );
00107     //create new column std::map
00108     NewColMap_ = new Epetra_Map( NewNumGlobalCols, NewNumMyCols,&Cols[0], DomainMap.IndexBase(), Comm );
00109 
00110     //New Graph
00111     std::vector<int> NumIndicesPerRow( NumMyRows );
00112     for( int i = 0; i < NumMyRows; ++i )
00113       NumIndicesPerRow[i] = orig.NumMyEntries(i);
00114     NewGraph_ = new Epetra_CrsGraph( Copy, RowMap, *NewColMap_, &NumIndicesPerRow[0] );
00115 
00116     int MaxNumEntries = orig.MaxNumEntries();
00117     int NumEntries;
00118     std::vector<int> Indices( MaxNumEntries );
00119     for( int i = 0; i < NumMyRows; ++i )
00120     {
00121       int RowGID = RowMap.GID(i);
00122       orig.Graph().ExtractGlobalRowCopy( RowGID, MaxNumEntries, NumEntries, &Indices[0] );
00123       NewGraph_->InsertGlobalIndices( RowGID, NumEntries, &Indices[0] );
00124     }
00125     const Epetra_Map & RangeMap = orig.RangeMap();
00126     NewGraph_->FillComplete(DomainMap,RangeMap);
00127 
00128     //intial construction of matrix 
00129     Epetra_CrsMatrix * NewMatrix = new Epetra_CrsMatrix( View, *NewGraph_ );
00130 
00131     //insert views of row values
00132     int * myIndices;
00133     double * myValues;
00134     int indicesCnt;
00135     int numMyRows = NewMatrix->NumMyRows();
00136     for( int i = 0; i < numMyRows; ++i )
00137     {
00138       orig.ExtractMyRowView( i, indicesCnt, myValues, myIndices );
00139       NewGraph_->ExtractMyRowView( i, indicesCnt, myIndices );
00140 
00141       NewMatrix->InsertMyValues( i, indicesCnt, myValues, myIndices );
00142     }
00143 
00144     NewMatrix->FillComplete(DomainMap,RangeMap);
00145 
00146     newObj_ = NewMatrix;
00147   }
00148 
00149   return *newObj_;
00150 }
00151 
00152 } // namespace EpetraExt
00153 
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines