diff --git a/Makefile b/Makefile
index ecdaa97..7ebe329 100644
--- a/Makefile
+++ b/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
diff --git a/libmatrix.cpp b/libmatrix.cpp
index bdb19d4..8abb266 100644
--- a/libmatrix.cpp
+++ b/libmatrix.cpp
@@ -1,10 +1,7 @@
-#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>
@@ -12,9 +9,18 @@
 #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}
diff --git a/libmatrix.py b/libmatrix.py
index 48881e6..7b38e13 100644
--- a/libmatrix.py
+++ b/libmatrix.py
@@ -20,6 +20,10 @@
 token_t  = numpy.dtype( _info.pop( 'token_t'  ) )
 char_t   = numpy.dtype( 'str' )
 
+# from BelosTypes.h
+Warnings, IterationDetails, OrthoDetails, FinalSummary, TimingDetails, StatusTestDetails, Debug = 2**numpy.arange(7)
+General, Brief = range(2)
+
 def cacheprop( f ):
   name = f.func_name
   def wrapped( self ):
@@ -113,6 +117,12 @@ def vector_new( self, map_handle ):
     self.bcast( [ vec_handle, map_handle ], handle_t )
     return vec_handle
 
+  @bcast_token
+  def vector_abs( self, vec_handle ):
+    abs_handle = self.claim_handle()
+    self.bcast( [ vec_handle, abs_handle ], handle_t )
+    return abs_handle
+
   @bcast_token
   def vector_copy( self, orig_handle ):
     copy_handle = self.claim_handle()
@@ -159,7 +169,12 @@ def vector_norm( self, handle ):
   @bcast_token
   def vector_sum( self, handle ):
     self.bcast( handle, handle_t )
-    return self.gather( scalar_t ).sum()
+    return self.reduce( scalar_t, numpy.sum )
+
+  @bcast_token
+  def vector_max( self, handle ):
+    self.bcast( handle, handle_t )
+    return self.reduce( scalar_t, numpy.max )
 
   @bcast_token
   def vector_dot( self, handle1, handle2 ):
@@ -211,6 +226,11 @@ def precon_new( self, matrix_handle, precontype_handle, preconparams_handle ):
     self.bcast( [ precon_handle, matrix_handle, precontype_handle, preconparams_handle ], handle_t )
     return precon_handle
 
+  @bcast_token
+  def matrix_condest( self, operator_handle, params_handle, solvername_handle, preconname_handle ):
+    self.bcast( [operator_handle,params_handle,solvername_handle,preconname_handle], handle_t )
+    return self.gather_equal( scalar_t )
+
   @bcast_token
   def matrix_add( self, mat1_handle, mat2_handle, scale_mat1, scale_mat2 ):
     sum_handle = self.claim_handle()
@@ -347,19 +367,23 @@ def __init__( self, comm, items={} ):
   #</ParameterList>
 
   @staticmethod
-  def toxml( items ):
-    s = '<ParameterList>\n'
-    for key, value in items.items():
+  def toxml( items, name='ANONYMOUS' ):
+    s = '<ParameterList name="%s">\n' % name
+    for key,value in items.items():
       assert isinstance( key, str ), 'parameter keys should be strings'
-      if isinstance( value, bool ):  
-        dtype = 'bool'
-      elif isinstance( value, int ):
-        dtype = 'int'
-      elif isinstance( value, float ):
-        dtype = 'double'
+
+      if isinstance(value,dict):
+        s+=ParameterList.toxml(value,name=key)
       else:
-        raise Exception, 'invalid value %r' % value
-      s += '<Parameter isUsed="true" name="%s" type="%s" value="%s"/>\n' % ( key, dtype, value )
+        if isinstance( value, bool ):  
+          dtype = 'bool'
+        elif isinstance( value, int ):
+          dtype = 'int'
+        elif isinstance( value, float ):
+          dtype = 'double'
+        else:
+          raise Exception, 'invalid value %r' % value
+        s += '<Parameter isUsed="true" name="%s" type="%s" value="%s"/>\n' % ( key, dtype, value )
     s += '</ParameterList>'
     return s
 
@@ -479,9 +503,15 @@ def __len__ ( self ):
   def norm( self ):
     return self.comm.vector_norm( self.handle )
 
+  def abs ( self ):
+    return Vector( self.map, self.comm.vector_abs( self.handle ) )
+
   def sum( self ):
     return self.comm.vector_sum( self.handle )
 
+  def max( self ):
+    return self.comm.vector_max( self.handle )
+
   def dot( self, other ):
     assert isinstance( other, Vector )
     assert self.shape == other.shape
