EpetraExt Development
EpetraExt_BlockCrsMatrix.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_BlockCrsMatrix.h"
00043 #include "EpetraExt_BlockUtility.h"
00044 #include "Epetra_Map.h"
00045 
00046 namespace EpetraExt {
00047 
00048 using std::vector;
00049 
00050 //==============================================================================
00051 #ifndef EPETRA_NO_32BIT_GLOBAL_INDICES
00052 BlockCrsMatrix::BlockCrsMatrix(
00053         const Epetra_CrsGraph & BaseGraph,
00054         const vector<int> & RowStencil,
00055         int RowIndex,
00056         const Epetra_Comm & GlobalComm  ) 
00057   : Epetra_CrsMatrix( Copy, *(BlockUtility::GenerateBlockGraph( BaseGraph, vector< vector<int> >(1,RowStencil), vector<int>(1,RowIndex), GlobalComm )) ),
00058     BaseGraph_( BaseGraph ),
00059     RowStencil_int_( vector< vector<int> >(1,RowStencil) ),
00060     RowIndices_int_( vector<int>(1,RowIndex) ),
00061     ROffset_(BlockUtility::CalculateOffset64(BaseGraph.RowMap())),
00062     COffset_(BlockUtility::CalculateOffset64(BaseGraph.ColMap()))
00063 {
00064 }
00065 #endif
00066 
00067 #ifndef EPETRA_NO_64BIT_GLOBAL_INDICES
00068 BlockCrsMatrix::BlockCrsMatrix(
00069         const Epetra_CrsGraph & BaseGraph,
00070         const vector<long long> & RowStencil,
00071         long long RowIndex,
00072         const Epetra_Comm & GlobalComm  ) 
00073   : Epetra_CrsMatrix( Copy, *(BlockUtility::GenerateBlockGraph( BaseGraph, vector< vector<long long> >(1,RowStencil), vector<long long>(1,RowIndex), GlobalComm )) ),
00074     BaseGraph_( BaseGraph ),
00075     RowStencil_LL_( vector< vector<long long> >(1,RowStencil) ),
00076     RowIndices_LL_( vector<long long>(1,RowIndex) ),
00077     ROffset_(BlockUtility::CalculateOffset64(BaseGraph.RowMap())),
00078     COffset_(BlockUtility::CalculateOffset64(BaseGraph.ColMap()))
00079 {
00080 }
00081 #endif
00082 
00083 //==============================================================================
00084 #ifndef EPETRA_NO_32BIT_GLOBAL_INDICES
00085 BlockCrsMatrix::BlockCrsMatrix(
00086         const Epetra_CrsGraph & BaseGraph,
00087         const vector< vector<int> > & RowStencil,
00088         const vector<int> & RowIndices,
00089         const Epetra_Comm & GlobalComm  ) 
00090   : Epetra_CrsMatrix( Copy, *(BlockUtility::GenerateBlockGraph( BaseGraph, RowStencil, RowIndices, GlobalComm )) ),
00091     BaseGraph_( BaseGraph ),
00092     RowStencil_int_( RowStencil ),
00093     RowIndices_int_( RowIndices ),
00094     ROffset_(BlockUtility::CalculateOffset64(BaseGraph.RowMap())),
00095     COffset_(BlockUtility::CalculateOffset64(BaseGraph.ColMap()))
00096 {
00097 }
00098 #endif
00099 #ifndef EPETRA_NO_64BIT_GLOBAL_INDICES
00100 BlockCrsMatrix::BlockCrsMatrix(
00101         const Epetra_CrsGraph & BaseGraph,
00102         const vector< vector<long long> > & RowStencil,
00103         const vector<long long> & RowIndices,
00104         const Epetra_Comm & GlobalComm  ) 
00105   : Epetra_CrsMatrix( Copy, *(BlockUtility::GenerateBlockGraph( BaseGraph, RowStencil, RowIndices, GlobalComm )) ),
00106     BaseGraph_( BaseGraph ),
00107     RowStencil_LL_( RowStencil ),
00108     RowIndices_LL_( RowIndices ),
00109     ROffset_(BlockUtility::CalculateOffset64(BaseGraph.RowMap())),
00110     COffset_(BlockUtility::CalculateOffset64(BaseGraph.ColMap()))
00111 {
00112 }
00113 #endif
00114 
00115 //==============================================================================
00116 BlockCrsMatrix::BlockCrsMatrix(
00117         const Epetra_CrsGraph & BaseGraph,
00118         const Epetra_CrsGraph & LocalBlockGraph,
00119         const Epetra_Comm & GlobalComm  ) 
00120   : Epetra_CrsMatrix( Copy, *(BlockUtility::GenerateBlockGraph( BaseGraph, LocalBlockGraph, GlobalComm )) ),
00121     BaseGraph_( BaseGraph ),
00122 #ifndef EPETRA_NO_32BIT_GLOBAL_INDICES
00123     RowStencil_int_( ),
00124     RowIndices_int_( ),
00125 #endif
00126 #ifndef EPETRA_NO_64BIT_GLOBAL_INDICES
00127     RowStencil_LL_( ),
00128     RowIndices_LL_( ),
00129 #endif
00130     ROffset_(BlockUtility::CalculateOffset64(BaseGraph.RowMap())),
00131     COffset_(BlockUtility::CalculateOffset64(BaseGraph.ColMap()))
00132 {
00133 #ifndef EPETRA_NO_32BIT_GLOBAL_INDICES
00134   if(Epetra_CrsMatrix::RowMatrixRowMap().GlobalIndicesInt() && LocalBlockGraph.RowMap().GlobalIndicesInt())
00135     BlockUtility::GenerateRowStencil(LocalBlockGraph, RowIndices_int_, RowStencil_int_);
00136   else
00137 #endif
00138 #ifndef EPETRA_NO_64BIT_GLOBAL_INDICES
00139   if(Epetra_CrsMatrix::RowMatrixRowMap().GlobalIndicesLongLong() && LocalBlockGraph.RowMap().GlobalIndicesLongLong())
00140     BlockUtility::GenerateRowStencil(LocalBlockGraph, RowIndices_LL_, RowStencil_LL_);
00141   else
00142 #endif
00143     throw "EpetraExt::BlockCrsMatrix::BlockCrsMatrix: Error, Global indices unknown.";
00144 }
00145 
00146 //==============================================================================
00147 #ifndef EPETRA_NO_32BIT_GLOBAL_INDICES
00148 BlockCrsMatrix::BlockCrsMatrix(
00149         const Epetra_RowMatrix & BaseMatrix,
00150         const vector< vector<int> > & RowStencil,
00151         const vector<int> & RowIndices,
00152         const Epetra_Comm & GlobalComm  ) 
00153   : Epetra_CrsMatrix( Copy, *(BlockUtility::GenerateBlockGraph( BaseMatrix, RowStencil, RowIndices, GlobalComm )) ),
00154     BaseGraph_( Copy, BaseMatrix.RowMatrixRowMap(), 1 ), //Junk to satisfy constructor
00155     RowStencil_int_( RowStencil ),
00156     RowIndices_int_( RowIndices ),
00157     ROffset_(BlockUtility::CalculateOffset64(BaseMatrix.RowMatrixRowMap())),
00158     COffset_(BlockUtility::CalculateOffset64(BaseMatrix.RowMatrixColMap()))
00159 {
00160 }
00161 #endif
00162 
00163 #ifndef EPETRA_NO_64BIT_GLOBAL_INDICES
00164 BlockCrsMatrix::BlockCrsMatrix(
00165         const Epetra_RowMatrix & BaseMatrix,
00166         const vector< vector<long long> > & RowStencil,
00167         const vector<long long> & RowIndices,
00168         const Epetra_Comm & GlobalComm  ) 
00169   : Epetra_CrsMatrix( Copy, *(BlockUtility::GenerateBlockGraph( BaseMatrix, RowStencil, RowIndices, GlobalComm )) ),
00170     BaseGraph_( Copy, BaseMatrix.RowMatrixRowMap(), 1 ), //Junk to satisfy constructor
00171     RowStencil_LL_( RowStencil ),
00172     RowIndices_LL_( RowIndices ),
00173     ROffset_(BlockUtility::CalculateOffset64(BaseMatrix.RowMatrixRowMap())),
00174     COffset_(BlockUtility::CalculateOffset64(BaseMatrix.RowMatrixColMap()))
00175 {
00176 }
00177 #endif
00178 
00179 //==============================================================================
00180 BlockCrsMatrix::BlockCrsMatrix( const BlockCrsMatrix & Matrix ) 
00181   : Epetra_CrsMatrix( dynamic_cast<const Epetra_CrsMatrix &>( Matrix ) ),
00182     BaseGraph_( Matrix.BaseGraph_ ),
00183 #ifndef EPETRA_NO_32BIT_GLOBAL_INDICES
00184     RowStencil_int_( Matrix.RowStencil_int_ ),
00185     RowIndices_int_( Matrix.RowIndices_int_ ),
00186 #endif
00187 #ifndef EPETRA_NO_64BIT_GLOBAL_INDICES
00188     RowStencil_LL_( Matrix.RowStencil_LL_ ),
00189     RowIndices_LL_( Matrix.RowIndices_LL_ ),
00190 #endif
00191     ROffset_( Matrix.ROffset_ ),
00192     COffset_( Matrix.COffset_ )
00193 {
00194 }
00195 
00196 //==============================================================================
00197 BlockCrsMatrix::~BlockCrsMatrix()
00198 {
00199 }
00200 
00201 //==============================================================================
00202 template<typename int_type>
00203 void BlockCrsMatrix::TLoadBlock(const Epetra_RowMatrix & BaseMatrix, const int_type Row, const int_type Col)
00204 {
00205   std::vector<int_type>& RowIndices_ = TRowIndices<int_type>();
00206   std::vector< std::vector<int_type> >& RowStencil_ = TRowStencil<int_type>();
00207   int_type RowOffset = RowIndices_[(std::size_t)Row] * ROffset_;
00208   int_type ColOffset = (RowIndices_[(std::size_t)Row] + RowStencil_[(std::size_t)Row][(std::size_t)Col]) * COffset_;
00209 
00210 //  const Epetra_CrsGraph & BaseGraph = BaseMatrix.Graph();
00211   const Epetra_BlockMap & BaseMap = BaseMatrix.RowMatrixRowMap();
00212   const Epetra_BlockMap & BaseColMap = BaseMatrix.RowMatrixColMap();
00213 
00214   // This routine copies entries of a BaseMatrix into big  BlockCrsMatrix
00215   // It performs the following operation on the global IDs row-by-row
00216   // this->val[i+rowOffset][j+ColOffset] = BaseMatrix.val[i][j]
00217 
00218   int MaxIndices = BaseMatrix.MaxNumEntries();
00219   vector<int> Indices_local(MaxIndices);
00220   vector<int_type> Indices_global(MaxIndices);
00221   vector<double> Values(MaxIndices);
00222   int NumIndices;
00223   int ierr=0;
00224 
00225   for (int i=0; i<BaseMap.NumMyElements(); i++) {
00226     BaseMatrix.ExtractMyRowCopy( i, MaxIndices, NumIndices, &Values[0], &Indices_local[0] );
00227 
00228     // Convert to BlockMatrix Global numbering scheme
00229     for( int l = 0; l < NumIndices; ++l )
00230        Indices_global[l] = ColOffset +  (int_type) BaseColMap.GID64(Indices_local[l]);
00231 
00232     int_type BaseRow = (int_type) BaseMap.GID64(i);
00233     ierr = this->ReplaceGlobalValues(BaseRow + RowOffset, NumIndices, &Values[0], &Indices_global[0]); 
00234     if (ierr != 0) std::cout << "WARNING BlockCrsMatrix::LoadBlock ReplaceGlobalValues err = " << ierr <<
00235       "\n\t  Row " << BaseRow + RowOffset << "Col start" << Indices_global[0] << std::endl;
00236 
00237   }
00238 }
00239 
00240 #ifndef EPETRA_NO_32BIT_GLOBAL_INDICES
00241 void BlockCrsMatrix::LoadBlock(const Epetra_RowMatrix & BaseMatrix, const int Row, const int Col)
00242 {
00243   if(Epetra_CrsMatrix::RowMatrixRowMap().GlobalIndicesInt() && BaseMatrix.RowMatrixRowMap().GlobalIndicesInt())
00244   return TLoadBlock<int>(BaseMatrix, Row, Col);
00245   else
00246     throw "EpetraExt::BlockCrsMatrix::LoadBlock: Global indices not int";
00247 }
00248 #endif
00249 
00250 #ifndef EPETRA_NO_64BIT_GLOBAL_INDICES
00251 void BlockCrsMatrix::LoadBlock(const Epetra_RowMatrix & BaseMatrix, const long long Row, const long long Col)
00252 {
00253   if(Epetra_CrsMatrix::RowMatrixRowMap().GlobalIndicesLongLong() && BaseMatrix.RowMatrixRowMap().GlobalIndicesLongLong())
00254   return TLoadBlock<long long>(BaseMatrix, Row, Col);
00255   else
00256     throw "EpetraExt::BlockCrsMatrix::LoadBlock: Global indices not long long";
00257 }
00258 #endif
00259 
00260 //==============================================================================
00261 template<typename int_type>
00262 void BlockCrsMatrix::TSumIntoBlock(double alpha, const Epetra_RowMatrix & BaseMatrix, const int_type Row, const int_type Col)
00263 {
00264   std::vector<int_type>& RowIndices_ = TRowIndices<int_type>();
00265   std::vector< std::vector<int_type> >& RowStencil_ = TRowStencil<int_type>();
00266   int_type RowOffset = RowIndices_[(std::size_t)Row] * ROffset_;
00267   int_type ColOffset = (RowIndices_[(std::size_t)Row] + RowStencil_[(std::size_t)Row][(std::size_t)Col]) * COffset_;
00268 
00269 //  const Epetra_CrsGraph & BaseGraph = BaseMatrix.Graph();
00270   const Epetra_BlockMap & BaseMap = BaseMatrix.RowMatrixRowMap();
00271   const Epetra_BlockMap & BaseColMap = BaseMatrix.RowMatrixColMap();
00272 
00273   // This routine copies entries of a BaseMatrix into big  BlockCrsMatrix
00274   // It performs the following operation on the global IDs row-by-row
00275   // this->val[i+rowOffset][j+ColOffset] = BaseMatrix.val[i][j]
00276 
00277   int MaxIndices = BaseMatrix.MaxNumEntries();
00278   vector<int> Indices_local(MaxIndices);
00279   vector<int_type> Indices_global(MaxIndices);
00280   vector<double> Values(MaxIndices);
00281   int NumIndices;
00282   int ierr=0;
00283 
00284   for (int i=0; i<BaseMap.NumMyElements(); i++) {
00285     BaseMatrix.ExtractMyRowCopy( i, MaxIndices, NumIndices, &Values[0], &Indices_local[0] );
00286 
00287     // Convert to BlockMatrix Global numbering scheme
00288     for( int l = 0; l < NumIndices; ++l ) {
00289        Indices_global[l] = ColOffset +  (int_type) BaseColMap.GID64(Indices_local[l]);
00290        Values[l] *= alpha;
00291     }
00292 
00293     int_type BaseRow = (int_type) BaseMap.GID64(i);
00294     ierr = this->SumIntoGlobalValues(BaseRow + RowOffset, NumIndices, &Values[0], &Indices_global[0]); 
00295     if (ierr != 0) std::cout << "WARNING BlockCrsMatrix::SumIntoBlock SumIntoGlobalValues err = " << ierr <<
00296       "\n\t  Row " << BaseRow + RowOffset << "Col start" << Indices_global[0] << std::endl;
00297 
00298   }
00299 }
00300 
00301 #ifndef EPETRA_NO_32BIT_GLOBAL_INDICES
00302 void BlockCrsMatrix::SumIntoBlock(double alpha, const Epetra_RowMatrix & BaseMatrix, const int Row, const int Col)
00303 {
00304   if(Epetra_CrsMatrix::RowMatrixRowMap().GlobalIndicesInt() && BaseMatrix.RowMatrixRowMap().GlobalIndicesInt())
00305   return TSumIntoBlock<int>(alpha, BaseMatrix, Row, Col);
00306   else
00307     throw "EpetraExt::BlockCrsMatrix::SumIntoBlock: Global indices not int";
00308 }
00309 #endif
00310 
00311 #ifndef EPETRA_NO_64BIT_GLOBAL_INDICES
00312 void BlockCrsMatrix::SumIntoBlock(double alpha, const Epetra_RowMatrix & BaseMatrix, const long long Row, const long long Col)
00313 {
00314   if(Epetra_CrsMatrix::RowMatrixRowMap().GlobalIndicesLongLong() && BaseMatrix.RowMatrixRowMap().GlobalIndicesLongLong())
00315   return TSumIntoBlock<long long>(alpha, BaseMatrix, Row, Col);
00316   else
00317     throw "EpetraExt::BlockCrsMatrix::SumIntoBlock: Global indices not long long";
00318 }
00319 #endif
00320 
00321 //==============================================================================
00322 template<typename int_type>
00323 void BlockCrsMatrix::TSumIntoGlobalBlock(double alpha, const Epetra_RowMatrix & BaseMatrix, const int_type Row, const int_type Col)
00324 {
00325   int_type RowOffset = Row * ROffset_;
00326   int_type ColOffset = Col * COffset_;
00327 
00328 //  const Epetra_CrsGraph & BaseGraph = BaseMatrix.Graph();
00329   const Epetra_BlockMap & BaseMap = BaseMatrix.RowMatrixRowMap();
00330   const Epetra_BlockMap & BaseColMap = BaseMatrix.RowMatrixColMap();
00331 
00332   // This routine copies entries of a BaseMatrix into big  BlockCrsMatrix
00333   // It performs the following operation on the global IDs row-by-row
00334   // this->val[i+rowOffset][j+ColOffset] = BaseMatrix.val[i][j]
00335 
00336   int MaxIndices = BaseMatrix.MaxNumEntries();
00337   vector<int> Indices_local(MaxIndices);
00338   vector<int_type> Indices_global(MaxIndices);
00339   vector<double> Values(MaxIndices);
00340   int NumIndices;
00341   int ierr=0;
00342 
00343   for (int i=0; i<BaseMap.NumMyElements(); i++) {
00344     BaseMatrix.ExtractMyRowCopy( i, MaxIndices, NumIndices, &Values[0], &Indices_local[0] );
00345 
00346     // Convert to BlockMatrix Global numbering scheme
00347     for( int l = 0; l < NumIndices; ++l ) {
00348        Indices_global[l] = ColOffset +  (int_type) BaseColMap.GID64(Indices_local[l]);
00349        Values[l] *= alpha;
00350     }
00351 
00352     int_type BaseRow = (int_type) BaseMap.GID64(i);
00353     ierr = this->SumIntoGlobalValues(BaseRow + RowOffset, NumIndices, &Values[0], &Indices_global[0]); 
00354     if (ierr != 0) std::cout << "WARNING BlockCrsMatrix::SumIntoBlock SumIntoGlobalValues err = " << ierr <<
00355       "\n\t  Row " << BaseRow + RowOffset << "Col start" << Indices_global[0] << std::endl;
00356 
00357   }
00358 }
00359 
00360 #ifndef EPETRA_NO_32BIT_GLOBAL_INDICES
00361 void BlockCrsMatrix::SumIntoGlobalBlock(double alpha, const Epetra_RowMatrix & BaseMatrix, const int Row, const int Col)
00362 {
00363   if(Epetra_CrsMatrix::RowMatrixRowMap().GlobalIndicesInt() && BaseMatrix.RowMatrixRowMap().GlobalIndicesInt())
00364   return TSumIntoGlobalBlock<int>(alpha, BaseMatrix, Row, Col);
00365   else
00366     throw "EpetraExt::BlockCrsMatrix::SumIntoGlobalBlock: Global indices not int";
00367 }
00368 #endif
00369 
00370 #ifndef EPETRA_NO_64BIT_GLOBAL_INDICES
00371 void BlockCrsMatrix::SumIntoGlobalBlock(double alpha, const Epetra_RowMatrix & BaseMatrix, const long long Row, const long long Col)
00372 {
00373   if(Epetra_CrsMatrix::RowMatrixRowMap().GlobalIndicesLongLong() && BaseMatrix.RowMatrixRowMap().GlobalIndicesLongLong())
00374   return TSumIntoGlobalBlock<long long>(alpha, BaseMatrix, Row, Col);
00375   else
00376     throw "EpetraExt::BlockCrsMatrix::SumIntoGlobalBlock: Global indices not long long";
00377 }
00378 #endif
00379 
00380 
00381 //==============================================================================
00382 template<typename int_type>
00383 void BlockCrsMatrix::TBlockSumIntoGlobalValues(const int_type BaseRow, int NumIndices,
00384      double* Values, const int_type* Indices, const int_type Row, const int_type Col)
00385 //All arguments could be const, except some were not set as const in CrsMatrix
00386 {
00387   std::vector<int_type>& RowIndices_ = TRowIndices<int_type>();
00388   std::vector< std::vector<int_type> >& RowStencil_ = TRowStencil<int_type>();
00389   int_type RowOffset = RowIndices_[(std::size_t)Row] * ROffset_;
00390   int_type ColOffset = (RowIndices_[(std::size_t)Row] + RowStencil_[(std::size_t)Row][(std::size_t)Col]) * COffset_;
00391 
00392   // Convert to BlockMatrix Global numbering scheme
00393   vector<int_type> OffsetIndices(NumIndices);
00394   for( int l = 0; l < NumIndices; ++l ) OffsetIndices[l] = Indices[l] + ColOffset;
00395 
00396   int ierr = this->SumIntoGlobalValues(BaseRow + RowOffset, NumIndices,
00397                                    Values, &OffsetIndices[0]); 
00398 
00399   if (ierr != 0) std::cout << "WARNING BlockCrsMatrix::BlockSumIntoGlobalValues err = "
00400      << ierr << "\n\t  Row " << BaseRow + RowOffset << "Col start" << Indices[0] << std::endl;
00401 }
00402 
00403 #ifndef EPETRA_NO_32BIT_GLOBAL_INDICES
00404 void BlockCrsMatrix::BlockSumIntoGlobalValues(const int BaseRow, int NumIndices,
00405      double* Values, const int* Indices, const int Row, const int Col)
00406 {
00407   if(Epetra_CrsMatrix::RowMatrixRowMap().GlobalIndicesInt())
00408   return TBlockSumIntoGlobalValues<int>(BaseRow, NumIndices, Values, Indices, Row, Col);
00409   else
00410     throw "EpetraExt::BlockCrsMatrix::BlockSumIntoGlobalValues: Global indices not int";
00411 }
00412 #endif
00413 
00414 #ifndef EPETRA_NO_64BIT_GLOBAL_INDICES
00415 void BlockCrsMatrix::BlockSumIntoGlobalValues(const long long BaseRow, int NumIndices,
00416      double* Values, const long long* Indices, const long long Row, const long long Col)
00417 {
00418   if(Epetra_CrsMatrix::RowMatrixRowMap().GlobalIndicesLongLong() && Epetra_CrsMatrix::RowMatrixRowMap().GlobalIndicesLongLong())
00419   return TBlockSumIntoGlobalValues<long long>(BaseRow, NumIndices, Values, Indices, Row, Col);
00420   else
00421     throw "EpetraExt::BlockCrsMatrix::BlockSumIntoGlobalValues: Global indices not long long";
00422 }
00423 #endif
00424 
00425 //==============================================================================
00426 template<typename int_type>
00427 void BlockCrsMatrix::TBlockReplaceGlobalValues(const int_type BaseRow, int NumIndices,
00428      double* Values, const int_type* Indices, const int_type Row, const int_type Col)
00429 //All arguments could be const, except some were not set as const in CrsMatrix
00430 {
00431   std::vector<int_type>& RowIndices_ = TRowIndices<int_type>();
00432   std::vector< std::vector<int_type> >& RowStencil_ = TRowStencil<int_type>();
00433   int_type RowOffset = RowIndices_[(std::size_t)Row] * ROffset_;
00434   int_type ColOffset = (RowIndices_[(std::size_t)Row] + RowStencil_[(std::size_t)Row][(std::size_t)Col]) * COffset_;
00435 
00436   // Convert to BlockMatrix Global numbering scheme
00437   vector<int_type> OffsetIndices(NumIndices);
00438   for( int l = 0; l < NumIndices; ++l ) OffsetIndices[l] = Indices[l] + ColOffset;
00439 
00440   int ierr = this->ReplaceGlobalValues(BaseRow + RowOffset, NumIndices,
00441                                    Values, &OffsetIndices[0]); 
00442 
00443   if (ierr != 0) std::cout << "WARNING BlockCrsMatrix::BlockReplaceGlobalValues err = "
00444      << ierr << "\n\t  Row " << BaseRow + RowOffset << "Col start" << Indices[0] << std::endl;
00445 }
00446 
00447 #ifndef EPETRA_NO_32BIT_GLOBAL_INDICES
00448 void BlockCrsMatrix::BlockReplaceGlobalValues(const int BaseRow, int NumIndices,
00449      double* Values, const int* Indices, const int Row, const int Col)
00450 {
00451   if(Epetra_CrsMatrix::RowMatrixRowMap().GlobalIndicesInt())
00452   return TBlockReplaceGlobalValues<int>(BaseRow, NumIndices, Values, Indices, Row, Col);
00453   else
00454     throw "EpetraExt::BlockCrsMatrix::BlockReplaceGlobalValues: Global indices not int";
00455 }
00456 #endif
00457 
00458 #ifndef EPETRA_NO_64BIT_GLOBAL_INDICES
00459 void BlockCrsMatrix::BlockReplaceGlobalValues(const long long BaseRow, int NumIndices,
00460      double* Values, const long long* Indices, const long long Row, const long long Col)
00461 {
00462   if(Epetra_CrsMatrix::RowMatrixRowMap().GlobalIndicesLongLong())
00463   return TBlockReplaceGlobalValues<long long>(BaseRow, NumIndices, Values, Indices, Row, Col);
00464   else
00465     throw "EpetraExt::BlockCrsMatrix::BlockReplaceGlobalValues: Global indices not long long";
00466 }
00467 #endif
00468 
00469 //==============================================================================
00470 template<typename int_type>
00471 void BlockCrsMatrix::TBlockExtractGlobalRowView(const int_type BaseRow, 
00472                  int& NumEntries, 
00473                  double*& Values, 
00474                  const int_type Row, 
00475                  const int_type Col)
00476 //All arguments could be const, except some were not set as const in CrsMatrix
00477 {
00478   std::vector<int_type>& RowIndices_ = TRowIndices<int_type>();
00479   std::vector< std::vector<int_type> >& RowStencil_ = TRowStencil<int_type>();
00480   int_type RowOffset = RowIndices_[(std::size_t)Row] * ROffset_;
00481   int_type ColOffset = (RowIndices_[(std::size_t)Row] + RowStencil_[(std::size_t)Row][(std::size_t)Col]) * COffset_;
00482 
00483   // Get the whole row
00484   int ierr = this->ExtractGlobalRowView(BaseRow + RowOffset, NumEntries,
00485           Values); 
00486 
00487   // Adjust for just this block column
00488   Values += ColOffset;
00489   NumEntries -= ColOffset;
00490 
00491   if (ierr != 0) std::cout << "WARNING BlockCrsMatrix::BlockExtractGlobalRowView err = "
00492      << ierr << "\n\t  Row " << BaseRow + RowOffset << "Col " << Col+ColOffset << std::endl;
00493 }
00494 
00495 #ifndef EPETRA_NO_32BIT_GLOBAL_INDICES
00496 void BlockCrsMatrix::BlockExtractGlobalRowView(const int BaseRow, int& NumEntries,
00497      double*& Values, const int Row, const int Col)
00498 {
00499   if(Epetra_CrsMatrix::RowMatrixRowMap().GlobalIndicesInt())
00500   return TBlockExtractGlobalRowView<int>(BaseRow, NumEntries, Values, Row, Col);
00501   else
00502     throw "EpetraExt::BlockCrsMatrix::BlockExtractGlobalRowView: Global indices not int";
00503 }
00504 #endif
00505 
00506 #ifndef EPETRA_NO_64BIT_GLOBAL_INDICES
00507 void BlockCrsMatrix::BlockExtractGlobalRowView(const long long BaseRow, int& NumEntries,
00508      double*& Values, const long long Row, const long long Col)
00509 {
00510   if(Epetra_CrsMatrix::RowMatrixRowMap().GlobalIndicesLongLong())
00511   return TBlockExtractGlobalRowView<long long>(BaseRow, NumEntries, Values, Row, Col);
00512   else
00513     throw "EpetraExt::BlockCrsMatrix::BlockExtractGlobalRowView: Global indices not long long";
00514 }
00515 #endif
00516 
00517 //==============================================================================
00518 
00519 template<typename int_type>
00520 void BlockCrsMatrix::TExtractBlock(Epetra_CrsMatrix & BaseMatrix, const int_type Row, const int_type Col)
00521 {
00522   std::vector<int_type>& RowIndices_ = TRowIndices<int_type>();
00523   std::vector< std::vector<int_type> >& RowStencil_ = TRowStencil<int_type>();
00524   int_type RowOffset = RowIndices_[(std::size_t)Row] * ROffset_;
00525   int_type ColOffset = (RowIndices_[(std::size_t)Row] + RowStencil_[(std::size_t)Row][(std::size_t)Col]) * COffset_;
00526 
00527 //  const Epetra_CrsGraph & BaseGraph = BaseMatrix.Graph();
00528   const Epetra_BlockMap & BaseMap = BaseMatrix.RowMatrixRowMap();
00529   //const Epetra_BlockMap & BaseColMap = BaseMatrix.RowMatrixColMap();
00530 
00531   // This routine extracts entries of a BaseMatrix from a big  BlockCrsMatrix
00532   // It performs the following operation on the global IDs row-by-row
00533   // BaseMatrix.val[i][j] = this->val[i+rowOffset][j+ColOffset] 
00534 
00535   int MaxIndices = BaseMatrix.MaxNumEntries();
00536   vector<int_type> Indices(MaxIndices);
00537   vector<double> Values(MaxIndices);
00538   int NumIndices;
00539   int_type indx,icol;
00540   double* BlkValues;
00541   int *BlkIndices;
00542   int BlkNumIndices;
00543   int ierr=0;
00544   (void) ierr; // Forestall compiler warning for unused variable.
00545 
00546   for (int i=0; i<BaseMap.NumMyElements(); i++) {
00547 
00548     // Get pointers to values and indices of whole block matrix row
00549     int_type BaseRow = (int_type) BaseMap.GID64(i);
00550     int myBlkBaseRow = this->RowMatrixRowMap().LID(BaseRow + RowOffset);
00551     ierr = this->ExtractMyRowView(myBlkBaseRow, BlkNumIndices, BlkValues, BlkIndices); 
00552 
00553     NumIndices = 0;
00554     // Grab columns with global indices in correct range for this block
00555     for( int l = 0; l < BlkNumIndices; ++l ) {
00556        icol = (int_type) this->RowMatrixColMap().GID64(BlkIndices[l]);
00557        indx = icol - ColOffset;
00558        if (indx >= 0 && indx < COffset_) {
00559          Indices[NumIndices] = indx;
00560          Values[NumIndices] = BlkValues[l];
00561    NumIndices++;
00562        }
00563     }
00564 
00565     //Load this row into base matrix
00566     BaseMatrix.ReplaceGlobalValues(BaseRow, NumIndices, &Values[0], &Indices[0] );
00567 
00568   }
00569 }
00570 
00571 #ifndef EPETRA_NO_32BIT_GLOBAL_INDICES
00572 void BlockCrsMatrix::ExtractBlock(Epetra_CrsMatrix & BaseMatrix, const int Row, const int Col)
00573 {
00574   if(Epetra_CrsMatrix::RowMatrixRowMap().GlobalIndicesInt() && BaseMatrix.RowMatrixRowMap().GlobalIndicesInt())
00575   return TExtractBlock<int>(BaseMatrix, Row, Col);
00576   else
00577     throw "EpetraExt::BlockCrsMatrix::ExtractBlock: Global indices not int";
00578 }
00579 #endif
00580 
00581 #ifndef EPETRA_NO_64BIT_GLOBAL_INDICES
00582 void BlockCrsMatrix::ExtractBlock(Epetra_CrsMatrix & BaseMatrix, const long long Row, const long long Col)
00583 {
00584   if(Epetra_CrsMatrix::RowMatrixRowMap().GlobalIndicesLongLong() && BaseMatrix.RowMatrixRowMap().GlobalIndicesLongLong())
00585   return TExtractBlock<long long>(BaseMatrix, Row, Col);
00586   else
00587     throw "EpetraExt::BlockCrsMatrix::ExtractBlock: Global indices not long long";
00588 }
00589 #endif
00590 
00591 } //namespace EpetraExt
00592 
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines