--- a/include/bvp.hpp
+++ b/include/bvp.hpp
@@ -1,5 +1,9 @@
 /*! \file bvp.hpp
  * \brief Classes domain, BVP, and linear_BVP2 are declared here.
+ *
+ *  \file bvp.cpp
+ *  \brief Implementations and instantiations of the functions and
+ *  classes declared in bvp.hpp
  */
 
 #ifndef __BVP_HPP__
@@ -24,8 +28,9 @@
   using linalg::matrix;
   using boost::shared_ptr;
 
-  /// A domain is a bunch of interior points, a bunch of boundary points,
-  /// and a normal vector at each boundary point.
+  /*! \brief A domain is a bunch of interior points, a bunch of boundary points,
+   *  and a normal vector at each boundary point.
+   */
   class domain{
   public:
     ///Allocate empty domain of given dimension.
@@ -98,10 +103,15 @@
     ///Can't create a domain without at least specifying its
     ///dimension. 
     domain();
+    ///Domain's dimension
     size_t dim;
+    /*! \name The domain's ingredients
+     */
+    //@{
     set<point> interior;
     set<point> boundary;
     map<point, vector> normals;
+    //@}
   };
 
   /*! \brief A boundary value problem 

--- a/include/ddm.hpp
+++ b/include/ddm.hpp
@@ -1,3 +1,9 @@
+/*! \file ddm.hpp
+ *  \brief Domain decomposition methods declared here.
+ *  \file ddm.cpp
+ *  \brief Implementations and instantiations of the functions and
+ *  classes declared in ddm.hpp
+ */
 #ifndef __DDM_HPP__
 #define __DDM_HPP__
 
@@ -15,12 +21,15 @@
   using std::set;
   using boost::shared_ptr;
 
-  //Boundary differential operators that may be specified in peculiar
-  //ways depending on the ddm we're working on.
-  //
-  //On interior boundary points as specified by the domain
-  //decomposition, the alternate Bprime operator is applied instead of the
-  //B operator.
+  /*! \brief A bdry_diff_op for domain decomposition methods.
+   *
+   * Boundary differential operators that may be specified in peculiar
+   * ways depending on the ddm we're working on.
+   *
+   * On interior boundary points as specified by the domain
+   * decomposition, the alternate Bprime operator is applied instead
+   * of the B operator.
+   */
   class ddm_bdry_diff_op : public bdry_diff_op{
   public:
     ddm_bdry_diff_op(shared_ptr<const bdry_diff_op> B_in, 
@@ -36,14 +45,15 @@
 
   //************************ ddm ********************************************
 
-  //A generic domain decomposition method
+  //FIXME: derive this from realfunc
+  ///A generic domain decomposition method
   class ddm{
   public:
-    //What are the domains of this ddm, and what's the bvp?
+    ///What are the domains of this ddm, and what's the bvp?
     ddm(const set<shared_ptr<const domain> >& ds, 
 	shared_ptr<const BVP> thebvp);
         
-    //Relative tolerance. Defaults to 1e-5.
+    ///Relative tolerance. Defaults to 1e-5.
     void set_tolerance(double tol);
     virtual double at(const point& p) const = 0;
     double operator()(const point& p) const;
@@ -51,11 +61,14 @@
     virtual ~ddm();
 
   protected:
+    /// The BVP for this ddm
     shared_ptr<const BVP> bvp;
+    /// The domain decomposition
     set<shared_ptr<const domain> > domains;
+    /// The tolerance until convergence
     double tolerance;
 
-    //Each specific ddm has its own way of doing things.
+    ///Each specific ddm has its own way of doing things.
     virtual void solve() = 0;
   };
 
@@ -63,7 +76,10 @@
   
   //***************** additive_schwarz_ddm ***********************************
 
-  //This method described in Li and Chen [2003].
+  /*! \brief Additive Sschwarz domain decomposition method.
+   *
+   *  This method described in Li and Chen [2003].
+   */
   template <typename RBF>
   class additive_schwarz_ddm : public ddm{
   public:
@@ -71,73 +87,113 @@
 			 shared_ptr<const linear_BVP2> bvp);
     double at(const point& p) const;
   private:
+    /// Each domain has an associated interpolator.n
     map<shared_ptr<const overlapping_domain>, 
 	shared_ptr<interpolator<RBF> > > phis;
 
+    /*! \brief Say which interpolators are at work at given point.
+     *
+     * \param p - The point
+     * \return A set of pointers to the interpolators  governing this point.
+     */
     set<shared_ptr<const interpolator<RBF> > > 
     which_interps(const point& p) const;
 
+    /// Actual algorithm implemented here, called by constructor.
     void solve();
 
-    //Evaluate using the average of relevant interpolators.
+    ///Evaluate using the average of relevant interpolators.
     double avg_interp(set<shared_ptr<const interpolator<RBF> > > 
 		      relevant_interpolators,
 		      const point& p) const;
 
 
+    /*! \brief Return a vector of the current estimated evaluated at all
+     * artificial boundary points.
+    */
     vector at_all_points() const;
     
   };
 
-  //An overlapping domain's constructor, besides needing the same
-  //arguments as required for an ordinary domain, also needs to know
-  //which domains it will overlap with, and on which domain its
-  //boundary points lie. These two extra pieces of information are
-  //given by
-  // 
-  //        set<const shared_ptr<overlapping_domain> >
-  //        map<point, shared_ptr<const overlapping_domain> >
-  // 
-  //respectively. By default, if these arguments are omitted, it is
-  //assumed that the overlapping domain doesn't overlap with anything
-  //at all. If a point in the domain's boundary doesn't overlap with
-  //another domain, then it shouldn't be specified in the map<point,
-  //overlapping_domain> data structure (or it can quite harmlessly be
-  //assigned  to itself).
-  //
-  //Alternatively, the set_overlapper_info(...) function can handle
-  //all of the details.
+  /*! \brief Overlapping domains for the Scwharz DDMs.
+   *
+   * An overlapping domain's constructor, besides needing the same
+   * arguments as required for an ordinary domain, also needs to know
+   * which domains it will overlap with, and on which domain its
+   * boundary points lie. These two extra pieces of information are
+   * given by
+   *  
+   *         set<const shared_ptr<overlapping_domain> >
+   *         map<point, shared_ptr<const overlapping_domain> >
+   *  
+   * respectively. By default, if these arguments are omitted, it is
+   * assumed that the overlapping domain doesn't overlap with anything
+   * at all. If a point in the domain's boundary doesn't overlap with
+   * another domain, then it shouldn't be specified in the map<point,
+   * overlapping_domain> data structure (or it can quite harmlessly be
+   * assigned  to itself).
+   * 
+   * Alternatively, the set_overlapper_info(...) function can handle
+   * all of the details.
+  */
   class overlapping_domain : public domain{
   public:
+    /// A domain that in fact doesn't overlap with anything at all.
     overlapping_domain(string intr, string bdry, string ns);
+    /*! \brief Specify with what other domains it overlaps.
+     *
+     * An overlapping domain's constructor, besides needing the same
+     * arguments as required for an ordinary domain, also needs to know
+     * which domains it will overlap with, and on which domain its
+     * boundary points lie. These two extra pieces of information are
+     * given by
+     *  
+     *         set<const shared_ptr<overlapping_domain> >
+     *         map<point, shared_ptr<const overlapping_domain> >
+     *  
+     * respectively. 
+     */
     overlapping_domain(string intr, string bdry, string ns, 
 		       const set<shared_ptr<const overlapping_domain> >& ols,
 		       const map<point, 
 		       shared_ptr<const overlapping_domain> >& bdry_asst);
 
+    ///Allocate an empty non-overlapping domain of a given dimension
     overlapping_domain(size_t dimension);
+    ///Same, but specifying interior, boundary, and normals.
     overlapping_domain(size_t dimension, set<point> intr, 
 	   set<point> bdry, map<point, vector> ns);
       
+    /// Give the domains with which this domain overlaps.
     set<shared_ptr<const overlapping_domain> > get_domains() const;
     
-    //With what domain this point overlap? If none, return 0.
+    /// With what domain this point overlap? If none, return 0.
     shared_ptr<const overlapping_domain> which_domain(const point& p) const; 
     
-    //Once domains are defined via ordinary domain constructors, this
-    //function will grab a set of domains and say which domain
-    //overlaps with which one. 
     friend void set_overlapper_info(set<shared_ptr<overlapping_domain> > 
 				    domains);
 
-    //Set overlapper info manually.
+    ///Set overlapper info manually.
     void set_overlapper_info(const point& p, 
 			     const shared_ptr<overlapping_domain> o);
 
   private:
+    /// The domains with which this one overlaps.
     set<shared_ptr<const overlapping_domain> > overlappers;
+    //FIXME: make this a multimap
+    /// The domains with which this point overlaps.
     map<point, shared_ptr<const overlapping_domain> > boundary_assignments;
   };
+
+  /*! \brief Automate the initialisation of setting domain
+   * information.
+   *
+   * Once domains are defined via ordinary domain constructors, this
+   * function will grab a set of domains and say which domain overlaps
+   * with which one.
+   */
+  void set_overlapper_info(set<shared_ptr<overlapping_domain> > 
+			   domains);
 }
 
 #endif // __DDM_HPP__

