EpetraExt_CrsSingletonFilter_LinearProblem.h

Go to the documentation of this file.

//@HEADER
// ***********************************************************************
//
//     EpetraExt: Epetra Extended - 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
                                                                                                    
#ifndef _EpetraExt_LINEARPROBLEM_CRSSINGLETONFILTER_H_
#define _EpetraExt_LINEARPROBLEM_CRSSINGLETONFILTER_H_

#include "Epetra_ConfigDefs.h"
#include "Epetra_Object.h"
#include "Epetra_CrsMatrix.h"
#include "Epetra_MapColoring.h"
#include "Epetra_SerialDenseVector.h"

#include "EpetraExt_Transform.h"
#include "Teuchos_RCP.hpp"

class Epetra_LinearProblem;
class Epetra_Map;
class Epetra_MultiVector;
class Epetra_Import;
class Epetra_Export;
class Epetra_IntVector;

namespace EpetraExt {


class LinearProblem_CrsSingletonFilter : public SameTypeTransform<Epetra_LinearProblem> {
      
 public:


  LinearProblem_CrsSingletonFilter( bool verbose = false );

  virtual ~LinearProblem_CrsSingletonFilter();

  NewTypeRef operator()( OriginalTypeRef orig );

  bool analyze( OriginalTypeRef orig );

  NewTypeRef construct();

  bool fwd();

  bool rvs();



  int Analyze(Epetra_RowMatrix * FullMatrix);

  bool SingletonsDetected() const {if (!AnalysisDone_) return(false); else return(NumSingletons()>0);};



  int ConstructReducedProblem(Epetra_LinearProblem * Problem);


  int UpdateReducedProblem(Epetra_LinearProblem * Problem);




  int ComputeFullSolution();


  int NumRowSingletons() const {return(NumGlobalRowSingletons_);};

  int NumColSingletons() const {return(NumGlobalColSingletons_);};


  int NumSingletons() const {return(NumColSingletons()+NumRowSingletons());};

  double RatioOfDimensions() const {return(RatioOfDimensions_);};

  double RatioOfNonzeros() const {return(RatioOfNonzeros_);};



  Epetra_LinearProblem * FullProblem() const {return(FullProblem_);};

  Epetra_LinearProblem * ReducedProblem() const {return(ReducedProblem_.get());};

  Epetra_RowMatrix * FullMatrix() const {return(FullMatrix_);};

  Epetra_CrsMatrix * ReducedMatrix() const {return(ReducedMatrix_.get());};

  Epetra_MapColoring * RowMapColors() const {return(RowMapColors_);};

  Epetra_MapColoring * ColMapColors() const {return(ColMapColors_);};

  Epetra_Map * ReducedMatrixRowMap() const {return(ReducedMatrixRowMap_);};

  Epetra_Map * ReducedMatrixColMap() const {return(ReducedMatrixColMap_);};

  Epetra_Map * ReducedMatrixDomainMap() const {return(ReducedMatrixDomainMap_);};

  Epetra_Map * ReducedMatrixRangeMap() const {return(ReducedMatrixRangeMap_);};

 protected:

 

  //  This pointer will be zero if full matrix is not a CrsMatrix.
  Epetra_CrsMatrix * FullCrsMatrix() const {return(FullCrsMatrix_);};

  const Epetra_Map & FullMatrixRowMap() const {return(FullMatrix()->RowMatrixRowMap());};
  const Epetra_Map & FullMatrixColMap() const {return(FullMatrix()->RowMatrixColMap());};
  const Epetra_Map & FullMatrixDomainMap() const {return((FullMatrix()->OperatorDomainMap()));};
  const Epetra_Map & FullMatrixRangeMap() const {return((FullMatrix()->OperatorRangeMap()));};
  void InitializeDefaults();
  int ComputeEliminateMaps();
  int Setup(Epetra_LinearProblem * Problem);
  int InitFullMatrixAccess();
  int GetRow(int Row, int & NumIndices, int * & Indices);
  int GetRow(int Row, int & NumIndices, double * & Values, int * & Indices);
#ifndef EPETRA_NO_32BIT_GLOBAL_INDICES
  int GetRowGCIDs(int Row, int & NumIndices, double * & Values, int * & GlobalIndices);
#endif
#ifndef EPETRA_NO_64BIT_GLOBAL_INDICES
  int GetRowGCIDs(int Row, int & NumIndices, double * & Values, long long * & GlobalIndices);
#endif
  int CreatePostSolveArrays(const Epetra_IntVector & RowIDs,
          const Epetra_MapColoring & RowMapColors,
          const Epetra_IntVector & ColProfiles,
          const Epetra_IntVector & NewColProfiles,
          const Epetra_IntVector & ColHasRowWithSingleton);
  
  int ConstructRedistributeExporter(Epetra_Map * SourceMap, Epetra_Map * TargetMap,
            Epetra_Export * & RedistributeExporter,
            Epetra_Map * & RedistributeMap);
  
  Epetra_LinearProblem * FullProblem_;
  Teuchos::RCP<Epetra_LinearProblem> ReducedProblem_;
  Epetra_RowMatrix * FullMatrix_;
  Epetra_CrsMatrix * FullCrsMatrix_;
  Teuchos::RCP<Epetra_CrsMatrix> ReducedMatrix_;
  Epetra_MultiVector * ReducedRHS_;
  Epetra_MultiVector * ReducedLHS_;
  
  Epetra_Map * ReducedMatrixRowMap_;
  Epetra_Map * ReducedMatrixColMap_;
  Epetra_Map * ReducedMatrixDomainMap_;
  Epetra_Map * ReducedMatrixRangeMap_;
  Epetra_Map * OrigReducedMatrixDomainMap_;
  Epetra_Import * Full2ReducedRHSImporter_;
  Epetra_Import * Full2ReducedLHSImporter_;
  Epetra_Export * RedistributeDomainExporter_;
  
  int * ColSingletonRowLIDs_;
  int * ColSingletonColLIDs_;
  int * ColSingletonPivotLIDs_;
  double * ColSingletonPivots_;
  
  
  int AbsoluteThreshold_;
  double RelativeThreshold_;

  int NumMyRowSingletons_;
  int NumMyColSingletons_;
  int NumGlobalRowSingletons_;
  int NumGlobalColSingletons_;
  double RatioOfDimensions_;
  double RatioOfNonzeros_;
  
  bool HaveReducedProblem_;
  bool UserDefinedEliminateMaps_;
  bool AnalysisDone_;
  bool SymmetricElimination_;
  
  Epetra_MultiVector * tempExportX_;
  Epetra_MultiVector * tempX_;
  Epetra_MultiVector * tempB_;
  Epetra_MultiVector * RedistributeReducedLHS_;
  int * Indices_int_;
#ifndef EPETRA_NO_64BIT_GLOBAL_INDICES
  long long * Indices_LL_;
#endif
  Epetra_SerialDenseVector Values_;
  
  Epetra_MapColoring * RowMapColors_;
  Epetra_MapColoring * ColMapColors_;
  bool FullMatrixIsCrsMatrix_;
  int MaxNumMyEntries_;

  bool verbose_;
  
  
 private:
  LinearProblem_CrsSingletonFilter(const LinearProblem_CrsSingletonFilter & Problem){};

  template<typename int_type>
  int TConstructReducedProblem(Epetra_LinearProblem * Problem);

  template<typename int_type>
  int TUpdateReducedProblem(Epetra_LinearProblem * Problem);

  template<typename int_type>
  int TConstructRedistributeExporter(Epetra_Map * SourceMap, Epetra_Map * TargetMap,
            Epetra_Export * & RedistributeExporter,
            Epetra_Map * & RedistributeMap);
};

} //namespace EpetraExt

#endif /* _EpetraExt_LINEARPROBLEM_CRSSINGLETONFILTER_H_ */