EpetraExt_BlockUtility.cpp

Go to the documentation of this file.
00001 //@HEADER
00002 // ************************************************************************
00003 // 
00004 //               EpetraExt: Extended Linear Algebra Services Package 
00005 //                 Copyright (2001) Sandia Corporation
00006 // 
00007 // Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
00008 // license for use of this work by or on behalf of the U.S. Government.
00009 // 
00010 // This library is free software; you can redistribute it and/or modify
00011 // it under the terms of the GNU Lesser General Public License as
00012 // published by the Free Software Foundation; either version 2.1 of the
00013 // License, or (at your option) any later version.
00014 //  
00015 // This library is distributed in the hope that it will be useful, but
00016 // WITHOUT ANY WARRANTY; without even the implied warranty of
00017 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
00018 // Lesser General Public License for more details.
00019 //  
00020 // You should have received a copy of the GNU Lesser General Public
00021 // License along with this library; if not, write to the Free Software
00022 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307
00023 // USA
00024 // Questions? Contact Michael A. Heroux (maherou@sandia.gov) 
00025 // 
00026 // ************************************************************************
00027 //@HEADER
00028 
00029 #include "EpetraExt_BlockUtility.h"
00030 #include "Epetra_Map.h"
00031 #include "Epetra_Comm.h"
00032 
00033 namespace EpetraExt {
00034 
00035 using std::vector;
00036 
00037 //==============================================================================
00038 Epetra_CrsGraph * BlockUtility::GenerateBlockGraph(
00039         const Epetra_CrsGraph & BaseGraph,
00040         const vector< vector<int> > & RowStencil,
00041         const vector<int> & RowIndices,
00042         const Epetra_Comm & GlobalComm ) 
00043 {
00044 
00045   const Epetra_BlockMap & BaseMap = BaseGraph.RowMap();
00046   int BaseIndex = BaseMap.IndexBase();
00047   int Offset = BlockUtility::CalculateOffset(BaseMap);
00048 
00049   //Get Base Global IDs
00050   int NumBlockRows = RowIndices.size();
00051   int Size = BaseMap.NumMyElements();
00052   int TotalSize = NumBlockRows * Size;
00053   vector<int> GIDs(Size);
00054   BaseMap.MyGlobalElements( &GIDs[0] );
00055 
00056   vector<int> GlobalGIDs( TotalSize );
00057   for( int i = 0; i < NumBlockRows; ++i )
00058   {
00059     for( int j = 0; j < Size; ++j )
00060       GlobalGIDs[i*Size+j] = GIDs[j] + RowIndices[i] * Offset;
00061   }
00062 
00063   int GlobalSize;
00064   GlobalComm.SumAll( &TotalSize, &GlobalSize, 1 );
00065 
00066   Epetra_Map GlobalMap( GlobalSize, TotalSize, &GlobalGIDs[0], BaseIndex, GlobalComm );
00067 
00068   int MaxIndices = BaseGraph.MaxNumIndices();
00069   vector<int> Indices(MaxIndices);
00070   int NumIndices;
00071 
00072   Epetra_CrsGraph * GlobalGraph = new Epetra_CrsGraph( Copy, 
00073                                dynamic_cast<Epetra_BlockMap&>(GlobalMap),
00074                                0 );
00075 
00076   for( int i = 0; i < NumBlockRows; ++i )
00077   {
00078     int StencilSize = RowStencil[i].size();
00079     for( int j = 0; j < Size; ++j )
00080     {
00081       int BaseRow = BaseMap.GID(j);
00082       int GlobalRow = GlobalMap.GID(j+i*Size);
00083 
00084       BaseGraph.ExtractGlobalRowCopy( BaseRow, MaxIndices, NumIndices, &Indices[0] );
00085       for( int k = 0; k < StencilSize; ++k )
00086       {
00087         int ColOffset = (RowIndices[i]+RowStencil[i][k]) * Offset;
00088         if( k > 0 ) ColOffset -= (RowIndices[i]+RowStencil[i][k-1]) * Offset;
00089 
00090         for( int l = 0; l < NumIndices; ++l )
00091           Indices[l] += ColOffset;
00092 
00093         GlobalGraph->InsertGlobalIndices( GlobalRow, NumIndices, &Indices[0] );
00094       }
00095     }
00096   }
00097 
00098   GlobalGraph->FillComplete();
00099 
00100   return GlobalGraph;
00101 }
00102 
00103 //==============================================================================
00104 Epetra_CrsGraph * BlockUtility::GenerateBlockGraph(
00105         const Epetra_RowMatrix & BaseMatrix,
00106         const vector< vector<int> > & RowStencil,
00107         const vector<int> & RowIndices,
00108         const Epetra_Comm & GlobalComm ) 
00109 {
00110 
00111   const Epetra_BlockMap & BaseMap = BaseMatrix.RowMatrixRowMap();
00112   const Epetra_BlockMap & BaseColMap = BaseMatrix.RowMatrixColMap();
00113   int BaseIndex = BaseMap.IndexBase();
00114   int Offset = BlockUtility::CalculateOffset(BaseMap);
00115 
00116   //Get Base Global IDs
00117   int NumBlockRows = RowIndices.size();
00118   int Size = BaseMap.NumMyElements();
00119   int TotalSize = NumBlockRows * Size;
00120   vector<int> GIDs(Size);
00121   BaseMap.MyGlobalElements( &GIDs[0] );
00122 
00123   vector<int> GlobalGIDs( TotalSize );
00124   for( int i = 0; i < NumBlockRows; ++i )
00125   {
00126     for( int j = 0; j < Size; ++j )
00127       GlobalGIDs[i*Size+j] = GIDs[j] + RowIndices[i] * Offset;
00128   }
00129 
00130   int GlobalSize;
00131   GlobalComm.SumAll( &TotalSize, &GlobalSize, 1 );
00132 
00133   Epetra_Map GlobalMap( GlobalSize, TotalSize, &GlobalGIDs[0], BaseIndex, GlobalComm );
00134 
00135   int MaxIndices = BaseMatrix.MaxNumEntries();
00136   vector<int> Indices(MaxIndices);
00137   vector<double> Values(MaxIndices); 
00138   int NumIndices;
00139 
00140   Epetra_CrsGraph * GlobalGraph = new Epetra_CrsGraph( Copy, 
00141                                dynamic_cast<Epetra_BlockMap&>(GlobalMap),
00142                                0 );
00143 
00144   for( int i = 0; i < NumBlockRows; ++i )
00145   {
00146     int StencilSize = RowStencil[i].size();
00147     for( int j = 0; j < Size; ++j )
00148     {
00149       int GlobalRow = GlobalMap.GID(j+i*Size);
00150 
00151       BaseMatrix.ExtractMyRowCopy( j, MaxIndices, NumIndices, &Values[0], &Indices[0] );
00152       for( int l = 0; l < NumIndices; ++l ) Indices[l] = BaseColMap.GID(Indices[l]);
00153 
00154       for( int k = 0; k < StencilSize; ++k )
00155       {
00156         int ColOffset = (RowIndices[i]+RowStencil[i][k]) * Offset;
00157         if( k > 0 ) ColOffset -= (RowIndices[i]+RowStencil[i][k-1]) * Offset;
00158 
00159         for( int l = 0; l < NumIndices; ++l )
00160           Indices[l] += ColOffset;
00161 
00162         GlobalGraph->InsertGlobalIndices( GlobalRow, NumIndices, &Indices[0] );
00163       }
00164     }
00165   }
00166 
00167   GlobalGraph->FillComplete();
00168 
00169   return GlobalGraph;
00170 }
00171 
00172 //==============================================================================
00173 int BlockUtility::CalculateOffset(const Epetra_BlockMap & BaseMap)
00174 {
00175   int MaxGID = BaseMap.MaxAllGID();
00176 
00177   int Offset = 1;
00178   while( Offset < MaxGID ) Offset *= 10;
00179 
00180   return Offset;
00181 }
00182 
00183 
00184 } //namespace EpetraExt

Generated on Thu Sep 18 12:31:43 2008 for EpetraExt by doxygen 1.3.9.1