--- a/include/diff_op.hpp
+++ b/include/diff_op.hpp
@@ -2,6 +2,9 @@
  *
  *  \brief Differential operators for boundary-value problems (BVP)
  *  are declared here.
+ *  \file diff_op.cpp
+ *  \brief Implementations and instantiations of the functions and
+ *  classes declared in diff_op.hpp
  */
 
 #ifndef __DIFF_OP_HPP__

--- a/include/error.hpp
+++ b/include/error.hpp
@@ -2,10 +2,14 @@
  * \brief Throw exceptions instead of using GSL error handler function which
  * prefers to call abort().
  *
-
+ * 
  * Remember to put `gsl_set_error_handler(&errorHandler);' in the
  * main() loops when including this header file; otherwise it's
  * useless!
+ *
+ *  \file error.cpp
+ *  \brief Implementations and instantiations of the functions and
+ *  classes declared in error.hpp
  */
 
 #ifndef __ERROR_HPP__

--- a/include/func.hpp
+++ b/include/func.hpp
@@ -5,6 +5,10 @@
  * to \f$\mathbb{R}\f$ instead of using naked function pointers, meant
  * to be a class that can be derived from. Wrappers for GSL-style
  * function pointers are also provided.
+ *
+ *  \file func.cpp
+ *  \brief Implementations and instantiations of the functions and
+ *  classes declared in func.hpp
  */
 #ifndef __FUNC_HPP__
 #define __FUNC_HPP__
@@ -69,21 +73,39 @@
   };
 
 
-  //A function wrapper for calling GSL derivatives. 
+  ///A function wrapper for calling GSL derivatives. 
   class gsl_function_wrapper{
   public:
+    /*! \brief Wraps a realfunc for the GSL
+     *
+     * Turns a realfunc into a one-dimensional GSL function from
+     * \f$\mathbb{R}\f$ to \f$\mathbb{R}\f$ by localising the function
+     * at a point along a given direction.
+     *
+     * \param f - The realfunc to wrap.
+     * \param p - The point on which to localise.
+     * \param idx - The direction along with to localise.
+     */
     gsl_function_wrapper(const realfunc &f, point p, size_t idx);
+    /// Same as the above constructor
     void set_params( const realfunc &f, point p, size_t idx);
+    /// Pointer to the GSL-style function
     gsl_function* get_gsl_function() const;
 
+    /// The actual function that's passed to the GSL
     static double takemyaddress(double xi, void* nothing);
   private:
+    /// No nontrivial construction!
     gsl_function_wrapper();
     
+    /// The point at which the function is localised
     static point x;
+    /// The direction along which the function is localised
     static size_t index;    
+    /// The function being localised
     static realfunc myfunc;
-
+    
+    /// The GSL function formed with the above data.
     static gsl_function* f;
   };
 }

--- a/include/interpolator.hpp
+++ b/include/interpolator.hpp
@@ -4,6 +4,9 @@
  * \lambda_\xi \varphi_\xi(x) \f$, templated by RBFs, and its helper
  * classes are declared here.
  * 
+ *  \file interpolator.cpp
+ *  \brief Implementations and instantiations of the functions and
+ *  classes declared in interpolator.hpp
  */
 #ifndef __INTERPOLATOR_HPP__
 #define __INTERPOLATOR_HPP__
@@ -251,6 +254,9 @@
     interp_values operator-(const interp_values& w) const;
     interp_values operator*(const interp_values& w) const;
     interp_values operator/(const interp_values& w) const;
+
+    interp_values operator+(const double a) const;
+    interp_values operator-(const double a) const;
     interp_values operator*(const double a) const;
     interp_values operator/(const double a) const;
     //@}
@@ -271,17 +277,34 @@
 
     template<typename RBF> friend class interpolator;
     
-    /// For comfortable syntax
-    friend interp_values operator*(double a, const interp_values& w); 
+    /*! \name Arithmetic operators
+     * For comfortable syntax.
+     */
+    //@{
+    friend interp_values operator+(double a, const interp_values& v); 
+    friend interp_values operator-(double a, const interp_values& v); 
+    friend interp_values operator*(double a, const interp_values& v); 
+    friend interp_values operator/(double a, const interp_values& v); 
+    //@}
   };
 
-  ///For comfortable syntax
-  //@{
+
+
+  /// For comfortable syntax
   template <typename RBF>
   interpolator<RBF> operator*(double a, const interpolator<RBF>& u)
   {
     return u*a;
   }
+
+  /*! \name Arithmetic operators
+   * For comfortable syntax.
+   */
+  //@{
+  interp_values operator+(double a, const interp_values& v); 
+  interp_values operator-(double a, const interp_values& v); 
+  interp_values operator*(double a, const interp_values& v); 
+  interp_values operator/(double a, const interp_values& v); 
   //@}
 }
 

--- a/include/linalg.hpp
+++ b/include/linalg.hpp
@@ -4,6 +4,10 @@
  * This header file puts some C++ wrappers around the GSL vector and
  * matrix structures plus providing a C++ interface to some pertinent
  * BLAS and general linear algebra routines.
+ * 
+ *  \file linalg.cpp
+ *  \brief Implementations and instantiations of the functions and
+ *  classes declared in linalg.hpp
  */
 
 #ifndef __LINALG_HPP__
@@ -392,8 +396,8 @@
 
 namespace linalg{ //Non-member functions.
   
-  // I/O
-
+  /// I/O
+  //@{
   ///Stream insertion operator
   std::ostream& operator<<(std::ostream& os, const vector &v);
   ///Stream extraction operator
@@ -402,8 +406,10 @@
   std::ostream& operator<<(std::ostream& os, const matrix &M);
   ///Stream extraction operator
   matrix operator>>(std::istream& is, matrix& v);
+  //@}
 
-  //Some arithmetic functions for comfortable syntax.
+  ///Some arithmetic functions for comfortable syntax.
+  //@{
 
   /// Scale a vector.
   vector operator*(double a, const vector& v);
@@ -422,7 +428,7 @@
   double det(matrix& A);        
   /// L2 condition number, computed with SVD.
   double cond(matrix& A);
-
+  //@}
 }
 
 //Inlined functions

--- a/include/rbf.hpp
+++ b/include/rbf.hpp
@@ -1,5 +1,9 @@
 /*! \file  rbf.hpp
  * \brief Radial basis functions declared here.
+ *
+ *  \file rbf.cpp
+ *  \brief Implementations and pinstantiations of the functions and
+ *  classes declared in rbf.hpp
  */
 #ifndef __RBF_HPP__
 #define __RBF_HPP__

--- a/include/utils.hpp
+++ b/include/utils.hpp
@@ -1,5 +1,9 @@
 /*! \file utils.hpp
  * \brief Defines some miscellaneous functions.
+ *
+ *  \file utils.cpp
+ *  \brief Implementations and instantiations of the functions and
+ *  classes declared in utils.hpp
  */
 
 #ifndef __UTILS_HPP__ 

--- a/interpolator.cpp
+++ b/interpolator.cpp
@@ -675,6 +675,20 @@
   }
 
   interp_values
+  interp_values::operator+(const double a) const
+  {
+    vector w(v.size(),a);
+    return interp_values(v+w,rbfs_hash);
+  }
+  
+  interp_values
+  interp_values::operator-(const double a) const
+  {
+    vector w(v.size(),a);
+    return interp_values(v-w, rbfs_hash);
+  }
+
+  interp_values
   interp_values::operator*(const double a) const
   {
     return interp_values(v*a,rbfs_hash);
@@ -687,9 +701,27 @@
   }
 
   //*************** Non-member global friend functions ******************
+
+  interp_values operator+(double a, const interp_values&  v)
+  { 
+    return v+a;
+  } 
+
+  interp_values operator-(double a, const interp_values&  v)
+  { 
+    linalg::vector w(v.v.size(),a);
+    return interp_values(w-v.v, v.rbfs_hash);
+  } 
+
   interp_values operator*(double a, const interp_values&  v)
   { 
-    return interp_values(a*v.v, v.rbfs_hash);
+    return v*a;
+  } 
+
+  interp_values operator/(double a, const interp_values&  v)
+  { 
+    linalg::vector w(v.v.size(),a);
+    return interp_values(w/v.v, v.rbfs_hash);
   } 
 
   //Instantiations

--- a/linalg.cpp
+++ b/linalg.cpp
@@ -601,9 +601,8 @@
     }
 
     for(size_t i = 0; i < x -> size; i++)
-      //if(gsl_fcmp(gsl_vector_get(x,i), gsl_vector_get(w.x,i), eps) != 0)
-      if (gsl_vector_get(x,i) != gsl_vector_get(w.x,i))
-	return false;
+      if(gsl_fcmp(gsl_vector_get(x,i), gsl_vector_get(w.x,i), eps) != 0)
+      	return false;
     
     return true;
   }

--- a/utils.cpp
+++ b/utils.cpp
@@ -100,10 +100,10 @@
   void show_exception(error_handling::error exc){
     using namespace std;
 
-    cout << "Caught an exception!" << endl;
-    cout << exc.reason << endl;
-    cout << "On line " << exc.line << endl;
-    cout << "From file " << exc.file << endl;
+    cerr << "Caught an exception!" << endl;
+    cerr << exc.reason << endl;
+    cerr << "On line " << exc.line << endl;
+    cerr << "From file " << exc.file << endl;
   }
 
 }