diff --git a/src/LinAlg/DenseMatrix.h b/src/LinAlg/DenseMatrix.h
index f8566f22..2c979123 100644
--- a/src/LinAlg/DenseMatrix.h
+++ b/src/LinAlg/DenseMatrix.h
@@ -138,6 +138,9 @@ public:
   bool solveEig(DenseMatrix& B, RealArray& eigVal, Matrix& eigVec, int nev,
 		real shift = 0.0);
 
+  //! \brief Returns the L-infinity norm of the matrix.
+  virtual real Linfnorm() const { return myMat.normInf(); }
+
 protected:
   //! \brief Augments a dense matrix symmetrically to the current matrix.
   //! \param[in] B  The matrix to be augmented
diff --git a/src/LinAlg/PETScMatrix.C b/src/LinAlg/PETScMatrix.C
index f6196936..df4729a8 100644
--- a/src/LinAlg/PETScMatrix.C
+++ b/src/LinAlg/PETScMatrix.C
@@ -15,6 +15,9 @@
 #ifdef HAS_PETSC
 #include "SAM.h"
 #include "petscmg.h"
+#ifdef HAS_SLEPC
+#include "slepceps.h"
+#endif
 
 
 PETScVector::PETScVector()
@@ -151,20 +154,17 @@ real PETScVector::Linfnorm() const
 
 PETScMatrix::PETScMatrix(const LinSolParams& spar) : solParams(spar)
 {
-  // Create matrix object. By default the matrix type is AIJ
+  // Create matrix object, by default the matrix type is AIJ
   MatCreate(PETSC_COMM_WORLD,&A);
 
-  // Create linear solver object.
-  KSPCreate(PETSC_COMM_WORLD,&ksp); 
+  // Create linear solver object
+  KSPCreate(PETSC_COMM_WORLD,&ksp);
 
   // Create null space if any
   if (solParams.getNullSpace() == CONSTANT) {
     MatNullSpaceCreate(PETSC_COMM_WORLD,PETSC_TRUE,0,0,&nsp);
     KSPSetNullSpace(ksp,nsp);
   }
-
-  // Create eigenvalue solver object.
-  //EPSCreate(PETSC_COMM_WORLD,&eps);
 }
 
 
@@ -182,17 +182,11 @@ PETScMatrix::PETScMatrix (const PETScMatrix& B) : solParams(B.solParams)
     MatNullSpaceCreate(PETSC_COMM_WORLD,PETSC_TRUE,0,0,&nsp);
     KSPSetNullSpace(ksp,nsp);
   }
-
-  // Create eigenvalue solver object.
-  //EPSCreate(PETSC_COMM_WORLD,&eps);
 }
 
 
 PETScMatrix::~PETScMatrix ()
 {
-  // Deallocation of eigenvalue solver object.
-  //EPSDestroy(eps);
-
   // Deallocation of null space
   if (solParams.getNullSpace() == CONSTANT) 
     MatNullSpaceDestroy(nsp);
@@ -205,20 +199,21 @@ PETScMatrix::~PETScMatrix ()
 }
 
 