@@ -726,6 +756,32 @@ def build_precon( self, precontype, fill=5., absthreshold=0., relthreshold=1., r
     myhandle = self.comm.precon_new( self.handle, _precons.index(precontype), preconparams.handle )
     return Operator( myhandle, self.rangemap, self.domainmap )
 
+  def condest( self, tol, symmetric=False, precon=None, maxiter=1000, outfreq=10, fill=5., absthreshold=0., relthreshold=1., relax=0. ):
+    solverparams = {
+      'Verbosity': Warnings | FinalSummary,
+      'Maximum Iterations': maxiter,
+      'Output Style': Brief,
+      'Convergence Tolerance': tol,
+      'Output Frequency': outfreq,
+    }
+    params = ParameterList( self.comm, {
+      'Symmetric': symmetric,
+      'Solver Parameters': solverparams,
+    })
+    if precon:
+      params['Preconditioner Parameters'] = {
+        'fact: %s level-of-fill' % precon.lower(): float(fill),
+        'fact: absolute threshold': float(absthreshold),
+        'fact: relative threshold': float(relthreshold),
+        'fact: relax value': float(relax),
+      }
+      precon_handle = _precons.index( precon )
+    else:  
+      precon_handle = numpy.array(-1,dtype=handle_t)
+      
+    solver_handle = _solvers.index( 'CG' if symmetric else 'GMRES' )
+    return self.comm.matrix_condest( self.handle, params.handle, solver_handle, precon_handle )
+
 
 class MatrixBuilder( Object ):
 
@@ -801,9 +857,6 @@ def set_precon( self, precon, right=False ):
     self.comm.linearproblem_set_precon( self.handle, precon.handle, right )
 
   def solve( self, name, tol, maxiter=1000, outfreq=10 ):
-    # from BelosTypes.h
-    Warnings, IterationDetails, OrthoDetails, FinalSummary, TimingDetails, StatusTestDetails, Debug = 2**numpy.arange(7)
-    General, Brief = range(2)
     solverparams = ParameterList( self.comm, {
       'Verbosity': StatusTestDetails | FinalSummary,
       'Maximum Iterations': maxiter,
@@ -873,5 +926,6 @@ def ArrayBuilder( shape ):
   if len( shape ) == 1:
     return VectorBuilder( shape[0] )
 
+  raise NotImplementedError()  
 
 # vim:shiftwidth=2:foldmethod=indent:foldnestmax=2
diff --git a/mpi.py b/mpi.py
index 8e7c3bb..d520806 100644
--- a/mpi.py
+++ b/mpi.py
@@ -1,7 +1,6 @@
 import numpy
 from mpi4py import MPI
 
-
 class InterComm( object ):
   'generic MPI communicator wrapper'
 
@@ -23,11 +22,13 @@ def gather( self, dtype ):
     array = numpy.empty( self.nprocs, dtype=dtype )
     self.__comm.Gather( None, [ array, MPI.BYTE ], root=MPI.ROOT )
     return array
-  
+
   def gather_equal( self, dtype ):
+    return self.reduce( dtype, assert_equal )
+
+  def reduce( self, dtype, op ):
     array = self.gather( dtype )
-    assert numpy.all( array[1:] == array[0] )
-    return array[0]
+    return op(array)
   
   def scatter( self, array, dtype ):
     array = numpy.ascontiguousarray( array, dtype=dtype )
@@ -72,3 +73,11 @@ def disconnect( self ):
 
   def __del__( self ):
     self.disconnect()
+
+###############################
+# Utility functions           #
+###############################
+
+def assert_equal( array ):
+  assert numpy.all( array[1:] == array[0] )
+  return array[0]
diff --git a/test.py b/test.py
index a622c22..f001278 100755
--- a/test.py
+++ b/test.py
@@ -66,7 +66,7 @@ def create_overlapping_rowmap():
 @withcomm()
 def distributed_vector(comm):
   '''
-  Create a uncompletely filled vector "v" with entries i^2 if i is odd and 0 otherwise
+  Create a uncompletely filled vector "v" with entries 4*i*(n-i-1)/(n-1)**2-0.5 if i is even and 0 otherwise
   and a vector "w" with all entries equal to i
   '''
 
@@ -77,12 +77,13 @@ def distributed_vector(comm):
   rowmap = libmatrix.Map( comm, rowdofmap, ndofs )
 
   v   = libmatrix.VectorBuilder( rowmap )
-  for idx in range(1,ndofs+1,2):
-    v.add_global( idx=[numpy.array([idx])], data=[numpy.array([idx**2],dtype=float)] )
+  for idx in range(0,ndofs,2):
+    v.add_global( idx=[numpy.array([idx])], data=[numpy.array([(4./(ndofs-1)**2)*idx*(ndofs-1-idx)-0.5],dtype=float)] )
   v = v.complete()
-  
-  v_npy = numpy.arange(0,ndofs,dtype=float)**2
-  v_npy[0::2] = 0. 
+ 
+  v_npy = numpy.arange(0,ndofs,dtype=float)
+  v_npy = (4./(ndofs-1.)**2)*v_npy*(ndofs-1-v_npy)-0.5
+  v_npy[1::2] = 0. 
 
   numpy.testing.assert_almost_equal( v.toarray(), v_npy, err_msg='Vector "v" was not filled correctly' )
 
@@ -97,6 +98,7 @@ def distributed_vector(comm):
   numpy.testing.assert_almost_equal( (v*2.).toarray(), 2.*v_npy, err_msg='Scalar multiplication not computed correctly' )
   numpy.testing.assert_almost_equal( (v/2.).toarray(), v_npy/2., err_msg='Scalar division not computed correctly' )
   numpy.testing.assert_almost_equal( (2./(v+1.)).toarray(), 2./(v_npy+1.), err_msg='Scalar division not computed correctly' )
+  numpy.testing.assert_almost_equal( v.abs().toarray(), numpy.abs(v_npy), err_msg='Vector "v" absolute value failed' )
 
   #Create the vector w based on the reversed rowdofmap
   w = libmatrix.VectorBuilder( rowmap )
@@ -341,6 +343,13 @@ def solve_laplace(comm):
   x = Ac.solve( rhs=vc, precon=None, symmetric=True, tol=1e-10 )
   numpy.testing.assert_almost_equal( x.toarray(), x_npy )
 
+  #Maximum operation
+  numpy.testing.assert_almost_equal( x.max(), x_npy.max() )
+
+  #Compute the condition number estimate
+  condest = Ac.condest( tol=1e-10, symmetric=True, precon='ILUT' )
+  numpy.testing.assert_almost_equal( condest, numpy.linalg.cond( Ac_npy, p=1 ) )
+
 ## END UNIT TESTS ##
 
 unittest.exit()
diff --git a/typedefs.h b/typedefs.h
new file mode 100644
index 0000000..976b9fc
--- /dev/null
+++ b/typedefs.h
@@ -0,0 +1,10 @@
+#ifndef LIBMATRIX_TYPEDEFS_H
+#define LIBMATRIX_TYPEDEFS_H
+
+typedef double scalar_t;
+typedef int number_t;
+typedef int local_t;
+typedef long global_t;
+typedef bool bool_t;
+
+#endif
diff --git a/utils.cpp b/utils.cpp
new file mode 100644
index 0000000..ec0ffda
--- /dev/null
+++ b/utils.cpp
@@ -0,0 +1,61 @@
+#include "utils.h"
+
+#include <Tpetra_RowMatrixTransposer.hpp>
+
+/*
+Vector utilities
+*/
+
+Teuchos::RCP<vector_t> utils::vector::abs ( Teuchos::RCP<vector_t> vector )
+{
+  auto absvector = Teuchos::rcp(new vector_t( vector->getMap() ));
+  auto in = vector->getData();
+  auto out = absvector->getDataNonConst();
+  for( local_t irow=0 ; irow<absvector->getMap()->getNodeNumElements() ; irow++ )
+  {
+    out[irow] = std::abs(in[irow]);
+  }
+  return absvector;
+}
+
+
+/*
+Matrix utilities
+*/
+
+Teuchos::RCP<crsmatrix_t> utils::matrix::transpose ( const Teuchos::RCP<const crsmatrix_t> matrix )
+{
+  Tpetra::RowMatrixTransposer<scalar_t,local_t,global_t,node_t> transposer ( matrix );
+  return transposer.createTranspose();
+}
+
+
+scalar_t utils::matrix::norminf ( const Teuchos::RCP<const crsmatrix_t> matrix )
+{
+  scalar_t norm = 0.;
+  scalar_t rowsum;
+  Teuchos::ArrayView<const local_t> indices;
+  Teuchos::ArrayView<const scalar_t> values;
+
+  for( local_t irow=0 ; irow<matrix->getRowMap()->getNodeNumElements() ; irow++ )
+  {
+    matrix->getLocalRowView(irow, indices, values);
+    rowsum = 0.;
+    for( auto &v : values ) rowsum += std::abs(v);
+    if( rowsum > norm ) norm=rowsum;
+  }
+
+  if (matrix->isDistributed())
+  {
+    scalar_t lnorm = norm;
+    Teuchos::reduceAll(*matrix->getRowMap()->getComm(),Teuchos::REDUCE_MAX,lnorm,Teuchos::outArg(norm));
+  }
+
+  return norm;
+}
+
+
+scalar_t utils::matrix::norm1 ( const Teuchos::RCP<const crsmatrix_t> matrix )
+{
+  return utils::matrix::norminf( utils::matrix::transpose(matrix) );
+}
diff --git a/utils.h b/utils.h
new file mode 100644
index 0000000..b543b29
--- /dev/null
+++ b/utils.h
@@ -0,0 +1,32 @@
+#ifndef LIBMATRIX_UTILS_H
+#define LIBMATRIX_UTILS_H
+
+#include "typedefs.h"
+
+#include <Teuchos_RCP.hpp>
+#include <Tpetra_CrsMatrix.hpp>
+
+typedef Kokkos::DefaultNode::DefaultNodeType node_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;
+
+namespace utils
+{
+  namespace vector
+  {
+    Teuchos::RCP<vector_t> abs ( Teuchos::RCP<vector_t> );
+
+    //local_t argmax ( )
+  }
+
+  namespace matrix
+  {
+    Teuchos::RCP<crsmatrix_t> transpose ( const Teuchos::RCP<const crsmatrix_t>  );
+
+    scalar_t norminf ( const Teuchos::RCP<const crsmatrix_t> );
+
+    scalar_t norm1 ( const Teuchos::RCP<const crsmatrix_t> );
+  }
+}
+
+#endif