EpetraExt_BlockCrsMatrix.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_BlockCrsMatrix.h"
00030 #include "EpetraExt_BlockUtility.h"
00031 #include "Epetra_Map.h"
00032 
00033 namespace EpetraExt {
00034 
00035 using std::vector;
00036 
00037 //==============================================================================
00038 BlockCrsMatrix::BlockCrsMatrix(
00039         const Epetra_CrsGraph & BaseGraph,
00040         const vector<int> & RowStencil,
00041         int RowIndex,
00042         const Epetra_Comm & GlobalComm  ) 
00043   : Epetra_CrsMatrix( Copy, *(BlockUtility::GenerateBlockGraph( BaseGraph, vector< vector<int> >(1,RowStencil), vector<int>(1,RowIndex), GlobalComm )) ),
00044     BaseGraph_( BaseGraph ),
00045     RowStencil_( vector< vector<int> >(1,RowStencil) ),
00046     RowIndices_( vector<int>(1,RowIndex) ),
00047     Offset_(BlockUtility::CalculateOffset(BaseGraph.RowMap()))
00048 {
00049 }
00050 
00051 //==============================================================================
00052 BlockCrsMatrix::BlockCrsMatrix(
00053         const Epetra_CrsGraph & BaseGraph,
00054         const vector< vector<int> > & RowStencil,
00055         const vector<int> & RowIndices,
00056         const Epetra_Comm & GlobalComm  ) 
00057   : Epetra_CrsMatrix( Copy, *(BlockUtility::GenerateBlockGraph( BaseGraph, RowStencil, RowIndices, GlobalComm )) ),
00058     BaseGraph_( BaseGraph ),
00059     RowStencil_( RowStencil ),
00060     RowIndices_( RowIndices ),
00061     Offset_(BlockUtility::CalculateOffset(BaseGraph.RowMap()))
00062 {
00063 }
00064 
00065 //==============================================================================
00066 BlockCrsMatrix::BlockCrsMatrix(
00067         const Epetra_RowMatrix & BaseMatrix,
00068         const vector< vector<int> > & RowStencil,
00069         const vector<int> & RowIndices,
00070         const Epetra_Comm & GlobalComm  ) 
00071   : Epetra_CrsMatrix( Copy, *(BlockUtility::GenerateBlockGraph( BaseMatrix, RowStencil, RowIndices, GlobalComm )) ),
00072     BaseGraph_( Copy, BaseMatrix.RowMatrixRowMap(), 1 ), //Junk to satisfy constructor
00073     RowStencil_( RowStencil ),
00074     RowIndices_( RowIndices ),
00075     Offset_(BlockUtility::CalculateOffset(BaseMatrix.RowMatrixRowMap()))
00076 {
00077 }
00078 
00079 //==============================================================================
00080 BlockCrsMatrix::BlockCrsMatrix( const BlockCrsMatrix & Matrix ) 
00081   : Epetra_CrsMatrix( dynamic_cast<const Epetra_CrsMatrix &>( Matrix ) ),
00082     BaseGraph_( Matrix.BaseGraph_ ),
00083     RowStencil_( Matrix.RowStencil_ ),
00084     RowIndices_( Matrix.RowIndices_ ),
00085     Offset_( Matrix.Offset_ )
00086 {
00087 }
00088 
00089 //==============================================================================
00090 BlockCrsMatrix::~BlockCrsMatrix()
00091 {
00092 }
00093 
00094 //==============================================================================
00095 void BlockCrsMatrix::LoadBlock(const Epetra_RowMatrix & BaseMatrix, const int Row, const int Col)
00096 {
00097   int RowOffset = RowIndices_[Row] * Offset_;
00098   int ColOffset = (RowIndices_[Row] + RowStencil_[Row][Col]) * Offset_;
00099 
00100 //  const Epetra_CrsGraph & BaseGraph = BaseMatrix.Graph();
00101   const Epetra_BlockMap & BaseMap = BaseMatrix.RowMatrixRowMap();
00102   const Epetra_BlockMap & BaseColMap = BaseMatrix.RowMatrixColMap();
00103 
00104   // This routine copies entries of a BaseMatrix into big  BlockCrsMatrix
00105   // It performs the following operation on the global IDs row-by-row
00106   // this->val[i+rowOffset][j+ColOffset] = BaseMatrix.val[i][j]
00107 
00108   int MaxIndices = BaseMatrix.MaxNumEntries();
00109   vector<int> Indices(MaxIndices);
00110   vector<double> Values(MaxIndices);
00111   int NumIndices;
00112   int ierr=0;
00113 
00114   for (int i=0; i<BaseMap.NumMyElements(); i++) {
00115     BaseMatrix.ExtractMyRowCopy( i, MaxIndices, NumIndices, &Values[0], &Indices[0] );
00116 
00117     // Convert to BlockMatrix Global numbering scheme
00118     for( int l = 0; l < NumIndices; ++l )
00119        Indices[l] = ColOffset +  BaseColMap.GID(Indices[l]);
00120 
00121     int BaseRow = BaseMap.GID(i);
00122     ierr = this->ReplaceGlobalValues(BaseRow + RowOffset, NumIndices, &Values[0], &Indices[0]); 
00123     if (ierr != 0) cout << "WARNING BlockCrsMatrix::LoadBlock ReplaceGlobalValues err = " << ierr <<
00124       "\n\t  Row " << BaseRow + RowOffset << "Col start" << Indices[0] << endl;
00125 
00126   }
00127 }
00128 
00129 //==============================================================================
00130   void BlockCrsMatrix::SumIntoBlock(double alpha, const Epetra_RowMatrix & BaseMatrix, const int Row, const int Col)
00131 {
00132   int RowOffset = RowIndices_[Row] * Offset_;
00133   int ColOffset = (RowIndices_[Row] + RowStencil_[Row][Col]) * Offset_;
00134 
00135 //  const Epetra_CrsGraph & BaseGraph = BaseMatrix.Graph();
00136   const Epetra_BlockMap & BaseMap = BaseMatrix.RowMatrixRowMap();
00137   const Epetra_BlockMap & BaseColMap = BaseMatrix.RowMatrixColMap();
00138 
00139   // This routine copies entries of a BaseMatrix into big  BlockCrsMatrix
00140   // It performs the following operation on the global IDs row-by-row
00141   // this->val[i+rowOffset][j+ColOffset] = BaseMatrix.val[i][j]
00142 
00143   int MaxIndices = BaseMatrix.MaxNumEntries();
00144   vector<int> Indices(MaxIndices);
00145   vector<double> Values(MaxIndices);
00146   int NumIndices;
00147   int ierr=0;
00148 
00149   for (int i=0; i<BaseMap.NumMyElements(); i++) {
00150     BaseMatrix.ExtractMyRowCopy( i, MaxIndices, NumIndices, &Values[0], &Indices[0] );
00151 
00152     // Convert to BlockMatrix Global numbering scheme
00153     for( int l = 0; l < NumIndices; ++l ) {
00154        Indices[l] = ColOffset +  BaseColMap.GID(Indices[l]);
00155        Values[l] *= alpha;
00156     }
00157 
00158     int BaseRow = BaseMap.GID(i);
00159     ierr = this->SumIntoGlobalValues(BaseRow + RowOffset, NumIndices, &Values[0], &Indices[0]); 
00160     if (ierr != 0) cout << "WARNING BlockCrsMatrix::SumIntoBlock SumIntoGlobalValues err = " << ierr <<
00161       "\n\t  Row " << BaseRow + RowOffset << "Col start" << Indices[0] << endl;
00162 
00163   }
00164 }
00165 
00166 //==============================================================================
00167 void BlockCrsMatrix::BlockSumIntoGlobalValues(const int BaseRow, int NumIndices,
00168      double* Values, const int* Indices, const int Row, const int Col)
00169 //All arguments could be const, except some were not set as const in CrsMatrix
00170 {
00171   int RowOffset = RowIndices_[Row] * Offset_;
00172   int ColOffset = (RowIndices_[Row] + RowStencil_[Row][Col]) * Offset_;
00173 
00174   // Convert to BlockMatrix Global numbering scheme
00175   vector<int> OffsetIndices(NumIndices);
00176   for( int l = 0; l < NumIndices; ++l ) OffsetIndices[l] = Indices[l] + ColOffset;
00177 
00178   int ierr = this->SumIntoGlobalValues(BaseRow + RowOffset, NumIndices,
00179                                    Values, &OffsetIndices[0]); 
00180 
00181   if (ierr != 0) cout << "WARNING BlockCrsMatrix::BlockSumIntoGlobalValues err = "
00182      << ierr << "\n\t  Row " << BaseRow + RowOffset << "Col start" << Indices[0] << endl;
00183 }
00184 
00185 //==============================================================================
00186 void BlockCrsMatrix::BlockReplaceGlobalValues(const int BaseRow, int NumIndices,
00187      double* Values, const int* Indices, const int Row, const int Col)
00188 //All arguments could be const, except some were not set as const in CrsMatrix
00189 {
00190   int RowOffset = RowIndices_[Row] * Offset_;
00191   int ColOffset = (RowIndices_[Row] + RowStencil_[Row][Col]) * Offset_;
00192 
00193   // Convert to BlockMatrix Global numbering scheme
00194   vector<int> OffsetIndices(NumIndices);
00195   for( int l = 0; l < NumIndices; ++l ) OffsetIndices[l] = Indices[l] + ColOffset;
00196 
00197   int ierr = this->ReplaceGlobalValues(BaseRow + RowOffset, NumIndices,
00198                                    Values, &OffsetIndices[0]); 
00199 
00200   if (ierr != 0) cout << "WARNING BlockCrsMatrix::BlockReplaceGlobalValues err = "
00201      << ierr << "\n\t  Row " << BaseRow + RowOffset << "Col start" << Indices[0] << endl;
00202 }
00203 
00204 //==============================================================================
00205 void BlockCrsMatrix::BlockExtractGlobalRowView(const int BaseRow, 
00206                  int& NumEntries, 
00207                  double*& Values, 
00208                  const int Row, 
00209                  const int Col)
00210 //All arguments could be const, except some were not set as const in CrsMatrix
00211 {
00212   int RowOffset = RowIndices_[Row] * Offset_;
00213   int ColOffset = (RowIndices_[Row] + RowStencil_[Row][Col]) * Offset_;
00214 
00215   // Get the whole row
00216   int ierr = this->ExtractGlobalRowView(BaseRow + RowOffset, NumEntries,
00217           Values); 
00218 
00219   // Adjust for just this block column
00220   Values += ColOffset;
00221   NumEntries -= ColOffset;
00222 
00223   if (ierr != 0) cout << "WARNING BlockCrsMatrix::BlockExtractGlobalRowView err = "
00224      << ierr << "\n\t  Row " << BaseRow + RowOffset << "Col " << Col+ColOffset << endl;
00225 }
00226 
00227 //==============================================================================
00228 void BlockCrsMatrix::ExtractBlock(Epetra_CrsMatrix & BaseMatrix, const int Row, const int Col)
00229 {
00230   int RowOffset = RowIndices_[Row] * Offset_;
00231   int ColOffset = (RowIndices_[Row] + RowStencil_[Row][Col]) * Offset_;
00232 
00233 //  const Epetra_CrsGraph & BaseGraph = BaseMatrix.Graph();
00234   const Epetra_BlockMap & BaseMap = BaseMatrix.RowMatrixRowMap();
00235   //const Epetra_BlockMap & BaseColMap = BaseMatrix.RowMatrixColMap();
00236 
00237   // This routine extracts entries of a BaseMatrix from a big  BlockCrsMatrix
00238   // It performs the following operation on the global IDs row-by-row
00239   // BaseMatrix.val[i][j] = this->val[i+rowOffset][j+ColOffset] 
00240 
00241   int MaxIndices = BaseMatrix.MaxNumEntries();
00242   vector<int> Indices(MaxIndices);
00243   vector<double> Values(MaxIndices);
00244   int NumIndices;
00245   int indx,icol;
00246   double* BlkValues;
00247   int *BlkIndices;
00248   int BlkNumIndices;
00249   int ierr=0;
00250 
00251   for (int i=0; i<BaseMap.NumMyElements(); i++) {
00252 
00253     // Get pointers to values and indices of whole block matrix row
00254     int BaseRow = BaseMap.GID(i);
00255     int myBlkBaseRow = this->RowMatrixRowMap().LID(BaseRow + RowOffset);
00256     ierr = this->ExtractMyRowView(myBlkBaseRow, BlkNumIndices, BlkValues, BlkIndices); 
00257 
00258     NumIndices = 0;
00259     // Grab columns with global indices in correct range for this block
00260     for( int l = 0; l < BlkNumIndices; ++l ) {
00261        icol = this->RowMatrixColMap().GID(BlkIndices[l]);
00262        indx = icol - ColOffset;
00263        if (indx >= 0 && indx < Offset_) {
00264          Indices[NumIndices] = indx;
00265          Values[NumIndices] = BlkValues[l];
00266    NumIndices++;
00267        }
00268     }
00269 
00270     //Load this row into base matrix
00271     BaseMatrix.ReplaceGlobalValues(BaseRow, NumIndices, &Values[0], &Indices[0] );
00272 
00273   }
00274 }
00275 
00276 } //namespace EpetraExt

Generated on Wed May 12 21:24:45 2010 for EpetraExt by  doxygen 1.4.7