-static void assemPETScPara(const Matrix& eM, Mat SM, PETScVector& SV,
-			const std::vector<int>& meen, const int* meqn,
-			const int* mpmceq, const int* mmceq, const real* ttcc)
+static void assemPETScPara (const Matrix& eM, Mat SM, PETScVector& SV,
+			    const std::vector<int>& meen, const int* meqn,
+			    const int* mpmceq, const int* mmceq,
+			    const real* ttcc)
 {
-  real c0;
-  int i, j, jp, jceq;
+  real   c0;
+  size_t i, j;
+  int    jp, jceq;
 
   // Number of degrees of freedom for element
-  const int nedof = meen.size();
+  size_t nedof = meen.size();
 
-  // Convert meen to C array. Start column numbering at 0.
-  int *l2g;
-  l2g = new int[nedof];
-  for (i = 0;i < nedof;i++)
+  // Convert meen to 0-based C array
+  int* l2g = new int[nedof];
+  for (i = 0; i < nedof; i++)
     l2g[i] = meqn[meen[i]-1]-1;
 
   // Cast to non-constant Matrix to modify for Dirichlet BCs
@@ -270,24 +265,24 @@ static void assemPETScPara(const Matrix& eM, Mat SM, PETScVector& SV,
 }
 
 
-static void assemPETSc(const Matrix& eM, Mat SM, PETScVector& SV,
-		       const std::vector<int>& meen, const int* meqn,
-		       const int* mpmceq, const int* mmceq, const real* ttcc)
+static void assemPETSc (const Matrix& eM, Mat SM, PETScVector& SV,
+			const std::vector<int>& meen, const int* meqn,
+			const int* mpmceq, const int* mmceq, const real* ttcc)
 {
-  real c0;
-  int  i, j, ieq, jeq, ip, jp, iceq, jceq;
+  real   c0;
+  size_t i, j;
+  int    ieq, jeq, ip, jp, iceq, jceq;
 
   // Get C array
   PetscScalar* svec;
   VecGetArray(SV.getVector(),&svec); 
   
   // Number of degrees of freedom for element
-  const int nedof = meen.size();
+  size_t nedof = meen.size();
 
-  // Convert meen to C array. Start column numbering at 0. 
-  int *l2g;
-  l2g = new int[nedof];
-  for (i = 0;i < nedof; i++)
+  // Convert meen to 0-based C array
+  int* l2g = new int[nedof];
+  for (i = 0; i < nedof; i++)
     l2g[i] = meen[i]-1;
 
   // Cast to non-constant Matrix 
@@ -358,20 +353,20 @@ static void assemPETSc(const Matrix& eM, Mat SM, PETScVector& SV,
 }
 
 
-static void assemPETSc(const Matrix& eM, Mat SM, const std::vector<int>& meen, 
-		       const int* meqn, const int* mpmceq, const int* mmceq, 
-		       const real* ttcc)
+static void assemPETSc (const Matrix& eM, Mat SM, const std::vector<int>& meen,
+			const int* meqn, const int* mpmceq, const int* mmceq,
+			const real* ttcc)
 {
-  real c0;
-  int  i, j, ieq, jeq, ip, jp, iceq, jceq;
+  real   c0;
+  size_t i, j;
+  int    ieq, jeq, ip, jp, iceq, jceq;
 
   // Number of degrees of freedom for element
-  const int nedof = meen.size();
+  size_t nedof = meen.size();
 
-  // Convert meen to C array. Start column numbering at 0. 
-  int *l2g;
-  l2g = new int[nedof];
-  for (i = 0;i < nedof; i++) 
+  // Convert meen to 0-based C array
+  int* l2g = new int[nedof];
+  for (i = 0; i < nedof; i++)
     l2g[i] = meen[i]-1;
 
   // Cast to non-constant Matrix 
@@ -559,12 +554,10 @@ bool PETScMatrix::solve (SystemVector& B, bool newLHS)
 }
 
 
-bool PETScMatrix::solveEig(PETScMatrix& B, RealArray& val, 
-			   Matrix& vec, int nv, real shift, int iop)
+bool PETScMatrix::solveEig (PETScMatrix& B, RealArray& val,
+			    Matrix& vec, int nv, real shift, int iop)
 {
 #ifdef HAS_SLEPC
-  return false;
-  int         i;
   ST          st;
   PetscInt    m, n, nconv;
   PetscScalar kr, ki;
@@ -600,7 +593,7 @@ bool PETScMatrix::solveEig(PETScMatrix& B, RealArray& val,
 
   val.resize(nv);
   vec.resize(n,nv);
-  for ( i = 0; i < nv; i++) {
+  for (int i = 0; i < nv; i++) {
     EPSGetEigenpair(eps,i,&kr,&ki,xr,xi);
     VecGetArray(xr,&xrarr);
     
@@ -620,4 +613,12 @@ bool PETScMatrix::solveEig(PETScMatrix& B, RealArray& val,
   return false;
 }
 
+
+real PETScMatrix::Linfnorm () const
+{
+  PetscReal norm;
+  MatNorm(A,NORM_INFINITY,&norm);
+  return norm;
+}
+
 #endif // HAS_PETSC
diff --git a/src/LinAlg/PETScMatrix.h b/src/LinAlg/PETScMatrix.h
index b3dc634f..ba729df7 100644
--- a/src/LinAlg/PETScMatrix.h
+++ b/src/LinAlg/PETScMatrix.h
@@ -97,7 +97,6 @@ private:
   Vec x; //!< The actual PETSc vector.
 
 #else // dummy implementation when PETSc is not included
-  PETScVector() {}
   virtual SystemVector* copy() const { return 0; }
   virtual size_t dim() const { return 0; }
   virtual void redim(size_t) {}
@@ -203,6 +202,9 @@ public:
   bool solveEig(PETScMatrix& B, RealArray& eigVal, Matrix& eigVec, int nev,
 		real shift = real(0), int iop = 1);
 
+  //! \brief Returns the L-infinity norm of the matrix.
+  virtual real Linfnorm() const;
+
 private:
   Mat                 A;         //!< Linear system matrix
   KSP                 ksp;       //!< Linear solver
diff --git a/src/LinAlg/SPRMatrix.C b/src/LinAlg/SPRMatrix.C
index 233ac59a..87e1b1ca 100644
--- a/src/LinAlg/SPRMatrix.C
+++ b/src/LinAlg/SPRMatrix.C
@@ -403,3 +403,10 @@ bool SPRMatrix::solveEig (SPRMatrix& B, RealArray& val, Matrix& vec, int nv,
 #endif
   return false;
 }
+
+
+real SPRMatrix::Linfnorm () const
+{
+  std::cerr <<"SPRMatrix::Linfnorm not yet implemented"<< std::endl;
+  return real(0);
+}
diff --git a/src/LinAlg/SPRMatrix.h b/src/LinAlg/SPRMatrix.h
index 246e977f..280a262b 100644
--- a/src/LinAlg/SPRMatrix.h
+++ b/src/LinAlg/SPRMatrix.h
@@ -106,6 +106,9 @@ public:
   bool solveEig (SPRMatrix& B, RealArray& eigVal, Matrix& eigVec, int nev,
 		 real shift = 0.0, int iop = 1);
 
+  //! \brief Returns the L-infinity norm of the matrix.
+  virtual real Linfnorm() const;
+
 private:
   int mpar[NS]; //!< Matrix of sparse PARameters
   int* msica;   //!< Matrix of Storage Information for CA
diff --git a/src/LinAlg/SparseMatrix.C b/src/LinAlg/SparseMatrix.C
index f82b6edb..127dddcd 100644
--- a/src/LinAlg/SparseMatrix.C
+++ b/src/LinAlg/SparseMatrix.C
@@ -199,7 +199,7 @@ real& SparseMatrix::operator () (size_t r, size_t c)
 	      << nrow <<"x"<< ncol << std::endl;
   else if (editable) {
     IJPair index(r,c);
-    if (elem.find(index) == elem.end()) elem[index] = 0.0;
+    if (elem.find(index) == elem.end()) elem[index] = real(0);
     return elem[index];
   }
 
@@ -988,3 +988,25 @@ bool SparseMatrix::solveSAMG (bool isFirstRHS, Vector& B)
 #endif
   return ierr <= 0;
 }
+
+
+real SparseMatrix::Linfnorm () const
+{
+  RealArray sums(nrow,real(0));
+
+  if (editable)
+    for (ValueIter it = elem.begin(); it != elem.end(); ++it)
+      sums[it->first.first-1] += fabs(it->second);
+  else if (solver == SUPERLU)
+    // Column-oriented format with 0-based row-indices
+    for (size_t j = 1; j <= ncol; j++)
+      for (int i = IA[j-1]; i < IA[j]; i++)
+	sums[JA[i]] += fabs(A[i]);
+  else
+    // Row-oriented format with 1-based row-indices
+    for (size_t i = 1; i <= nrow; i++)
+      for (int j = IA[i-1]; j < IA[i]; j++)
+	sums[i-1] += fabs(A[j-1]);
+
+  return *std::max_element(sums.begin(),sums.end());
+}
diff --git a/src/LinAlg/SparseMatrix.h b/src/LinAlg/SparseMatrix.h
index 373f3331..6341cc9b 100644
--- a/src/LinAlg/SparseMatrix.h
+++ b/src/LinAlg/SparseMatrix.h
@@ -196,6 +196,9 @@ protected:
   //! \brief Writes the system matrix to the given output stream.
   virtual std::ostream& write(std::ostream& os) const;
 
+  //! \brief Returns the L-infinity norm of the matrix.
+  virtual real Linfnorm() const;
+
 public:
   static bool printSLUstat; //!< Print solution statistics for SuperLU?
 
diff --git a/src/LinAlg/SystemMatrix.h b/src/LinAlg/SystemMatrix.h
index 25adb6cc..a4dea93e 100644
--- a/src/LinAlg/SystemMatrix.h
+++ b/src/LinAlg/SystemMatrix.h
@@ -265,6 +265,9 @@ public:
   //! \param newLHS \e true if the left-hand-side matrix has been updated
   virtual bool solve(SystemVector&, bool newLHS = true) { return false; }
 
+  //! \brief Returns the L-infinity norm of the matrix.
+  virtual real Linfnorm() const = 0;
+
 protected:
   //! \brief Writes the system matrix to the given output stream.
   virtual std::ostream& write(std::ostream& os) const { return os; }
diff --git a/src/LinAlg/matrix.h b/src/LinAlg/matrix.h
index eb80c4e8..d23ef6c2 100644
--- a/src/LinAlg/matrix.h
+++ b/src/LinAlg/matrix.h
@@ -505,6 +505,18 @@ namespace utl //! General utility classes and functions.
       return true;
     }
 
+    //! \brief Return the infinite norm of the matrix.
+    T normInf() const
+    {
+      if (nrow == 0) return T(0);
+
+      // Compute row sums
+      vector<T> sums(nrow);
+      for (size_t i = 0; i < nrow; i++)
+	sums[i] = elem.asum(i,ncol);
+      return *std::max_element(sums.begin(),sums.end());
+    }
+
   private:
     //! \brief Check dimension compatibility for matrix-vector multiplication.
     bool compatible(const std::vector<T>& X, const std::vector<T>& Y,