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 BlockCrsMatrix::BlockCrsMatrix(
00052         const Epetra_CrsGraph & BaseGraph,
00053         const vector<int> & RowStencil,
00054         int RowIndex,
00055         const Epetra_Comm & GlobalComm  ) 
00056   : Epetra_CrsMatrix( Copy, *(BlockUtility::GenerateBlockGraph( BaseGraph, vector< vector<int> >(1,RowStencil), vector<int>(1,RowIndex), GlobalComm )) ),
00057     BaseGraph_( BaseGraph ),
00058     RowStencil_( vector< vector<int> >(1,RowStencil) ),
00059     RowIndices_( vector<int>(1,RowIndex) ),
00060     ROffset_(BlockUtility::CalculateOffset(BaseGraph.RowMap())),
00061     COffset_(BlockUtility::CalculateOffset(BaseGraph.ColMap()))
00062 {
00063 }
00064 
00065 //==============================================================================
00066 BlockCrsMatrix::BlockCrsMatrix(
00067         const Epetra_CrsGraph & BaseGraph,
00068         const vector< vector<int> > & RowStencil,
00069         const vector<int> & RowIndices,
00070         const Epetra_Comm & GlobalComm  ) 
00071   : Epetra_CrsMatrix( Copy, *(BlockUtility::GenerateBlockGraph( BaseGraph, RowStencil, RowIndices, GlobalComm )) ),
00072     BaseGraph_( BaseGraph ),
00073     RowStencil_( RowStencil ),
00074     RowIndices_( RowIndices ),
00075     ROffset_(BlockUtility::CalculateOffset(BaseGraph.RowMap())),
00076     COffset_(BlockUtility::CalculateOffset(BaseGraph.ColMap()))
00077 {
00078 }
00079 
00080 //==============================================================================
00081 BlockCrsMatrix::BlockCrsMatrix(
00082         const Epetra_CrsGraph & BaseGraph,
00083         const Epetra_CrsGraph & LocalBlockGraph,
00084         const Epetra_Comm & GlobalComm  ) 
00085   : Epetra_CrsMatrix( Copy, *(BlockUtility::GenerateBlockGraph( BaseGraph, LocalBlockGraph, GlobalComm )) ),
00086     BaseGraph_( BaseGraph ),
00087     RowStencil_( ),
00088     RowIndices_( ),
00089     ROffset_(BlockUtility::CalculateOffset(BaseGraph.RowMap())),
00090     COffset_(BlockUtility::CalculateOffset(BaseGraph.ColMap()))
00091 {
00092   BlockUtility::GenerateRowStencil(LocalBlockGraph, RowIndices_, RowStencil_);
00093 }
00094 
00095 //==============================================================================
00096 BlockCrsMatrix::BlockCrsMatrix(
00097         const Epetra_RowMatrix & BaseMatrix,
00098         const vector< vector<int> > & RowStencil,
00099         const vector<int> & RowIndices,
00100         const Epetra_Comm & GlobalComm  ) 
00101   : Epetra_CrsMatrix( Copy, *(BlockUtility::GenerateBlockGraph( BaseMatrix, RowStencil, RowIndices, GlobalComm )) ),
00102     BaseGraph_( Copy, BaseMatrix.RowMatrixRowMap(), 1 ), //Junk to satisfy constructor
00103     RowStencil_( RowStencil ),
00104     RowIndices_( RowIndices ),
00105     ROffset_(BlockUtility::CalculateOffset(BaseMatrix.RowMatrixRowMap())),
00106     COffset_(BlockUtility::CalculateOffset(BaseMatrix.RowMatrixColMap()))
00107 {
00108 }
00109 
00110 //==============================================================================
00111 BlockCrsMatrix::BlockCrsMatrix( const BlockCrsMatrix & Matrix ) 
00112   : Epetra_CrsMatrix( dynamic_cast<const Epetra_CrsMatrix &>( Matrix ) ),
00113     BaseGraph_( Matrix.BaseGraph_ ),
00114     RowStencil_( Matrix.RowStencil_ ),
00115     RowIndices_( Matrix.RowIndices_ ),
00116     ROffset_( Matrix.ROffset_ ),
00117     COffset_( Matrix.COffset_ )
00118 {
00119 }
00120 
00121 //==============================================================================
00122 BlockCrsMatrix::~BlockCrsMatrix()
00123 {
00124 }
00125 
00126 //==============================================================================
00127 void BlockCrsMatrix::LoadBlock(const Epetra_RowMatrix & BaseMatrix, const int Row, const int Col)
00128 {
00129   int RowOffset = RowIndices_[Row] * ROffset_;
00130   int ColOffset = (RowIndices_[Row] + RowStencil_[Row][Col]) * COffset_;
00131 
00132 //  const Epetra_CrsGraph & BaseGraph = BaseMatrix.Graph();
00133   const Epetra_BlockMap & BaseMap = BaseMatrix.RowMatrixRowMap();
00134   const Epetra_BlockMap & BaseColMap = BaseMatrix.RowMatrixColMap();
00135 
00136   // This routine copies entries of a BaseMatrix into big  BlockCrsMatrix
00137   // It performs the following operation on the global IDs row-by-row
00138   // this->val[i+rowOffset][j+ColOffset] = BaseMatrix.val[i][j]
00139 
00140   int MaxIndices = BaseMatrix.MaxNumEntries();
00141   vector<int> Indices(MaxIndices);
00142   vector<double> Values(MaxIndices);
00143   int NumIndices;
00144   int ierr=0;
00145 
00146   for (int i=0; i<BaseMap.NumMyElements(); i++) {
00147     BaseMatrix.ExtractMyRowCopy( i, MaxIndices, NumIndices, &Values[0], &Indices[0] );
00148 
00149     // Convert to BlockMatrix Global numbering scheme
00150     for( int l = 0; l < NumIndices; ++l )
00151        Indices[l] = ColOffset +  BaseColMap.GID(Indices[l]);
00152 
00153     int BaseRow = BaseMap.GID(i);
00154     ierr = this->ReplaceGlobalValues(BaseRow + RowOffset, NumIndices, &Values[0], &Indices[0]); 
00155     if (ierr != 0) cout << "WARNING BlockCrsMatrix::LoadBlock ReplaceGlobalValues err = " << ierr <<
00156       "\n\t  Row " << BaseRow + RowOffset << "Col start" << Indices[0] << endl;
00157 
00158   }
00159 }
00160 
00161 //==============================================================================
00162   void BlockCrsMatrix::SumIntoBlock(double alpha, const Epetra_RowMatrix & BaseMatrix, const int Row, const int Col)
00163 {
00164   int RowOffset = RowIndices_[Row] * ROffset_;
00165   int ColOffset = (RowIndices_[Row] + RowStencil_[Row][Col]) * COffset_;
00166 
00167 //  const Epetra_CrsGraph & BaseGraph = BaseMatrix.Graph();
00168   const Epetra_BlockMap & BaseMap = BaseMatrix.RowMatrixRowMap();
00169   const Epetra_BlockMap & BaseColMap = BaseMatrix.RowMatrixColMap();
00170 
00171   // This routine copies entries of a BaseMatrix into big  BlockCrsMatrix
00172   // It performs the following operation on the global IDs row-by-row
00173   // this->val[i+rowOffset][j+ColOffset] = BaseMatrix.val[i][j]
00174 
00175   int MaxIndices = BaseMatrix.MaxNumEntries();
00176   vector<int> Indices(MaxIndices);
00177   vector<double> Values(MaxIndices);
00178   int NumIndices;
00179   int ierr=0;
00180 
00181   for (int i=0; i<BaseMap.NumMyElements(); i++) {
00182     BaseMatrix.ExtractMyRowCopy( i, MaxIndices, NumIndices, &Values[0], &Indices[0] );
00183 
00184     // Convert to BlockMatrix Global numbering scheme
00185     for( int l = 0; l < NumIndices; ++l ) {
00186        Indices[l] = ColOffset +  BaseColMap.GID(Indices[l]);
00187        Values[l] *= alpha;
00188     }
00189 
00190     int BaseRow = BaseMap.GID(i);
00191     ierr = this->SumIntoGlobalValues(BaseRow + RowOffset, NumIndices, &Values[0], &Indices[0]); 
00192     if (ierr != 0) cout << "WARNING BlockCrsMatrix::SumIntoBlock SumIntoGlobalValues err = " << ierr <<
00193       "\n\t  Row " << BaseRow + RowOffset << "Col start" << Indices[0] << endl;
00194 
00195   }
00196 }
00197 
00198 //==============================================================================
00199   void BlockCrsMatrix::SumIntoGlobalBlock(double alpha, const Epetra_RowMatrix & BaseMatrix, const int Row, const int Col)
00200 {
00201   int RowOffset = Row * ROffset_;
00202   int ColOffset = Col * COffset_;
00203 
00204 //  const Epetra_CrsGraph & BaseGraph = BaseMatrix.Graph();
00205   const Epetra_BlockMap & BaseMap = BaseMatrix.RowMatrixRowMap();
00206   const Epetra_BlockMap & BaseColMap = BaseMatrix.RowMatrixColMap();
00207 
00208   // This routine copies entries of a BaseMatrix into big  BlockCrsMatrix
00209   // It performs the following operation on the global IDs row-by-row
00210   // this->val[i+rowOffset][j+ColOffset] = BaseMatrix.val[i][j]
00211 
00212   int MaxIndices = BaseMatrix.MaxNumEntries();
00213   vector<int> Indices(MaxIndices);
00214   vector<double> Values(MaxIndices);
00215   int NumIndices;
00216   int ierr=0;
00217 
00218   for (int i=0; i<BaseMap.NumMyElements(); i++) {
00219     BaseMatrix.ExtractMyRowCopy( i, MaxIndices, NumIndices, &Values[0], &Indices[0] );
00220 
00221     // Convert to BlockMatrix Global numbering scheme
00222     for( int l = 0; l < NumIndices; ++l ) {
00223        Indices[l] = ColOffset +  BaseColMap.GID(Indices[l]);
00224        Values[l] *= alpha;
00225     }
00226 
00227     int BaseRow = BaseMap.GID(i);
00228     ierr = this->SumIntoGlobalValues(BaseRow + RowOffset, NumIndices, &Values[0], &Indices[0]); 
00229     if (ierr != 0) cout << "WARNING BlockCrsMatrix::SumIntoBlock SumIntoGlobalValues err = " << ierr <<
00230       "\n\t  Row " << BaseRow + RowOffset << "Col start" << Indices[0] << endl;
00231 
00232   }
00233 }
00234 
00235 //==============================================================================
00236 void BlockCrsMatrix::BlockSumIntoGlobalValues(const int BaseRow, int NumIndices,
00237      double* Values, const int* Indices, const int Row, const int Col)
00238 //All arguments could be const, except some were not set as const in CrsMatrix
00239 {
00240   int RowOffset = RowIndices_[Row] * ROffset_;
00241   int ColOffset = (RowIndices_[Row] + RowStencil_[Row][Col]) * COffset_;
00242 
00243   // Convert to BlockMatrix Global numbering scheme
00244   vector<int> OffsetIndices(NumIndices);
00245   for( int l = 0; l < NumIndices; ++l ) OffsetIndices[l] = Indices[l] + ColOffset;
00246 
00247   int ierr = this->SumIntoGlobalValues(BaseRow + RowOffset, NumIndices,
00248                                    Values, &OffsetIndices[0]); 
00249 
00250   if (ierr != 0) cout << "WARNING BlockCrsMatrix::BlockSumIntoGlobalValues err = "
00251      << ierr << "\n\t  Row " << BaseRow + RowOffset << "Col start" << Indices[0] << endl;
00252 }
00253 
00254 //==============================================================================
00255 void BlockCrsMatrix::BlockReplaceGlobalValues(const int BaseRow, int NumIndices,
00256      double* Values, const int* Indices, const int Row, const int Col)
00257 //All arguments could be const, except some were not set as const in CrsMatrix
00258 {
00259   int RowOffset = RowIndices_[Row] * ROffset_;
00260   int ColOffset = (RowIndices_[Row] + RowStencil_[Row][Col]) * COffset_;
00261 
00262   // Convert to BlockMatrix Global numbering scheme
00263   vector<int> OffsetIndices(NumIndices);
00264   for( int l = 0; l < NumIndices; ++l ) OffsetIndices[l] = Indices[l] + ColOffset;
00265 
00266   int ierr = this->ReplaceGlobalValues(BaseRow + RowOffset, NumIndices,
00267                                    Values, &OffsetIndices[0]); 
00268 
00269   if (ierr != 0) cout << "WARNING BlockCrsMatrix::BlockReplaceGlobalValues err = "
00270      << ierr << "\n\t  Row " << BaseRow + RowOffset << "Col start" << Indices[0] << endl;
00271 }
00272 
00273 //==============================================================================
00274 void BlockCrsMatrix::BlockExtractGlobalRowView(const int BaseRow, 
00275                  int& NumEntries, 
00276                  double*& Values, 
00277                  const int Row, 
00278                  const int Col)
00279 //All arguments could be const, except some were not set as const in CrsMatrix
00280 {
00281   int RowOffset = RowIndices_[Row] * ROffset_;
00282   int ColOffset = (RowIndices_[Row] + RowStencil_[Row][Col]) * COffset_;
00283 
00284   // Get the whole row
00285   int ierr = this->ExtractGlobalRowView(BaseRow + RowOffset, NumEntries,
00286           Values); 
00287 
00288   // Adjust for just this block column
00289   Values += ColOffset;
00290   NumEntries -= ColOffset;
00291 
00292   if (ierr != 0) cout << "WARNING BlockCrsMatrix::BlockExtractGlobalRowView err = "
00293      << ierr << "\n\t  Row " << BaseRow + RowOffset << "Col " << Col+ColOffset << endl;
00294 }
00295 
00296 //==============================================================================
00297 void BlockCrsMatrix::ExtractBlock(Epetra_CrsMatrix & BaseMatrix, const int Row, const int Col)
00298 {
00299   int RowOffset = RowIndices_[Row] * ROffset_;
00300   int ColOffset = (RowIndices_[Row] + RowStencil_[Row][Col]) * COffset_;
00301 
00302 //  const Epetra_CrsGraph & BaseGraph = BaseMatrix.Graph();
00303   const Epetra_BlockMap & BaseMap = BaseMatrix.RowMatrixRowMap();
00304   //const Epetra_BlockMap & BaseColMap = BaseMatrix.RowMatrixColMap();
00305 
00306   // This routine extracts entries of a BaseMatrix from a big  BlockCrsMatrix
00307   // It performs the following operation on the global IDs row-by-row
00308   // BaseMatrix.val[i][j] = this->val[i+rowOffset][j+ColOffset] 
00309 
00310   int MaxIndices = BaseMatrix.MaxNumEntries();
00311   vector<int> Indices(MaxIndices);
00312   vector<double> Values(MaxIndices);
00313   int NumIndices;
00314   int indx,icol;
00315   double* BlkValues;
00316   int *BlkIndices;
00317   int BlkNumIndices;
00318   int ierr=0;
00319   (void) ierr; // Forestall compiler warning for unused variable.
00320 
00321   for (int i=0; i<BaseMap.NumMyElements(); i++) {
00322 
00323     // Get pointers to values and indices of whole block matrix row
00324     int BaseRow = BaseMap.GID(i);
00325     int myBlkBaseRow = this->RowMatrixRowMap().LID(BaseRow + RowOffset);
00326     ierr = this->ExtractMyRowView(myBlkBaseRow, BlkNumIndices, BlkValues, BlkIndices); 
00327 
00328     NumIndices = 0;
00329     // Grab columns with global indices in correct range for this block
00330     for( int l = 0; l < BlkNumIndices; ++l ) {
00331        icol = this->RowMatrixColMap().GID(BlkIndices[l]);
00332        indx = icol - ColOffset;
00333        if (indx >= 0 && indx < COffset_) {
00334          Indices[NumIndices] = indx;
00335          Values[NumIndices] = BlkValues[l];
00336    NumIndices++;
00337        }
00338     }
00339 
00340     //Load this row into base matrix
00341     BaseMatrix.ReplaceGlobalValues(BaseRow, NumIndices, &Values[0], &Indices[0] );
00342 
00343   }
00344 }
00345 
00346 } //namespace EpetraExt
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines