diff --git a/src/LinAlg/SparseMatrix.C b/src/LinAlg/SparseMatrix.C
index 36a4c6f2..e139ef48 100644
--- a/src/LinAlg/SparseMatrix.C
+++ b/src/LinAlg/SparseMatrix.C
@@ -819,7 +819,7 @@ bool SparseMatrix::solveSLU (Vector& B)
 			   &A.front(), &JA.front(), &IA.front(),
 			   SLU_NC, SLU_D, SLU_GE);
   }
-  else if (factored) 
+  else if (factored)
     slu->opts->Fact = FACTORED; // Re-use previous factorization
   else
  {
diff --git a/src/LinAlg/SystemMatrix.C b/src/LinAlg/SystemMatrix.C
index 308a800f..9b352657 100644
--- a/src/LinAlg/SystemMatrix.C
+++ b/src/LinAlg/SystemMatrix.C
@@ -37,28 +37,24 @@ SystemVector* SystemVector::create (Type vectorType)
 }
 
 
-void SystemVector::copy(const SystemVector& x)
+SystemVector& SystemVector::copy (const SystemVector& x)
 {
   this->redim(x.size());
-  SystemVector* xsv = const_cast<SystemVector*>(&x);
   real* vec = this->getPtr();
-  memcpy(vec,xsv->getPtr(),x.dim()*sizeof(real));
+  memcpy(vec,x.getRef(),x.dim()*sizeof(real));
   this->restore(vec);
+
+  return *this;
 }
 
 
 SystemMatrix* SystemMatrix::create (Type matrixType, const LinSolParams& spar)
 {
-  if (matrixType == PETSC)
 #ifdef HAS_PETSC
+  if (matrixType == PETSC)
     return new PETScMatrix(spar);
-#else
-  {
-    std::cerr << "PETSc support not compiled in, bailing" << std::endl;
-    exit(1);
-  }
 #endif
- 
+
   return SystemMatrix::create(matrixType);
 }
 
@@ -77,7 +73,8 @@ SystemMatrix* SystemMatrix::create (Type matrixType, int num_thread_SLU)
       static LinSolParams defaultPar;
       return new PETScMatrix(defaultPar);
 #else
-      std::cerr << "PETSc support not compiled in, bailing" << std::endl;
+      std::cerr <<"SystemMatrix::create: PETSc not compiled in, bailing out..."
+		<< std::endl;
       exit(1);
 #endif
     default:
@@ -87,10 +84,3 @@ SystemMatrix* SystemMatrix::create (Type matrixType, int num_thread_SLU)
 
   return 0;
 }
-
-
-bool SystemMatrix::solve (const SystemVector& b, SystemVector& x, bool newLHS)
-{
-  x.copy(b);
-  return this->solve(x,newLHS);
-}
diff --git a/src/LinAlg/SystemMatrix.h b/src/LinAlg/SystemMatrix.h
index b49511ee..a2ab83a4 100644
--- a/src/LinAlg/SystemMatrix.h
+++ b/src/LinAlg/SystemMatrix.h
@@ -76,8 +76,8 @@ public:
   //! \brief Initializes the vector assuming it is properly dimensioned.
   virtual void init(real value = real(0)) = 0;
 
-  //! \brief Copies entries from input vector
-  virtual void copy(const SystemVector& x);
+  //! \brief Copies entries from input vector into \a *this.
+  SystemVector& copy(const SystemVector& x);
 
   //! \brief Begins communication step needed in parallel vector assembly.
   virtual bool beginAssembly() { return true; }
@@ -265,21 +265,26 @@ public:
   virtual bool multiply(const SystemVector&, SystemVector&) { return false; }
 
   //! \brief Solves the linear system of equations for a given right-hand-side.
-  //! \param newLHS \e true if the left-hand-side matrix has been updated
-  virtual bool solve(SystemVector&, bool newLHS = true) { return false; }
+  //! \param b Right-hand-side vector on input, solution vector on output
+  //! \param[in] newLHS \e true if the left-hand-side matrix has been updated
+  virtual bool solve(SystemVector& b, bool newLHS = true) { return false; }
 
   //! \brief Solves the linear system of equations for a given right-hand-side.
-  //! \param b Right-hand-side vector
-  //! \param x Solution vector
-  //! \param newLHS \e true if the left-hand-side matrix has been updated
-  virtual bool solve(const SystemVector& b, SystemVector& x, bool newLHS = true);
-
+  //! \param[in] b Right-hand-side vector
+  //! \param[out] x Solution vector
+  //! \param[in] newLHS \e true if the left-hand-side matrix has been updated
+  virtual bool solve(const SystemVector& b, SystemVector& x, bool newLHS = true)
+  {
+    return this->solve(x.copy(b),newLHS);
+  }
   //! \brief Solves the linear system of equations for a given right-hand-side.
-  //! \param b Right-hand-side vector
+  //! \param b Right-hand-side vector on input, solution vector on output
   //! \param P Preconditioning matrix (if different than system matrix)
-  //! \param newLHS \e true if the left-hand-side matrix has been updated
+  //! \param[in] newLHS \e true if the left-hand-side matrix has been updated
   virtual bool solve(SystemVector& b, SystemMatrix& P, bool newLHS = true)
-  { return false; }
+  {
+    return false;
+  }
 
   //! \brief Returns the L-infinity norm of the matrix.
   virtual real Linfnorm() const = 0;