Epetra_RowMatrixTransposer.cpp

Go to the documentation of this file.

 
//@HEADER
// ************************************************************************
//
//               Epetra: Linear Algebra Services Package
//                 Copyright 2011 Sandia Corporation
//
// Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
// the U.S. Government retains certain rights in this software.
//
// Redistribution and use in source and binary forms, with or without
// modification, are permitted provided that the following conditions are
// met:
//
// 1. Redistributions of source code must retain the above copyright
// notice, this list of conditions and the following disclaimer.
//
// 2. Redistributions in binary form must reproduce the above copyright
// notice, this list of conditions and the following disclaimer in the
// documentation and/or other materials provided with the distribution.
//
// 3. Neither the name of the Corporation nor the names of the
// contributors may be used to endorse or promote products derived from
// this software without specific prior written permission.
//
// THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
// EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
// PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
// CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
// EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
// PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
// PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
// LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
// NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
// SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
//
// Questions? Contact Michael A. Heroux (maherou@sandia.gov)
//
// ************************************************************************
//@HEADER
 
#include "Epetra_ConfigDefs.h"
#include "Epetra_RowMatrixTransposer.h"
#include "Epetra_RowMatrix.h"
#include "Epetra_CrsMatrix.h"
#include "Epetra_CrsGraph.h"
#include "Epetra_Map.h"
#include "Epetra_Export.h"
 
#ifndef EPETRA_NO_32BIT_GLOBAL_INDICES // FIXME
// FIXME long long : whole file
 
//=============================================================================
Epetra_RowMatrixTransposer::Epetra_RowMatrixTransposer(Epetra_RowMatrix * OrigMatrix)
  : OrigMatrix_(OrigMatrix),
    TransposeMatrix_(0),
    TransposeExporter_(0),
    TransposeRowMap_(0),
    TransposeCreated_(false),
    MakeDataContiguous_(false),
    NumMyRows_(0),
    NumMyCols_(0),
    MaxNumEntries_(0),
    Indices_(NULL),
    Values_(NULL),
    TransNumNz_(NULL),
    TransIndices_(NULL),
    TransValues_(NULL),
    TransMyGlobalEquations_(NULL),
    OrigMatrixIsCrsMatrix_(false)
{
}
//=============================================================================
Epetra_RowMatrixTransposer::Epetra_RowMatrixTransposer(const Epetra_RowMatrixTransposer& Source)
  :OrigMatrix_(Source.OrigMatrix_),
   TransposeMatrix_(0),
   TransposeExporter_(0),
   TransposeRowMap_(0),
   TransposeCreated_(Source.TransposeCreated_),
   MakeDataContiguous_(Source.MakeDataContiguous_),
   NumMyRows_(0),
   NumMyCols_(0),
   MaxNumEntries_(0),
   Indices_(NULL),
   Values_(NULL),
   TransNumNz_(NULL),
   TransIndices_(NULL),
   TransValues_(NULL),
   TransMyGlobalEquations_(NULL),
   OrigMatrixIsCrsMatrix_(false)
{
  TransposeMatrix_ = new Epetra_CrsMatrix(*Source.TransposeMatrix_);
  if (MakeDataContiguous_) TransposeMatrix_->MakeDataContiguous();
  TransposeExporter_ = new Epetra_Export(*Source.TransposeExporter_);
}
//=========================================================================
Epetra_RowMatrixTransposer::~Epetra_RowMatrixTransposer(){
 
  DeleteData();
 
}
 
//=========================================================================
void Epetra_RowMatrixTransposer::DeleteData (){
 
  int i;
 
  if (TransposeExporter_!=0) {delete TransposeExporter_; TransposeExporter_=0;}
 
  // Delete any intermediate storage
 
  if (!OrigMatrixIsCrsMatrix_) {
    delete [] Indices_;
    delete [] Values_;
  }
 
 
  for(i=0; i<NumMyCols_; i++) {
    int NumIndices = TransNumNz_[i];
    if (NumIndices>0) {
      delete [] TransIndices_[i];
      delete [] TransValues_[i];
    }
  }
  delete [] TransNumNz_;
  delete [] TransIndices_;
  delete [] TransValues_;
  delete [] TransMyGlobalEquations_;
}
 
//=========================================================================
int Epetra_RowMatrixTransposer::CreateTranspose (const bool MakeDataContiguous,
             Epetra_CrsMatrix *& TransposeMatrix,
             Epetra_Map * TransposeRowMap_in) {
 
// FIXME long long
 
  int i, j;
 
  if (TransposeCreated_) DeleteData(); // Get rid of existing data first
 
  if (TransposeRowMap_in==0)
    TransposeRowMap_ = (Epetra_Map *) &(OrigMatrix_->OperatorDomainMap()); // Should be replaced with refcount =
  else
    TransposeRowMap_ = TransposeRowMap_in;
 
  // This routine will work for any RowMatrix object, but will attempt cast the matrix to a CrsMatrix if
  // possible (because we can then use a View of the matrix and graph, which is much cheaper).
 
  // First get the local indices to count how many nonzeros will be in the
  // transpose graph on each processor
 
 
  Epetra_CrsMatrix * OrigCrsMatrix = dynamic_cast<Epetra_CrsMatrix *>(OrigMatrix_);
 
  OrigMatrixIsCrsMatrix_ = (OrigCrsMatrix!=0); // If this pointer is non-zero, the cast to CrsMatrix worked
 
  NumMyRows_ = OrigMatrix_->NumMyRows();
  NumMyCols_ = OrigMatrix_->NumMyCols();
  NumMyRows_ = OrigMatrix_->NumMyRows();
  TransNumNz_ = new int[NumMyCols_];
  TransIndices_ = new int*[NumMyCols_];
  TransValues_ = new double*[NumMyCols_];
 
 
  int NumIndices;
 
  if (OrigMatrixIsCrsMatrix_) {
 
 
    const Epetra_CrsGraph & OrigGraph = OrigCrsMatrix->Graph(); // Get matrix graph
 
    for (i=0;i<NumMyCols_; i++) TransNumNz_[i] = 0;
    for (i=0; i<NumMyRows_; i++) {
      EPETRA_CHK_ERR(OrigGraph.ExtractMyRowView(i, NumIndices, Indices_)); // Get view of ith row
      for (j=0; j<NumIndices; j++) ++TransNumNz_[Indices_[j]];
    }
  }
  else { // OrigMatrix is not a CrsMatrix
 
    MaxNumEntries_ = 0;
    int NumEntries;
    for (i=0; i<NumMyRows_; i++) {
      OrigMatrix_->NumMyRowEntries(i, NumEntries);
      MaxNumEntries_ = EPETRA_MAX(MaxNumEntries_, NumEntries);
    }
    Indices_ = new int[MaxNumEntries_];
    Values_ = new double[MaxNumEntries_];
 
    for (i=0;i<NumMyCols_; i++) TransNumNz_[i] = 0;
    for (i=0; i<NumMyRows_; i++) {
      // Get ith row
      EPETRA_CHK_ERR(OrigMatrix_->ExtractMyRowCopy(i, MaxNumEntries_, NumIndices, Values_, Indices_));
      for (j=0; j<NumIndices; j++) ++TransNumNz_[Indices_[j]];
    }
  }
 
 
  // Most of remaining code is common to both cases
 
  for(i=0; i<NumMyCols_; i++) {
    NumIndices = TransNumNz_[i];
    if (NumIndices>0) {
      TransIndices_[i] = new int[NumIndices];
      TransValues_[i] = new double[NumIndices];
    }
  }
 
  // Now copy values and global indices into newly created transpose storage
 
  for (i=0;i<NumMyCols_; i++) TransNumNz_[i] = 0; // Reset transpose NumNz counter
  for (i=0; i<NumMyRows_; i++) {
    if (OrigMatrixIsCrsMatrix_) {
      EPETRA_CHK_ERR(OrigCrsMatrix->ExtractMyRowView(i, NumIndices, Values_, Indices_));
    }
    else {
      EPETRA_CHK_ERR(OrigMatrix_->ExtractMyRowCopy(i, MaxNumEntries_, NumIndices, Values_, Indices_));
    }
 
    int ii = OrigMatrix_->RowMatrixRowMap().GID64(i); // FIXME long long
    for (j=0; j<NumIndices; j++) {
      int TransRow = Indices_[j];
      int loc = TransNumNz_[TransRow];
      TransIndices_[TransRow][loc] = ii;
      TransValues_[TransRow][loc] = Values_[j];
      ++TransNumNz_[TransRow]; // increment counter into current transpose row
    }
  }
 
  //  Build Transpose matrix with some rows being shared across processors.
  //  We will use a view here since the matrix will not be used for anything else
 
  const Epetra_Map & TransMap = OrigMatrix_->RowMatrixColMap();
 
  Epetra_CrsMatrix TempTransA1(View, TransMap, TransNumNz_);
  TransMyGlobalEquations_ = new int[NumMyCols_];
 
  TransMap.MyGlobalElements(TransMyGlobalEquations_);
 
  /* Add  rows one-at-a-time */
 
  for (i=0; i<NumMyCols_; i++)
    {
      EPETRA_CHK_ERR(TempTransA1.InsertGlobalValues(TransMyGlobalEquations_[i],
                TransNumNz_[i], TransValues_[i], TransIndices_[i]));
    }
  // Note: The following call to FillComplete is currently necessary because
  //       some global constants that are needed by the Export () are computed in this routine
 
  const Epetra_Map& domain_map = OrigMatrix_->OperatorDomainMap();
  const Epetra_Map& range_map = OrigMatrix_->OperatorRangeMap();
 
  EPETRA_CHK_ERR(TempTransA1.FillComplete(range_map, domain_map, false));
 
  // Now that transpose matrix with shared rows is entered, create a new matrix that will
  // get the transpose with uniquely owned rows (using the same row distribution as A).
 
  TransposeMatrix_ = new Epetra_CrsMatrix(Copy, *TransposeRowMap_,0);
 
  // Create an Export object that will move TempTransA around
 
  TransposeExporter_ = new Epetra_Export(TransMap, *TransposeRowMap_);
 
  EPETRA_CHK_ERR(TransposeMatrix_->Export(TempTransA1, *TransposeExporter_, Add));
 
  EPETRA_CHK_ERR(TransposeMatrix_->FillComplete(range_map, domain_map));
 
  if (MakeDataContiguous) {
    EPETRA_CHK_ERR(TransposeMatrix_->MakeDataContiguous());
  }
 
  TransposeMatrix = TransposeMatrix_;
  TransposeCreated_ = true;
 
  return(0);
}
//=========================================================================
int Epetra_RowMatrixTransposer::UpdateTransposeValues(Epetra_RowMatrix * MatrixWithNewValues){
 
  int i, j, NumIndices;
 
  if (!TransposeCreated_) EPETRA_CHK_ERR(-1); // Transpose must be already created
 
  // Sanity check of incoming matrix.  Perform some tests to see if it is compatible with original input matrix
  if (OrigMatrix_!=MatrixWithNewValues) { // Check if pointer of new matrix is same as previous input matrix
    OrigMatrix_ = MatrixWithNewValues; // Reset this pointer if not, then check for other attributes
    if (NumMyRows_ != OrigMatrix_->NumMyRows() ||
  NumMyCols_ != OrigMatrix_->NumMyCols() ||
  NumMyRows_ != OrigMatrix_->NumMyRows()) {
      EPETRA_CHK_ERR(-2); // New matrix not compatible with previous
    }
  }
 
  Epetra_CrsMatrix * OrigCrsMatrix = dynamic_cast<Epetra_CrsMatrix *>(MatrixWithNewValues);
 
 
  OrigMatrixIsCrsMatrix_ = (OrigCrsMatrix!=0); // If this pointer is non-zero, the cast to CrsMatrix worked
 
 
  // Now copy values and global indices into newly create transpose storage
 
  for (i=0;i<NumMyCols_; i++) TransNumNz_[i] = 0; // Reset transpose NumNz counter
  for (i=0; i<NumMyRows_; i++) {
    if (OrigMatrixIsCrsMatrix_) {
      EPETRA_CHK_ERR(OrigCrsMatrix->ExtractMyRowView(i, NumIndices, Values_, Indices_));
    }
    else {
      EPETRA_CHK_ERR(OrigMatrix_->ExtractMyRowCopy(i, MaxNumEntries_, NumIndices, Values_, Indices_));
    }
 
    int ii = OrigMatrix_->RowMatrixRowMap().GID64(i); // FIXME long long
    for (j=0; j<NumIndices; j++) {
      int TransRow = Indices_[j];
      int loc = TransNumNz_[TransRow];
      TransIndices_[TransRow][loc] = ii;
      TransValues_[TransRow][loc] = Values_[j];
      ++TransNumNz_[TransRow]; // increment counter into current transpose row
    }
  }
 
  //  Build Transpose matrix with some rows being shared across processors.
  //  We will use a view here since the matrix will not be used for anything else
 
  const Epetra_Map & TransMap = OrigMatrix_->RowMatrixColMap();
 
  Epetra_CrsMatrix TempTransA1(View, TransMap, TransNumNz_);
  TransMap.MyGlobalElements(TransMyGlobalEquations_);
  /* Add  rows one-at-a-time */
 
  for (i=0; i<NumMyCols_; i++)
    {
      EPETRA_CHK_ERR(TempTransA1.InsertGlobalValues(TransMyGlobalEquations_[i],
                TransNumNz_[i], TransValues_[i], TransIndices_[i]));
    }
  // Note: The following call to FillComplete is currently necessary because
  //       some global constants that are needed by the Export () are computed in this routine
  const Epetra_Map& domain_map = OrigMatrix_->OperatorDomainMap();
  const Epetra_Map& range_map = OrigMatrix_->OperatorRangeMap();
 
  EPETRA_CHK_ERR(TempTransA1.FillComplete(range_map, domain_map, false));
 
  // Now that transpose matrix with shared rows is entered, update values of target transpose matrix
 
  TransposeMatrix_->PutScalar(0.0);  // Zero out all values of the matrix
 
  EPETRA_CHK_ERR(TransposeMatrix_->Export(TempTransA1, *TransposeExporter_, Add));
 
  return(0);
}
 
Epetra_RowMatrixTransposer&
Epetra_RowMatrixTransposer::operator=(const Epetra_RowMatrixTransposer& src)
{
  (void)src;//prevents unused variable compiler warning
 
  //not currently supported
  bool throw_error = true;
  if (throw_error) {
    std::cerr << std::endl
        << "Epetra_RowMatrixTransposer::operator= not supported."
        <<std::endl;
    throw -1;
  }
 
  return(*this);
}
 
#endif // EPETRA_NO_32BIT_GLOBAL_INDICES