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::ExtractBlock(Epetra_CrsMatrix & 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 extracts entries of a BaseMatrix from a big  BlockCrsMatrix
00140   // It performs the following operation on the global IDs row-by-row
00141   // BaseMatrix.val[i][j] = this->val[i+rowOffset][j+ColOffset] 
00142 
00143   int MaxIndices = BaseMatrix.MaxNumEntries();
00144   vector<int> Indices(MaxIndices);
00145   vector<double> Values(MaxIndices);
00146   int NumIndices;
00147   int indx,icol;
00148   double* BlkValues;
00149   int *BlkIndices;
00150   int BlkNumIndices;
00151   int ierr=0;
00152 
00153   for (int i=0; i<BaseMap.NumMyElements(); i++) {
00154 
00155     // Get pointers to values and indices of whole block matrix row
00156     int BaseRow = BaseMap.GID(i);
00157     int myBlkBaseRow = this->RowMatrixRowMap().LID(BaseRow + RowOffset);
00158     ierr = this->ExtractMyRowView(myBlkBaseRow, BlkNumIndices, BlkValues, BlkIndices); 
00159 
00160     NumIndices = 0;
00161     // Grab columns with global indices in correct range for this block
00162     for( int l = 0; l < BlkNumIndices; ++l ) {
00163        icol = this->RowMatrixColMap().GID(BlkIndices[l]);
00164        indx = icol - ColOffset;
00165        if (indx >= 0 && indx < Offset_) {
00166          Indices[NumIndices] = indx;
00167          Values[NumIndices] = BlkValues[l];
00168    NumIndices++;
00169        }
00170     }
00171 
00172     //Load this row into base matrix
00173     BaseMatrix.ReplaceGlobalValues(BaseRow, NumIndices, &Values[0], &Indices[0] );
00174 
00175   }
00176 }
00177 
00178 } //namespace EpetraExt

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