Condest et al.

Makefile

@@ -38,11 +38,14 @@ print_info:
 	@echo "   Trilinos_BUILD_SHARED_LIBS = $(Trilinos_BUILD_SHARED_LIBS)"
 	@echo "End of Trilinos details\n"
 
-libmatrix.mpi: libmatrix.cpp
-	$(CXX) $(CXX_FLAGS) $< -o $@ $(LINK_FLAGS) $(INCLUDE_DIRS) $(LIBRARY_DIRS) $(LIBRARIES)
+%.o: %.cpp
+	$(CXX) -c $(CXX_FLAGS) $< -o $@ $(INCLUDE_DIRS)
+
+libmatrix.mpi: libmatrix.o utils.o
+	$(CXX) $(CXX_FLAGS) $^ -o $@ $(LINK_FLAGS) $(INCLUDE_DIRS) $(LIBRARY_DIRS) $(LIBRARIES)
 
 clean:
-	rm -f libmatrix.mpi *.pyc
+	rm -f libmatrix.mpi *.o *.pyc
 
 test: libmatrix.mpi
 	python test.py

libmatrix.cpp

@@ -1,20 +1,26 @@
-#include <Tpetra_DefaultPlatform.hpp>
-#include <Tpetra_Vector.hpp>
-#include <Tpetra_Version.hpp>
-#include <Tpetra_CrsMatrix.hpp>
+#include "typedefs.h"
+#include "utils.h"
 
-#include <BelosTpetraAdapter.hpp>
-#include <BelosSolverFactory.hpp>
+#include <mpi.h>
 
 #include <Teuchos_XMLParameterListCoreHelpers.hpp>
 #include <Teuchos_GlobalMPISession.hpp>
 #include <Teuchos_oblackholestream.hpp>
 #include <Teuchos_DefaultMpiComm.hpp>
 #include <Teuchos_oblackholestream.hpp>
 
+#include <Tpetra_DefaultPlatform.hpp>
+#include <Tpetra_Vector.hpp>
+#include <Tpetra_Version.hpp>
+#include <Tpetra_CrsMatrix.hpp>
+
+#include <BelosTpetraAdapter.hpp>
+#include <BelosSolverFactory.hpp>
+
 #include <Ifpack2_Factory.hpp>
 
-#include <mpi.h>
+typedef int handle_t;
+typedef uint8_t token_t;
 
 #define ASSERT( cond ) if ( ! (cond) ) throw( "assertion failed: " #cond );
 
@@ -53,26 +59,19 @@ inline bool contains( T iterable, V value ) {
   return false;
 }
 
-typedef double scalar_t;
-typedef int number_t;
-typedef int handle_t;
-typedef int local_t;
-typedef long global_t;
-typedef uint8_t token_t;
-typedef bool bool_t;
-
 typedef Kokkos::DefaultNode::DefaultNodeType node_t;
 typedef Tpetra::Map<local_t,global_t,node_t> map_t;
-typedef Tpetra::Vector<scalar_t,local_t,global_t,node_t> vector_t;
 typedef Tpetra::Operator<scalar_t,local_t,global_t,node_t> operator_t;
 typedef Tpetra::MultiVector<scalar_t,local_t,global_t,node_t> multivector_t;
 typedef Tpetra::RowMatrix<scalar_t,local_t,global_t,node_t> rowmatrix_t;
 typedef Tpetra::CrsMatrix<scalar_t,local_t,global_t,node_t> crsmatrix_t;
+typedef Tpetra::Vector<scalar_t,local_t,global_t,node_t> vector_t;
 typedef Tpetra::CrsGraph<local_t,global_t,node_t> crsgraph_t;
 typedef Tpetra::RowGraph<local_t,global_t,node_t> rowgraph_t;
 typedef Belos::SolverFactory<scalar_t,multivector_t,operator_t> solverfactory_t;
 typedef Teuchos::ScalarTraits<scalar_t>::magnitudeType magnitude_t;
 typedef Tpetra::Export<local_t,global_t,node_t> export_t;
+typedef Ifpack2::Preconditioner<scalar_t,local_t,global_t,node_t> precon_t;
 
 inline Teuchos::Array<local_t> get_none_nan_entries( Teuchos::RCP<const vector_t>& nanvec )
 {
@@ -566,6 +565,23 @@ class LibMatrix : public Intercomm {
 
     gather( &norm );
   }
+
+  void vector_abs() /* compute abs(vec)
+
+       -> broadcast HANDLE handle.{vec,abs}
+  */{
+
+    struct { handle_t vec, abs; } handle;
+    bcast( &handle );
+
+    auto vec = objects.get<vector_t>( handle.vec, out(DEBUG) );
+
+    out(INFO) << "creating vector #" << handle.abs << " from vector #" << handle.vec << std::endl;
+
+    auto absvec = utils::vector::abs( vec );
+
+    objects.set( handle.abs, absvec, out(DEBUG) );
+  }
 
   void vector_sum() /* compute the sum of a vector
 
@@ -585,6 +601,25 @@ class LibMatrix : public Intercomm {
     gather( &sum );
   }
 
+  void vector_max() /* compute the max of a vector
+
+       -> broadcast HANDLE handle.vector
+      <-  gather SCALAR max
+  */{
+
+    struct { handle_t vector; } handle;
+    bcast( &handle );
+    auto vector = objects.get<vector_t>( handle.vector, out(DEBUG) );
+
+    scalar_t max = NAN;
+    for ( auto const &v : vector->getData() )
+    {
+      if( std::isnan(max) || v>max ) max=v;
+    }
+
+    gather( &max );
+  }
+
   void vector_complete() /* export vector
 
        -> broadcast HANDLE handle.{vector,builder,exporter}
@@ -901,6 +936,164 @@ class LibMatrix : public Intercomm {
     matrix->apply( *rhs, *lhs );
   }
 
+  void matrix_condest() /* compute the operator's condition number
+
+      -> broadcast HANDLE handle.operator
+      -> gather SCALAR condest
+  */{
+
+    struct { handle_t oper, params, solvertype, precontype; } handle;
+    bcast( &handle );
+
+    auto B = objects.get<const crsmatrix_t>( handle.oper, out(DEBUG) );
+
+    auto params = objects.get<Teuchos::ParameterList>( handle.params, out(DEBUG) );
+
+    //Create the solver
+    solverfactory_t factory;
+    auto solverparams = Teuchos::sublist(params,"Solver Parameters");
+    auto solvermgr = factory.create( factory.supportedSolverNames()[handle.solvertype], solverparams );
+
+    //Create the preconditioner
+    Teuchos::RCP<precon_t> precon;
+    Ifpack2::Factory preconfactory;
+    auto preconparams = Teuchos::sublist(params,"Preconditioner Parameters");
+    if( handle.precontype!=-1 )
+    {
+      precon = preconfactory.create( supportedPreconNames[handle.precontype], B );
+      precon->setParameters( *preconparams );
+      precon->initialize();
+      precon->compute();
+    }
+
+    ASSERT( B->getRangeMap()== B->getDomainMap() );
+
+    //Compute the transpose of the operator
+    auto BT = B;
+    auto preconT = precon;
+    if( !params->get("Symmetric",false) )
+    {
+      BT = utils::matrix::transpose( B );
+
+      preconT = preconfactory.create( supportedPreconNames[handle.precontype], BT );
+      preconT->setParameters( *preconparams );
+      preconT->initialize();
+      preconT->compute();
+    }
+
+    ASSERT( BT->getRangeMap()== B->getRangeMap() );
+    ASSERT( BT->getRangeMap()== B->getRangeMap() );
+    ASSERT( BT->getRowMap()  == B->getRangeMap() );
+
+  /*+++++++++++++++++++++++
+    + Matrix 1-norm       +
+    +++++++++++++++++++++++*/
+
+    auto norm1B = utils::matrix::norm1( B );
+
+  /*+++++++++++++++++++++++
+    + Hager's algorithm   +
+    +++++++++++++++++++++++*/
+
+    //Construct the initial vectors for
+    auto x  = Teuchos::rcp( new vector_t( B->getRangeMap()  ) );
+    auto y  = Teuchos::rcp( new vector_t( B->getDomainMap() ) );
+    auto z  = Teuchos::rcp( new vector_t( B->getRangeMap()  ) );
+    auto xi = Teuchos::rcp( new vector_t( B->getDomainMap() ) );
+
+    x->putScalar( 1./x->getMap()->getGlobalNumElements() );
+
+    auto linprob = Teuchos::rcp( new linearproblem_t( B , y, x  ) );
+    auto adjprob = Teuchos::rcp( new linearproblem_t( BT, z, xi ) );
+
+    if( handle.precontype != -1 )
+    {
+      linprob->setLeftPrec( precon  );
+      adjprob->setLeftPrec( preconT );
+    }
+
+    if( params->get("Symmetric",false) )
+    {
+      linprob->setHermitian();
+      adjprob->setHermitian();
+    }
+
+    global_t iiter;
+    scalar_t gamma = NAN;
+
+    for( iiter=0 ; iiter<=x->getMap()->getGlobalNumElements(); iiter++ )
+    {
+      y->putScalar( 0. );
+      linprob->setProblem( y, x );
+
+      solvermgr->setProblem( linprob );
+      Belos::ReturnType result = solvermgr->solve();
+      ASSERT( result == Belos::Converged );
+
+      //Compute the sign vector
+      auto xi_vals = xi->getDataNonConst();
+      auto y_vals  = y->getData();
+
+      ASSERT( xi_vals.size()==y_vals.size() );
+
+      xi->putScalar(1.);
+      for( local_t irow=0; irow<xi_vals.size() ; irow++ )
+      {
+        if( y_vals[irow] < 0. ) xi_vals[irow] = -1.;
+      }
+      z->putScalar( 0. );
+      adjprob->setProblem( z, xi );
+
+      solvermgr->setProblem( adjprob );
+      result = solvermgr->solve();
+      ASSERT( result == Belos::Converged );
+
+      if( z->normInf()<=z->dot(*x) )
+      {
+        gamma = y->norm1();
+        break;
+      }
+
+      scalar_t max    = NAN;
+      local_t  argmax = -1;
+      for( local_t irow=0 ; irow<z->getMap()->getNodeNumElements() ; irow++ )
+      {
+        if( std::isnan(max) || std::abs(z->getData()[irow])>max )
+        {
+           max=std::abs(z->getData()[irow]);
+           argmax = irow;
+        }
+      }
+
+      Teuchos::Array<scalar_t> allmax ( nprocs, 0. );
+      Teuchos::Array<local_t> allargmax ( nprocs, -1 );
+      Teuchos::gatherAll(*x->getMap()->getComm(),1,&max,nprocs,allmax.getRawPtr());
+      Teuchos::gatherAll(*x->getMap()->getComm(),1,&argmax,nprocs,allargmax.getRawPtr());
+
+      int maxrank = 0;
+      max = allmax[0];
+      for( int rank = 1 ; rank < nprocs ; rank++ )
+      {
+        if( allmax[rank]>max )
+        {
+          maxrank = rank;
+          max     = allmax[rank];
+        }
+      }
+
+      x->putScalar(0.);
+      if( maxrank==myrank ) x->replaceLocalValue( allargmax[maxrank], 1. );
+    }
+
+    ASSERT( !std::isnan(gamma) );
+
+    out(INFO) << "Hager's algorithm converged in " << iiter << " iterations." << std::endl;
+
+    scalar_t condest = norm1B*gamma;
+
+    gather( &condest );
+  }
+
   void matrix_toarray() /* matrix vector multiplication
 
        -> broadcast HANDLE handle.matrix
@@ -1092,12 +1285,9 @@ class LibMatrix : public Intercomm {
     precon->initialize();
     precon->compute();
 
-    magnitude_t condest = precon->computeCondEst( Ifpack2::Cheap );
-    out(INFO) << "Ifpack2 preconditioner's estimated condition number: " << condest << std::endl;
-
     objects.set( handle.precon, precon, out(DEBUG) );
   }
-  
+
   void export_new() /* create new exporter
 
        -> broadcast HANDLE handle.{exporter,srcmap,dstmap}