EpetraExt Package Browser (Single Doxygen Collection) Development
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_Map * BlockUtility::GenerateBlockMap(
00039   const Epetra_BlockMap & BaseMap,
00040         const vector<int> & RowIndices,
00041         const Epetra_Comm & GlobalComm ) 
00042 {
00043   int BaseIndex = BaseMap.IndexBase();
00044   int Offset = BlockUtility::CalculateOffset(BaseMap);
00045 
00046   //Get Base Global IDs
00047   int NumBlockRows = RowIndices.size();
00048   int Size = BaseMap.NumMyElements();
00049   int TotalSize = NumBlockRows * Size;
00050   vector<int> GIDs(Size);
00051   BaseMap.MyGlobalElements( &GIDs[0] );
00052 
00053   vector<int> GlobalGIDs( TotalSize );
00054   for( int i = 0; i < NumBlockRows; ++i )
00055   {
00056     for( int j = 0; j < Size; ++j )
00057       GlobalGIDs[i*Size+j] = GIDs[j] + RowIndices[i] * Offset;
00058   }
00059 
00060   int GlobalSize;
00061   GlobalComm.SumAll( &TotalSize, &GlobalSize, 1 );
00062 
00063   Epetra_Map *GlobalMap = 
00064     new Epetra_Map( GlobalSize, TotalSize, &GlobalGIDs[0], BaseIndex, 
00065         GlobalComm );
00066 
00067   return GlobalMap;
00068 }
00069 
00070 //==============================================================================
00071 Epetra_CrsGraph * BlockUtility::GenerateBlockGraph(
00072         const Epetra_CrsGraph & BaseGraph,
00073         const vector< vector<int> > & RowStencil,
00074         const vector<int> & RowIndices,
00075         const Epetra_Comm & GlobalComm ) 
00076 {
00077 
00078   const Epetra_BlockMap & BaseMap = BaseGraph.RowMap();
00079   int BaseIndex = BaseMap.IndexBase();
00080   int Offset = BlockUtility::CalculateOffset(BaseMap);
00081 
00082   //Get Base Global IDs
00083   int NumBlockRows = RowIndices.size();
00084   int Size = BaseMap.NumMyElements();
00085   int TotalSize = NumBlockRows * Size;
00086   vector<int> GIDs(Size);
00087   BaseMap.MyGlobalElements( &GIDs[0] );
00088 
00089   vector<int> GlobalGIDs( TotalSize );
00090   for( int i = 0; i < NumBlockRows; ++i )
00091   {
00092     for( int j = 0; j < Size; ++j )
00093       GlobalGIDs[i*Size+j] = GIDs[j] + RowIndices[i] * Offset;
00094   }
00095 
00096   int GlobalSize;
00097   GlobalComm.SumAll( &TotalSize, &GlobalSize, 1 );
00098 
00099   Epetra_Map GlobalMap( GlobalSize, TotalSize, &GlobalGIDs[0], BaseIndex, GlobalComm );
00100 
00101   int MaxIndices = BaseGraph.MaxNumIndices();
00102   vector<int> Indices(MaxIndices);
00103   int NumIndices;
00104 
00105   Epetra_CrsGraph * GlobalGraph = new Epetra_CrsGraph( Copy, 
00106                                dynamic_cast<Epetra_BlockMap&>(GlobalMap),
00107                                0 );
00108 
00109   for( int i = 0; i < NumBlockRows; ++i )
00110   {
00111     int StencilSize = RowStencil[i].size();
00112     for( int j = 0; j < Size; ++j )
00113     {
00114       int BaseRow = BaseMap.GID(j);
00115       int GlobalRow = GlobalMap.GID(j+i*Size);
00116 
00117       BaseGraph.ExtractGlobalRowCopy( BaseRow, MaxIndices, NumIndices, &Indices[0] );
00118       for( int k = 0; k < StencilSize; ++k )
00119       {
00120         int ColOffset = (RowIndices[i]+RowStencil[i][k]) * Offset;
00121         if( k > 0 ) ColOffset -= (RowIndices[i]+RowStencil[i][k-1]) * Offset;
00122 
00123         for( int l = 0; l < NumIndices; ++l )
00124           Indices[l] += ColOffset;
00125 
00126         GlobalGraph->InsertGlobalIndices( GlobalRow, NumIndices, &Indices[0] );
00127       }
00128     }
00129   }
00130 
00131   GlobalGraph->FillComplete();
00132 
00133   return GlobalGraph;
00134 }
00135 
00136 //==============================================================================
00137 Epetra_CrsGraph * BlockUtility::GenerateBlockGraph(
00138         const Epetra_RowMatrix & BaseMatrix,
00139         const vector< vector<int> > & RowStencil,
00140         const vector<int> & RowIndices,
00141         const Epetra_Comm & GlobalComm ) 
00142 {
00143 
00144   const Epetra_BlockMap & BaseMap = BaseMatrix.RowMatrixRowMap();
00145   const Epetra_BlockMap & BaseColMap = BaseMatrix.RowMatrixColMap();
00146   int BaseIndex = BaseMap.IndexBase();
00147   int Offset = BlockUtility::CalculateOffset(BaseMap);
00148 
00149   //Get Base Global IDs
00150   int NumBlockRows = RowIndices.size();
00151   int Size = BaseMap.NumMyElements();
00152   int TotalSize = NumBlockRows * Size;
00153   vector<int> GIDs(Size);
00154   BaseMap.MyGlobalElements( &GIDs[0] );
00155 
00156   vector<int> GlobalGIDs( TotalSize );
00157   for( int i = 0; i < NumBlockRows; ++i )
00158   {
00159     for( int j = 0; j < Size; ++j )
00160       GlobalGIDs[i*Size+j] = GIDs[j] + RowIndices[i] * Offset;
00161   }
00162 
00163   int GlobalSize;
00164   GlobalComm.SumAll( &TotalSize, &GlobalSize, 1 );
00165 
00166   Epetra_Map GlobalMap( GlobalSize, TotalSize, &GlobalGIDs[0], BaseIndex, GlobalComm );
00167 
00168   int MaxIndices = BaseMatrix.MaxNumEntries();
00169   vector<int> Indices(MaxIndices);
00170   vector<double> Values(MaxIndices); 
00171   int NumIndices;
00172 
00173   Epetra_CrsGraph * GlobalGraph = new Epetra_CrsGraph( Copy, 
00174                                dynamic_cast<Epetra_BlockMap&>(GlobalMap),
00175                                0 );
00176 
00177   for( int i = 0; i < NumBlockRows; ++i )
00178   {
00179     int StencilSize = RowStencil[i].size();
00180     for( int j = 0; j < Size; ++j )
00181     {
00182       int GlobalRow = GlobalMap.GID(j+i*Size);
00183 
00184       BaseMatrix.ExtractMyRowCopy( j, MaxIndices, NumIndices, &Values[0], &Indices[0] );
00185       for( int l = 0; l < NumIndices; ++l ) Indices[l] = BaseColMap.GID(Indices[l]);
00186 
00187       for( int k = 0; k < StencilSize; ++k )
00188       {
00189         int ColOffset = (RowIndices[i]+RowStencil[i][k]) * Offset;
00190         if( k > 0 ) ColOffset -= (RowIndices[i]+RowStencil[i][k-1]) * Offset;
00191 
00192         for( int l = 0; l < NumIndices; ++l )
00193           Indices[l] += ColOffset;
00194 
00195         GlobalGraph->InsertGlobalIndices( GlobalRow, NumIndices, &Indices[0] );
00196       }
00197     }
00198   }
00199 
00200   GlobalGraph->FillComplete();
00201 
00202   return GlobalGraph;
00203 }
00204 
00205 //==============================================================================
00206 int BlockUtility::CalculateOffset(const Epetra_BlockMap & BaseMap)
00207 {
00208   int MaxGID = BaseMap.MaxAllGID();
00209 
00210 //   int Offset = 1;
00211 //   while( Offset <= MaxGID ) Offset *= 10;
00212 
00213 //   return Offset;
00214 
00215   return MaxGID+1;
00216 }
00217 
00218 
00219 } //namespace EpetraExt
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines