matrix.h

Go to the documentation of this file.

/* 
 * Scythe Statistical Library Copyright (C) 2000-2002 Andrew D. Martin
 * and Kevin M. Quinn; 2002-present Andrew D. Martin, Kevin M. Quinn,
 * and Daniel Pemstein.  All Rights Reserved.
 *
 * This program is free software; you can redistribute it and/or
 * modify under the terms of the GNU General Public License as
 * published by Free Software Foundation; either version 2 of the
 * License, or (at your option) any later version.  See the text files
 * COPYING and LICENSE, distributed with this source code, for further
 * information.
 * --------------------------------------------------------------------
 *  scythe's/matrix.h
 *
 */

#ifndef SCYTHE_MATRIX_H
#define SCYTHE_MATRIX_H

#include <climits>
#include <iostream>
#include <iomanip>
#include <string>
#include <sstream>
#include <fstream>
#include <iterator>
#include <algorithm>
//#include <numeric>
#include <functional>
#include <list>

#ifdef SCYTHE_COMPILE_DIRECT
#include "defs.h"
#include "algorithm.h"
#include "error.h"
#include "datablock.h"
#include "matrix_random_access_iterator.h"
#include "matrix_forward_iterator.h"
#include "matrix_bidirectional_iterator.h"
#ifdef SCYTHE_LAPACK
#include "lapack.h"
#endif
#else
#include "scythestat/defs.h"
#include "scythestat/algorithm.h"
#include "scythestat/error.h"
#include "scythestat/datablock.h"
#include "scythestat/matrix_random_access_iterator.h"
#include "scythestat/matrix_forward_iterator.h"
#include "scythestat/matrix_bidirectional_iterator.h"
#ifdef SCYTHE_LAPACK
#include "scythestat/lapack.h"
#endif
#endif

namespace scythe {

  namespace { // make the uint typedef local to this file
    /* Convenience typedefs */
    typedef unsigned int uint;
  }

  /* Forward declare the matrix multiplication functions for *= to use
   * within Matrix proper .
   */

  template <typename T_type, matrix_order ORDER, matrix_style STYLE,
            matrix_order L_ORDER, matrix_style L_STYLE,
            matrix_order R_ORDER, matrix_style R_STYLE>
  Matrix<T_type, ORDER, STYLE>
  operator* (const Matrix<T_type, L_ORDER, L_STYLE>& lhs,
             const Matrix<T_type, R_ORDER, R_STYLE>& rhs);


  template <matrix_order L_ORDER, matrix_style L_STYLE,
            matrix_order R_ORDER, matrix_style R_STYLE, typename T_type>
  Matrix<T_type, L_ORDER, Concrete>
  operator* (const Matrix<T_type, L_ORDER, L_STYLE>& lhs,
             const Matrix<T_type, R_ORDER, R_STYLE>& rhs);

  /* forward declaration of the matrix class */
  template <typename T_type, matrix_order ORDER, matrix_style STYLE>
  class Matrix;

  template<typename T_elem, typename T_iter, 
           matrix_order O, matrix_style S>
  class ListInitializer {
    // An unbound friend
    template <class T, matrix_order OO, matrix_style SS>
    friend class Matrix;
    
    public:
      ListInitializer (T_elem val, T_iter begin, T_iter end, 
                       Matrix<T_elem,O,S>* matrix)
        : vals_ (),
          iter_ (begin),
          end_ (end),
          matrix_ (matrix),
          populated_ (false)
      {
        vals_.push_back(val);
      }

      ~ListInitializer ()
      {
        if (! populated_)
          populate();
      }

      ListInitializer &operator, (T_elem x)
      {
        vals_.push_back(x);
        return *this;
      }

    private:
      void populate ()
      {
        typename std::list<T_elem>::iterator vi = vals_.begin();

        while (iter_ < end_) {
          if (vi == vals_.end())
            vi = vals_.begin();
          *iter_ = *vi;
          ++iter_; ++vi;
        }

        populated_ = true;
      }

      std::list<T_elem> vals_;
      T_iter iter_;
      T_iter end_;
      Matrix<T_elem, O, S>* matrix_;
      bool populated_;
  };
  
  template <matrix_order ORDER = Col, matrix_style STYLE = Concrete>
  class Matrix_base
  {
    protected:
      /**** CONSTRUCTORS ****/

      /* Default constructor */
      Matrix_base ()
        : rows_ (0),
          cols_ (0),
          rowstride_ (0),
          colstride_ (0),
          storeorder_ (ORDER)
      {}

      /* Standard constructor */
      Matrix_base (uint rows, uint cols)
        : rows_ (rows),
          cols_ (cols),
          storeorder_ (ORDER)
      {
        if (ORDER == Col) {
          rowstride_ = 1;
          colstride_ = rows;
        } else {
          rowstride_ = cols;
          colstride_ = 1;
        }
      }

      /* Copy constructors 
       *
       * The first version handles matrices of the same order and
       * style.  The second handles matrices with different
       * orders/styles.  The same- templates are more specific,
       * so they will always catch same- cases.
       *
       */

      Matrix_base (const Matrix_base &m)
        : rows_ (m.rows()),
          cols_ (m.cols()),
          rowstride_ (m.rowstride()),
          colstride_ (m.colstride())
      {
        if (STYLE == View)
          storeorder_ = m.storeorder();
        else
          storeorder_ = ORDER;
      }

      template <matrix_order O, matrix_style S>
      Matrix_base (const Matrix_base<O, S> &m)
        : rows_ (m.rows()),
          cols_ (m.cols())
      {
        if (STYLE == View) {
          storeorder_ = m.storeorder();
          rowstride_ = m.rowstride();
          colstride_ = m.colstride();
         } else {
          storeorder_ = ORDER;
          if (ORDER == Col) {
            rowstride_ = 1;
            colstride_ = rows_;
          } else {
            rowstride_ = cols_;
            colstride_ = 1;
          }
         }
      }

      /* Submatrix constructor */
      template <matrix_order O, matrix_style S>
      Matrix_base (const Matrix_base<O, S> &m,
          uint x1, uint y1, uint x2, uint y2)
        : rows_ (x2 - x1 + 1),
          cols_ (y2 - y1 + 1),
          rowstride_ (m.rowstride()),
          colstride_ (m.colstride()),
          storeorder_ (m.storeorder())
      {
        /* Submatrices always have to be views, but the whole
         * concrete-view thing is just a policy maintained by the
         * software.  Therefore, we need to ensure this constructor
         * only returns views.  There's no neat way to do it but this
         * is still a compile time check, even if it only reports at
         * run-time.  Of course, we should never get this far.
         */
        if (STYLE == Concrete) {
          SCYTHE_THROW(scythe_style_error, 
              "Tried to construct a concrete submatrix (Matrix_base)!");
        }
      }


      /**** DESTRUCTOR ****/

      ~Matrix_base ()
      {}

      /**** OPERRATORS ****/

      // I'm defining one just to make sure we don't get any subtle
      // bugs from defaults being called.
      Matrix_base& operator=(const Matrix_base& m)
      {
        SCYTHE_THROW(scythe_unexpected_default_error,
          "Unexpected call to Matrix_base default assignment operator");
      }

      /**** MODIFIERS ****/

      /* Make this Matrix_base an exact copy of the matrix passed in.
       * Just like an assignment operator but can be called from a derived
       * class w/out making the = optor public and doing something
       * like
       * *(dynamic_cast<Matrix_base *> (this)) = M;
       * in the derived class.
       *
       * Works across styles, but should be used with care
       */
      template <matrix_order O, matrix_style S>
      inline void mimic (const Matrix_base<O, S> &m)
      {
        rows_ = m.rows();
        cols_ = m.cols();
        rowstride_ = m.rowstride();
        colstride_ = m.colstride();
        storeorder_ = m.storeorder();
      }

      /* Reset the dimensions of this Matrix_base.
       *
       * TODO This function is a bit of an interface weakness.  It
       * assumes a resize always means a fresh matrix (concrete or
       * view) with no slack between dims and strides. This happens to
       * always be the case when it is called, but it tightly couples
       * Matrix_base and extending classes.  Not a big issue (since
       * Matrix is probably the only class that will ever extend this)
       * but something to maybe fix down the road.
       */
      inline void resize (uint rows, uint cols)
      {
        rows_ = rows;
        cols_ = cols;

        if (ORDER == Col) {
          rowstride_ = 1;
          colstride_ = rows;
        } else {
          rowstride_ = cols;
          colstride_ = 1;
        }

        storeorder_ = ORDER;
      }
      
    public:
      /**** ACCESSORS ****/

      inline uint size () const
      {
        return (rows() * cols());
      }

      inline uint rows() const
      {
        return rows_;
      }

      inline uint cols () const
      {
        return cols_;
      }

      inline matrix_order order() const
      {
        return ORDER;
      }

      inline matrix_style style() const
      {
        return STYLE;
      }

      inline matrix_order storeorder () const
      {
        return storeorder_;
      }

      inline uint rowstride () const
      {
        return rowstride_;
      }
      
      inline uint colstride () const
      {
        return colstride_;
      }

      inline uint max_size () const
      {
        return UINT_MAX;
      }

      inline bool isScalar () const
      {
        return (rows() == 1 && cols() == 1);
      }

      inline bool isRowVector () const
      {
        return (rows() == 1);
      }
      
      inline bool isColVector () const
      {
        return (cols() == 1);
      }

      inline bool isVector () const
      {
        return (cols() == 1 || rows() == 1);
      }
      
      inline bool isSquare () const
      {
        return (cols() == rows());
      }

      inline bool isNull () const
      {
        return (rows() == 0);
      }

      inline bool empty () const
      {
        return (rows() == 0);
      }
      

      /**** HELPERS ****/

      inline bool inRange (uint i) const
      {
        return (i < size());
      }

      inline bool inRange (uint i, uint j) const
      {
        return (i < rows() && j < cols());
      }

    protected:
      /* These methods take offsets into a matrix and convert them
       * into that actual position in the referenced memory block,
       * taking stride into account.  Protection is debatable.  They
       * could be useful to outside functions in the library but most
       * callers should rely on the public () operator in the derived
       * class or iterators.
       *
       * Note that these are very fast for concrete matrices but not
       * so great for views.  Of course, the type checks are done at
       * compile time.
       */
      
      /* Turn an index into a true offset into the data. */
      inline uint index (uint i) const
      {
        if (STYLE == View) {
          if (ORDER == Col) {
            uint col = i / rows();
            uint row = i % rows();
            return (index(row, col));
          } else {
            uint row = i / cols();
            uint col = i % cols();
            return (index(row, col));
          }
        } else
          return(i);
      }

      /* Turn an i, j into an index. */
      inline uint index (uint row, uint col) const
      {
        if (STYLE == Concrete) {
          if (ORDER == Col)
            return (col * rows() + row);
          else
            return (row * cols() + col);
        } else { // view
          if (storeorder_ == Col)
            return (col * colstride() + row);
          else
            return (row * rowstride() + col);
        }
      }

    
    /**** INSTANCE VARIABLES ****/
    protected:
      uint rows_;   // # of rows
      uint cols_;   // # of cols

    private:
      /* The derived class shouldn't have to worry about this
       * implementation detail.
       */
      uint rowstride_;   // the in-memory number of elements from the
      uint colstride_;   // beginning of one column(row) to the next
      matrix_order storeorder_; // The in-memory storage order of this
                                // matrix.  This will always be the
                                // same as ORDER for concrete
                                // matrices but views can look at
                                // matrices with storage orders that
                                // differ from their own.
                                // TODO storeorder is always known at
                                // compile time, so we could probably
                                // add a third template param to deal
                                // with this.  That would speed up
                                // views a touch.  Bit messy maybe.
  };

  /**** MATRIX CLASS ****/

  template <typename T_type = double, matrix_order ORDER = Col, 
            matrix_style STYLE = Concrete>
  class Matrix : public Matrix_base<ORDER, STYLE>,
                 public DataBlockReference<T_type>
  {
    public:
      /**** TYPEDEFS ****/

      /* Iterator types */

      typedef matrix_random_access_iterator<T_type, ORDER, ORDER, STYLE>
        iterator;

      typedef const_matrix_random_access_iterator<T_type, ORDER, ORDER,
                                                  STYLE> const_iterator;

      typedef typename 
        std::reverse_iterator<matrix_random_access_iterator<T_type, 
                              ORDER, ORDER, STYLE> > reverse_iterator;

      typedef typename
        std::reverse_iterator<const_matrix_random_access_iterator
                              <T_type, ORDER, ORDER, STYLE> >
        const_reverse_iterator;

      typedef matrix_forward_iterator<T_type, ORDER, ORDER, STYLE>
        forward_iterator;

      typedef const_matrix_forward_iterator<T_type, ORDER, ORDER, STYLE>
        const_forward_iterator;

      typedef matrix_bidirectional_iterator<T_type, ORDER, ORDER, STYLE>
        bidirectional_iterator;

      typedef const_matrix_bidirectional_iterator<T_type, ORDER, ORDER,
                                  STYLE> const_bidirectional_iterator;

      typedef typename 
        std::reverse_iterator<matrix_bidirectional_iterator<T_type, 
                ORDER, ORDER, STYLE> > reverse_bidirectional_iterator;

      typedef typename
        std::reverse_iterator<const_matrix_bidirectional_iterator
                              <T_type, ORDER, ORDER, STYLE> >
        const_reverse_bidirectional_iterator;

      typedef T_type ttype;
    
    private:
      /* Some convenience typedefs */
      typedef DataBlockReference<T_type> DBRef;
      typedef Matrix_base<ORDER, STYLE> Base;
      
    public:
      /**** CONSTRUCTORS ****/

      Matrix ()
        : Base (),
          DBRef ()
      {
      }

      Matrix (T_type element)
        : Base (1, 1),
          DBRef (1)
      {
        data_[Base::index(0)] = element;  // ALWAYS use index()
      }

      Matrix (uint rows, uint cols, bool fill = true,
          T_type fill_value = 0)
        : Base (rows, cols),
          DBRef (rows * cols)
      {
        // TODO Might use iterator here for abstraction.
        if (fill)
          for (uint i = 0; i < Base::size(); ++i)
            data_[Base::index(i)] = fill_value; // we know data contig
      }

      template <typename T_iterator>
      Matrix (uint rows, uint cols, T_iterator it)
        : Base (rows, cols),
          DBRef (rows * cols)
      {
        // TODO again, should probably use iterator
        for (uint i = 0; i < Base::size(); ++i) {
          data_[Base::index(i)] = *it; // we know data_ is contig
          ++it;
        }
      }

      Matrix (const std::string& path, bool oldstyle=false)
        : Base (),
          DBRef ()
      {
        std::ifstream file(path.c_str());
        SCYTHE_CHECK_10(! file, scythe_file_error,
            "Could not open file " << path);

        if (oldstyle) {
          uint rows, cols;
          file >> rows >> cols;
          resize(rows, cols);
          std::copy(std::istream_iterator<T_type> (file), 
                    std::istream_iterator<T_type>(), begin_f<Row>());
        } else {
          std::string row;

          unsigned int cols = -1;
          std::vector<std::vector<T_type> > vals;
          unsigned int rows = 0;
          while (std::getline(file, row)) {
            std::vector<T_type> column;
            std::istringstream rowstream(row);
            std::copy(std::istream_iterator<T_type> (rowstream),
                 std::istream_iterator<T_type>(),
                 std::back_inserter(column));

            if (cols == -1)
              cols = (unsigned int) column.size();

            SCYTHE_CHECK_10(cols != column.size(), scythe_file_error,
                "Row " << (rows + 1) << " of input file has "
                << column.size() << " elements, but should have " << cols);

            vals.push_back(column);
            rows++;
          }

          resize(rows, cols);
          for (unsigned int i = 0; i < rows; ++i)
            operator()(i, _) = Matrix<T_type>(1, cols, vals[i].begin());
        }
      }

      /* Copy constructors. Uses template args to set up correct
       * behavior for both concrete and view matrices.  The branches
       * are no-ops and get optimized away at compile time.
       *
       * We have to define this twice because we must explicitly
       * override the default copy constructor; otherwise it is the
       * most specific template in a lot of cases and causes ugliness.
       */

      Matrix (const Matrix& M)
        : Base (M), // this call deals with concrete-view conversions
          DBRef ()
      {
        if (STYLE == Concrete) {
          this->referenceNew(M.size());
          scythe::copy<ORDER,ORDER>(M, *this);
        } else // STYLE == View
          this->referenceOther(M);
      }

      template <matrix_order O, matrix_style S>
      Matrix (const Matrix<T_type, O, S> &M)
        : Base (M), // this call deals with concrete-view conversions
          DBRef ()
      {
        if (STYLE == Concrete) {
          this->referenceNew(M.size());
          scythe::copy<ORDER,ORDER> (M, *this);
        } else // STYLE == View
          this->referenceOther(M);
      }

      template <typename S_type, matrix_order O, matrix_style S>
      Matrix (const Matrix<S_type, O, S> &M)
        : Base(M), // this call deal with concrete-view conversions
          DBRef (M.size())
      {
        scythe::copy<ORDER,ORDER> (M, *this);
      }

      template <matrix_order O, matrix_style S>
      Matrix (const Matrix<T_type, O, S> &M,
          uint x1, uint y1, uint x2, uint y2)
        : Base(M, x1, y1, x2, y2),
          DBRef(M, Base::index(x1, y1))
      {
        /* Submatrices always have to be views, but the whole
         * concrete-view thing is just a policy maintained by the
         * software.  Therefore, we need to ensure this constructor
         * only returns views.  There's no neat way to do it but this
         * is still a compile time check, even if it only reports at
         * run-time.
         */
        if (STYLE == Concrete) {
          SCYTHE_THROW(scythe_style_error, 
              "Tried to construct a concrete submatrix (Matrix)!");
        }
      }

    public:
      /**** DESTRUCTOR ****/

      ~Matrix() {}

      /**** COPY/REFERENCE METHODS ****/

      /* Make this matrix a view of another's data. If this matrix's
       * previous datablock is not viewed by any other object it is
       * deallocated.  Concrete matrices cannot be turned into views
       * at run-time!  Therefore, we generate an error here if *this
       * is concrete.
       */

      template <matrix_order O, matrix_style S>
      inline void reference (const Matrix<T_type, O, S> &M)
      {
        if (STYLE == Concrete) {
          SCYTHE_THROW(scythe_style_error, 
              "Concrete matrices cannot reference other matrices");
        } else {
          this->referenceOther(M);
          this->mimic(M);
        }
      }

      inline Matrix<T_type, ORDER, Concrete> copy () const
      {
        Matrix<T_type, ORDER> res (Base::rows(), Base::cols(), false);
        std::copy(begin_f(), end_f(), res.begin_f());

        return res;
      }

      /* Make this matrix a copy of another. The matrix retains its
       * own order and style in this case, because that can't change
       * at run time.
       */
      template <matrix_order O, matrix_style S>
      inline void copy (const Matrix<T_type, O, S>& M)
      {
        resize2Match(M);
        scythe::copy<ORDER,ORDER> (M, *this);
      }

      /**** INDEXING OPERATORS ****/
      
      inline T_type& operator[] (uint i)
      {
        SCYTHE_CHECK_30 (! Base::inRange(i), scythe_bounds_error,
            "Index " << i << " out of range");

        return data_[Base::index(i)];
      }
      
      inline T_type& operator[] (uint i) const
      {
        SCYTHE_CHECK_30 (! Base::inRange(i), scythe_bounds_error,
            "Index " << i << " out of range");

        return data_[Base::index(i)];
      }

      inline T_type& operator() (uint i)
      {
        SCYTHE_CHECK_30 (! Base::inRange(i), scythe_bounds_error,
            "Index " << i << " out of range");

        return data_[Base::index(i)];
      }
      
      inline T_type& operator() (uint i) const
      {
        SCYTHE_CHECK_30 (! Base::inRange(i), scythe_bounds_error,
          "Index " << i << " out of range");

        return data_[Base::index(i)];
      }

      inline T_type& operator() (uint i, uint j)
      {
        SCYTHE_CHECK_30 (! Base::inRange(i, j), scythe_bounds_error,
            "Index (" << i << ", " << j << ") out of range");

        return data_[Base::index(i, j)];
      }
        
      inline T_type& operator() (uint i, uint j) const
      {
        SCYTHE_CHECK_30 (! Base::inRange(i, j), scythe_bounds_error,
            "Index (" << i << ", " << j << ") out of range");

        return data_[Base::index(i, j)];
      }

      /**** SUBMATRIX OPERATORS ****/


      /* Submatrices are always views.  An extra (but relatively
       * cheap) copy constructor call is made when mixing and matching
       * orders like
       *
       * Matrix<> A;
       * ...
       * Matrix<double, Row> B = A(2, 2, 4, 4);
       *
       * It is technically possible to get around this, by providing
       * templates of each function of the form
       * template <matrix_order O>
       * Matrix<T_type, O, View> operator() (...)
       *
       * but the syntax to call them (crappy return type inference):
       *
       * Matrix<double, Row> B = A.template operator()<Row>(2, 2, 4, 4)
       *
       * is such complete gibberish that I don't think this is worth
       * the slight optimization.
       */
      
      inline Matrix<T_type, ORDER, View> 
      operator() (uint x1, uint y1, uint x2, uint y2)
      {
        SCYTHE_CHECK_20 (! Base::inRange(x1, y1) 
            || ! Base::inRange(x2, y2)
            || x1 > x2 || y1 > y2,
            scythe_bounds_error,
            "Submatrix (" << x1 << ", " << y1 << ") ; ("
            << x2 << ", " << y2 << ") out of range or ill-formed");
        
        return (Matrix<T_type, ORDER, View>(*this, x1, y1, x2, y2));
      }
      
      inline Matrix<T_type, ORDER, View> 
      operator() (uint x1, uint y1, uint x2, uint y2) const
      {
        SCYTHE_CHECK_20 (! Base::inRange(x1, y1) 
            || ! Base::inRange(x2, y2)
            || x1 > x2 || y1 > y2,
            scythe_bounds_error,
            "Submatrix (" << x1 << ", " << y1 << ") ; ("
            << x2 << ", " << y2 << ") out of range or ill-formed");

        return (Matrix<T_type, ORDER, View>(*this, x1, y1, x2, y2));
      }

      inline Matrix<T_type, ORDER, View> 
      operator() (const all_elements a, uint j)
      {
        SCYTHE_CHECK_20 (j >= Base::cols(), scythe_bounds_error,
            "Column vector index " << j << " out of range");

        return (Matrix<T_type, ORDER, View>
           (*this, 0, j, Base::rows() - 1, j));
      }
      
      inline Matrix<T_type, ORDER, View> 
      operator() (const all_elements a, uint j) const
      {
        SCYTHE_CHECK_20 (j >= Base::cols(), scythe_bounds_error,
            "Column vector index " << j << " out of range");

        return (Matrix<T_type, ORDER, View>
           (*this, 0, j, Base::rows() - 1, j));
      }

      inline Matrix<T_type, ORDER, View> 
      operator() (uint i, const all_elements b)
      {
        SCYTHE_CHECK_20 (i >= Base::rows(), scythe_bounds_error,
            "Row vector index " << i << " out of range");

        return (Matrix<T_type, ORDER, View>
            (*this, i, 0, i, Base::cols() - 1));
      }
      
      inline Matrix<T_type, ORDER, View> 
      operator() (uint i, const all_elements b) const
      {
        SCYTHE_CHECK_20 (i >= Base::rows(), scythe_bounds_error,
            "Row vector index " << i << " out of range");
        return (Matrix<T_type, ORDER, View>
            (*this, i, 0, i, Base::cols() - 1));
      } 

      /**** ASSIGNMENT OPERATORS ****/

       /*
       * As with the copy constructor, we need to
       * explicitly define the same-order-same-style assignment
       * operator or the default operator will take over.
       *
       * TODO With views, it may be desirable to auto-grow (and,
       * technically, detach) views to the null matrix.  This means
       * you can write something like:
       *
       * Matrix<double, Col, View> X;
       * X = ...
       *
       * and not run into trouble because you didn't presize.  Still,
       * not sure this won't encourage silly mistakes...need to think
       * about it.
       */

      Matrix& operator= (const Matrix& M)
      {
        if (STYLE == Concrete) {
          resize2Match(M);
          scythe::copy<ORDER,ORDER> (M, *this);
        } else {
#ifndef SCYTHE_VIEW_ASSIGNMENT_RECYCLE
          SCYTHE_CHECK_10 (Base::size() != M.size(),
              scythe_conformation_error,
              "LHS has dimensions (" << Base::rows() 
              << ", " << Base::cols()
              << ") while RHS has dimensions (" << M.rows() << ", "
              << M.cols() << ")");

          scythe::copy<ORDER,ORDER> (M, *this);
#else
          copy_recycle<ORDER,ORDER>(M, *this);
#endif
        }

        return *this;
      }
      
      template <matrix_order O, matrix_style S>
      Matrix &operator= (const Matrix<T_type, O, S> &M)
      {
        if (STYLE == Concrete) {
          resize2Match(M);
          scythe::copy<ORDER,ORDER> (M, *this);
        } else {
#ifndef SCYTHE_VIEW_ASSIGNMENT_RECYCLE
          SCYTHE_CHECK_10 (Base::size() != M.size(),
              scythe_conformation_error,
              "LHS has dimensions (" << Base::rows() 
              << ", " << Base::cols()
              << ") while RHS has dimensions (" << M.rows() << ", "
              << M.cols() << ")");

          scythe::copy<ORDER,ORDER> (M, *this);
#else
          copy_recycle<ORDER,ORDER>(M, *this);
#endif
        }

        return *this;
      }
      
      /* List-wise initialization behavior is a touch complicated.
       * List needs to be less than or equal to matrix in size and it
       * is copied into the matrix with R-style recycling.
       *
       * The one issue is that, if you want true assignment of a
       * scalar to a concrete matrix (resize the matrix to a scalar)
       * you need to be explicit:
       *
       * Matrix<> A(2, 2);
       * double x = 3;
       * ...
       * A = Matrix<>(x); // -> 3
       *
       * A = x; // -> 3 3
       *        //    3 3
       */

      ListInitializer<T_type, iterator, ORDER, STYLE> 
      operator=(T_type x)
      {
        return (ListInitializer<T_type, iterator, ORDER, STYLE> 
          (x, begin(),end(), this));
      }

      template <typename ITERATOR, matrix_order O, matrix_style S>
      Matrix &operator=(ListInitializer<T_type, ITERATOR, O, S> li)
      {
        li.populate();
        *this = *(li.matrix_);
        return *this;
      }

      /**** ARITHMETIC OPERATORS ****/

    private:
      /* Reusable chunk of code for element-wise operator
       * assignments.  Updates are done in-place except for the 1x1 by
       * nXm case, which forces a resize.
       */
      template <typename OP, matrix_order O, matrix_style S>
      inline Matrix&
      elementWiseOperatorAssignment (const Matrix<T_type, O, S>& M, 
                                     OP op)
      {
        SCYTHE_CHECK_10 (Base::size() != 1 && M.size() != 1 && 
            (Base::rows () != M.rows() || Base::cols() != M.cols()),
            scythe_conformation_error,
            "Matrices with dimensions (" << Base::rows() 
            << ", " << Base::cols()
            << ") and (" << M.rows() << ", " << M.cols()
            << ") are not conformable");
        
        if (Base::size() == 1) { // 1x1 += nXm
          T_type tmp = (*this)(0);
          resize2Match(M);
          std::transform(M.begin_f<ORDER>(), M.end_f<ORDER>(), 
              begin_f(), std::bind1st(op, tmp));
        } else if (M.size() == 1) { // nXm += 1x1
          std::transform(begin_f(), end_f(), begin_f(),
              std::bind2nd(op, M(0)));
        } else { // nXm += nXm
            std::transform(begin_f(), end_f(), M.begin_f<ORDER>(), 
                begin_f(), op);
        }

        return *this;
      }

    public:
      template <matrix_order O, matrix_style S>
      inline Matrix& operator+= (const Matrix<T_type, O, S> &M)
      {
        return elementWiseOperatorAssignment(M, std::plus<T_type> ());
      }

      inline Matrix& operator+= (T_type x)
      {
        return elementWiseOperatorAssignment(Matrix(x), 
            std::plus<T_type> ());
      }
      
      template <matrix_order O, matrix_style S>
      inline Matrix& operator-= (const Matrix<T_type, O, S> &M)
      {
        return elementWiseOperatorAssignment(M, std::minus<T_type> ());
      }

      inline Matrix& operator-= (T_type x)
      {
        return elementWiseOperatorAssignment(Matrix(x), 
            std::minus<T_type> ());
      }
      
      template <matrix_order O, matrix_style S>
      inline Matrix& operator%= (const Matrix<T_type, O, S> &M)
      {
        return elementWiseOperatorAssignment(M, 
            std::multiplies<T_type> ());
      }

      inline Matrix& operator%= (T_type x)
      {
        return elementWiseOperatorAssignment(Matrix(x), 
            std::multiplies<T_type> ());
      }
      
      template <matrix_order O, matrix_style S>
      inline Matrix& operator/= (const Matrix<T_type, O, S> &M)
      {
        return elementWiseOperatorAssignment(M, std::divides<T_type>());
      }

      inline Matrix& operator/= (T_type x)
      {
        return elementWiseOperatorAssignment(Matrix(x), 
            std::divides<T_type> ());
      }

      template <matrix_order O, matrix_style S>
      inline Matrix& operator^= (const Matrix<T_type, O, S> &M)
      {
        return elementWiseOperatorAssignment(M, 
            exponentiate<T_type>());
      }

      inline Matrix& operator^= (T_type x)
      {
        return elementWiseOperatorAssignment(Matrix(x), 
            exponentiate<T_type> ());
      }

      /* Matrix mult always disengages views because it generally
       * requires a resize.  We force a disengage in the one place it
       * isn't absolutely necessary(this->size()==1), for consistency.
       */

      template <matrix_order O, matrix_style S>
      Matrix& operator*= (const Matrix<T_type, O, S>& M)
      {
        /* Farm out the work to the plain old * operator and make this
         * matrix a reference (the only reference) to the result.  We
         * always have to create a new matrix here, so there is no
         * speed-up from using *=.
         */
        
        /* This saves a copy over 
         * *this = (*this) * M;
         * if we're concrete
         */
        Matrix<T_type, ORDER> res = (*this) * M;
        this->referenceOther(res);
        this->mimic(res);

        return *this;
      }

      inline Matrix& operator*= (T_type x)
      {
        return elementWiseOperatorAssignment(Matrix(x), 
            std::multiplies<T_type> ());
      }

      template <matrix_order O, matrix_style S> Matrix& kronecker
        (const Matrix<T_type, O, S>& M) { uint totalrows =
          Base::rows() * M.rows(); uint totalcols = Base::cols() *
            M.cols();
        // Even if we're a view, make this guy concrete.
        Matrix<T_type,ORDER> res(totalrows, totalcols, false);

        /* TODO: This the most natural way to write this in scythe
         * (with a small optimization based on ordering) but probably
         * not the fastest because it uses submatrix assignments.
         * Optimizations should be considered.
         */
        forward_iterator it = begin_f();
        if (ORDER == Row) {
          for (uint row = 0; row < totalrows; row += M.rows()) {
            for (uint col = 0; col < totalcols; col += M.cols()){
              res(row, col, row + M.rows() - 1, col + M.cols() - 1)
                 = (*it) * M;
              it++;
            }
          }
        } else {
          for (uint col = 0; col < totalcols; col += M.cols()) {
            for (uint row = 0; row < totalrows; row += M.rows()){
              res(row, col, row + M.rows() - 1, col + M.cols() - 1)
                = (*it) * M;
              it++;
            }
          }
        }
       
        this->referenceOther(res);
        this->mimic(res);

        return *this;
      }
        
      inline Matrix& kronecker (T_type x)
      {
        return elementWiseOperatorAssignment(Matrix(x), 
            std::multiplies<T_type> ());
      }

      /* Logical assignment operators */

      template <matrix_order O, matrix_style S>
      inline Matrix& operator&= (const Matrix<T_type, O, S> &M)
      {
        return elementWiseOperatorAssignment(M, 
            std::logical_and<T_type>());
      }

      inline Matrix& operator&= (T_type x)
      {
        return elementWiseOperatorAssignment(Matrix(x), 
            std::logical_and<T_type> ());
      }

      template <matrix_order O, matrix_style S>
      inline Matrix& operator|= (const Matrix<T_type, O, S> &M)
      {
        return elementWiseOperatorAssignment(M, 
            std::logical_or<T_type>());
      }

      inline Matrix& operator|= (T_type x)
      {
        return elementWiseOperatorAssignment(Matrix(x), 
            std::logical_or<T_type> ());
      }

      /**** MODIFIERS ****/

      /* Resize a matrix view.  resize() takes dimensions as
       * parameters while resize2Match() takes a matrix reference and
       * uses its dimensions.
       */

      void resize (uint rows, uint cols, bool preserve=false)
      {
        if (preserve) {
          /* TODO Optimize this case.  It is rather clunky. */
          Matrix<T_type, ORDER, View> tmp(*this);
          this->referenceNew(rows * cols);
          Base::resize(rows, cols);
          uint min_cols = std::min(Base::cols(), tmp.cols());
          uint min_rows = std::min(Base::rows(), tmp.rows());

          // TODO use iterators here perhaps
          if (ORDER == Col) {
            for (uint j = 0; j < min_cols; ++j)
              for (uint i = 0; i < min_rows; ++i)
                (*this)(i, j) = tmp(i, j);
          } else {
            for (uint i = 0; i < min_rows; ++i)
              for (uint j = 0; j < min_cols; ++j)
                (*this)(i, j) = tmp(i, j);
          }
        } else {
          this->referenceNew(rows * cols);
          Base::resize(rows, cols);
        }
      }

      template <typename T, matrix_order O, matrix_style S>
      inline void resize2Match(const Matrix<T, O, S> &M,
                               bool preserve=false)
      {
        resize(M.rows(), M.cols(), preserve);
      }

      /* Copy this matrix to a new datablock in contiguous storage */
      inline void detach ()
      {
        resize2Match(*this, true);
      }

      /* Swap operator: sort of a dual copy constructor.  Part of the
       * standard STL container interface. We only support swaps
       * between matrices of like order and style because things get
       * hairy otherwise.  The behavior of this for concrete matrices
       * is a little hairy in any case.
       *
       * Matrix<> A, B;
       * ... // fill in A and B
       * Matrix<double, Col, View> v1 = A(_, 1);
       * A.swap(B);
       * Matrix<double, Col, View> v2 = B(_, 1);
       * 
       * v1 == v2; // evaluates to true
       *
       */

      inline void swap (Matrix &M)
      {
        Matrix tmp = *this;

        /* This are just reference() calls, but we do this explicitly
         * here to avoid throwing errors on the concrete case.  While
         * having a concrete matrix reference another matrix is
         * generally a bad idea, it is safe when the referenced matrix
         * is concrete, has the same order, and gets deallocated (or
         * redirected at another block) like here.
         */

        this->referenceOther(M);
        this->mimic(M);

        M.referenceOther(tmp);
        M.mimic(tmp);
      }

      /**** ACCESSORS ****/

      /* Accessors that don't access the data itself (that don't rely
       * on T_type) are in Matrix_base
       */

      /* Are all the elements of this Matrix == 0 */

      inline bool isZero () const
      {
        const_forward_iterator last = end_f();
        return (last == std::find_if(begin_f(), last, 
          std::bind1st(std::not_equal_to<T_type> (), 0)));
      }

      /* M(i,j) == 0 when i != j */
      inline bool isDiagonal() const
      {
        if (! Base::isSquare())
          return false;
        /* Always travel in order.  It would be nice to use iterators
         * here, but we'd need to take views and their iterators are
         * too slow at the moment.
         * TODO redo with views and iterators if optimized.
         */
        if (ORDER == Row) {
          for (uint i = 0; i < Base::rows(); ++i) {
            for (uint j = 0; j < Base::cols(); ++j) {
              if (i != j && (*this)(i, j) != 0)
                return false;
            }
          }
        } else { // ORDER == Col
          for (uint j = 0; j < Base::cols(); ++j) {
            for (uint i = 0; i < Base::rows(); ++i) {
              if (i != j && (*this)(i, j) != 0)
                return false;
            }
          }
        }
        return true;
      }

      /* M(I, j) == 0 when i!= j and 1 when i == j */
      inline bool isIdentity () const
      {
        if (! Base::isSquare())
          return false;
        // TODO redo with views and iterator if optimized
        if (ORDER == Row) {
          for (uint i = 0; i < Base::rows(); ++i) {
            for (uint j = 0; j < Base::cols(); ++j) {
              if (i != j) {
                if ((*this)(i,j) != 0)
                  return false;
              } else if ((*this)(i,j) != 1)
                return false;
            }
          }
        } else { // ORDER == Col
          for (uint j = 0; j < Base::rows(); ++j) {
            for (uint i = 0; i < Base::cols(); ++i) {
              if (i != j) {
                if ((*this)(i,j) != 0)
                  return false;
              } else if ((*this)(i,j) != 1)
                return false;
            }
          }
        }
        return true;
      }

      /* M(i,j) == 0 when i < j */
      inline bool isLowerTriangular () const
      {
        if (! Base::isSquare())
          return false;
        // TODO view+iterator if optimized
        if (ORDER == Row) {
          for (uint i = 0; i < Base::rows(); ++i)
            for (uint j = i + 1; j < Base::cols(); ++j)
              if ((*this)(i,j) != 0)
                return false;
        } else {
          for (uint j = 0; j < Base::cols(); ++j)
            for (uint i = 0; i < j; ++i)
              if ((*this)(i,j) != 0)
                return false;
       }
        return true;
      }

      /* M(i,j) == 0 when i > j */
      inline bool isUpperTriangular () const
      {
        if (! Base::isSquare())
          return false;
        // TODO view+iterator if optimized
        if (ORDER == Row) {
          for (uint i = 0; i < Base::rows(); ++i)
            for (uint j = 0; j < i; ++j)
              if ((*this)(i,j) != 0)
                return false;
        } else {
          for (uint j = 0; j < Base::cols(); ++j)
            for (uint i = j + 1; i < Base::rows(); ++i)
              if ((*this)(i,j) != 0)
                return false;
       }
        return true;
      }

      inline bool isSingular() const
      {
        if (! Base::isSquare() || Base::isNull())
          return false;
        if ((~(*this)) == (T_type) 0)
          return true;
        return false;
      }

      /* Square and t(M) = M(inv(M) * t(M) == I */
      inline bool isSymmetric () const
      {
        if (! Base::isSquare())
          return false;
        // No point in order optimizing here
        for (uint i = 1; i < Base::rows(); ++i)
          for (uint j = 0; j < i; ++j)
            if ((*this)(i, j) != (*this)(j, i))
              return false;

        return true;
      }

      /* The matrix is square and t(A) = -A */
      inline bool isSkewSymmetric () const
      {
        if (! Base::isSquare())
          return false;
        // No point in order optimizing here
        for (uint i = 1; i < Base::rows(); ++i)
          for (uint j = 0; j < i; ++j)
            if ((*this)(i, j) != 0 - (*this)(j, i))
              return false;

        return true;
      }

      template <matrix_order O, matrix_style S>
      inline bool
      equals(const Matrix<T_type, O, S>& M) const
      {
        if (data_ == M.getArray() && STYLE == Concrete && S == Concrete)
          return true;
        else if (data_ == M.getArray() && Base::rows() == M.rows() 
                 && Base::cols() == M.cols()) {
          return true;
        } else if (this->isNull() && M.isNull())
          return true;
        else if (Base::rows() != M.rows() || Base::cols() != M.cols())
          return false;

        return std::equal(begin_f(), end_f(),
            M.template begin_f<ORDER>());
      }

      inline bool
      equals(T_type x) const
      {
        const_forward_iterator last = end_f();
        return (last == std::find_if(begin_f(), last, 
          std::bind1st(std::not_equal_to<T_type> (), x)));
      }


      /**** OTHER UTILITIES ****/

      inline T_type* getArray () const
      {
        return data_;
      }

      inline void
      save (const std::string& path, const char flag = 'n',
            const bool header = false) const
      {
        std::ofstream out;
        if (flag == 'n') {
          std::fstream temp(path.c_str(), std::ios::in);
          if (! temp)
            out.open(path.c_str(), std::ios::out);
          else {
            temp.close();
            SCYTHE_THROW(scythe_file_error, "Cannot overwrite file "
                << path << " when flag = n");
          }
        } else if (flag == 'o')
          out.open(path.c_str(), std::ios::out | std::ios::trunc);
        else if (flag == 'a')
          out.open(path.c_str(), std::ios::out | std::ios::app);
        else
          SCYTHE_THROW(scythe_invalid_arg, "Invalid flag: " << flag);

        if (! out)
          SCYTHE_THROW(scythe_file_error, 
              "Could not open file " << path); 

        if (header) {
          out << Base::rows() << " " << Base::cols();
          for (uint i = 0; i < Base::size(); ++i)
            out << " " << (*this)[i];
          out << std::endl;
        } else {
          for (uint i = 0; i < Base::rows(); ++i) {
            for (uint j = 0; j < Base::cols(); ++j)
              out << (*this)(i,j) << " ";
            out << "\n";
          }
        }
        out.close();
      }


      /**** ITERATOR FACTORIES ****/

      /* TODO Write some cpp macro code to reduce this to something
       * manageable.
       */

      /* Random Access Iterator Factories */
      
      /* Generalized versions */

      template <matrix_order I_ORDER>
      inline 
      matrix_random_access_iterator<T_type, I_ORDER, ORDER, STYLE>
      begin ()
      {
        return matrix_random_access_iterator<T_type, I_ORDER, ORDER,
                                             STYLE>(*this);
      }
      
      template <matrix_order I_ORDER>
      inline
      const_matrix_random_access_iterator<T_type, I_ORDER, ORDER, STYLE>
      begin() const
      {
        return const_matrix_random_access_iterator<T_type, I_ORDER,
                                                   ORDER, STYLE>
          (*this);
      }

      template <matrix_order I_ORDER>
      inline 
      matrix_random_access_iterator<T_type, I_ORDER, ORDER, STYLE>
      end ()
      {
        return (begin<I_ORDER>() + Base::size());
      }
      
      template <matrix_order I_ORDER>
      inline 
      const_matrix_random_access_iterator<T_type, I_ORDER, ORDER, STYLE>
      end () const
      {
        return (begin<I_ORDER>() + Base::size());
      }

      template <matrix_order I_ORDER>
      inline std::reverse_iterator<matrix_random_access_iterator<T_type,
                                   I_ORDER, ORDER, STYLE> >
      rbegin()
      {
        return std::reverse_iterator<matrix_random_access_iterator
                                     <T_type, I_ORDER, ORDER, STYLE> > 
               (end<I_ORDER>());
      }

      template <matrix_order I_ORDER>
      inline 
      std::reverse_iterator<const_matrix_random_access_iterator
                            <T_type, I_ORDER, ORDER, STYLE> > 
      rbegin() const
      {
        return std::reverse_iterator<const_matrix_random_access_iterator
                                     <T_type, I_ORDER, ORDER, STYLE> > 
        (end<I_ORDER>());
      }

      template <matrix_order I_ORDER>
      inline std::reverse_iterator<matrix_random_access_iterator
                                   <T_type, I_ORDER, ORDER, STYLE> >
      rend()
      {
        return std::reverse_iterator<matrix_random_access_iterator
                                     <T_type, I_ORDER, ORDER, STYLE> > 
               (begin<I_ORDER>());
      }

      template <matrix_order I_ORDER>
      inline 
      std::reverse_iterator<const_matrix_random_access_iterator
                            <T_type, I_ORDER, ORDER, STYLE> > 
      rend() const
      {
        return std::reverse_iterator<const_matrix_random_access_iterator
                                     <T_type, I_ORDER, ORDER, STYLE> > 
          (begin<I_ORDER>());
      }

      /* Specific versions --- the generalized versions force you
       * choose the ordering explicitly.  These definitions set up
       * in-order iteration as a default */
      
      inline iterator begin ()
      {
        return iterator(*this);
      }
      
      inline const_iterator begin() const
      {
        return const_iterator (*this);
      }

      inline iterator end ()
      {
        return (begin() + Base::size());
      }
      
      inline 
      const_iterator end () const
      {
        return (begin() + Base::size());
      }

      inline reverse_iterator rbegin()
      {
        return reverse_iterator (end());
      }

      inline const_reverse_iterator rbegin() const
      {
        return const_reverse_iterator (end());
      }

      inline reverse_iterator rend()
      {
        return reverse_iterator (begin());
      }

      inline const_reverse_iterator rend() const
      {
        return const_reverse_iterator (begin());
      }

      /* Forward Iterator Factories */

      /* Generalized versions */

      template <matrix_order I_ORDER>
      inline 
      matrix_forward_iterator<T_type, I_ORDER, ORDER, STYLE>
      begin_f ()
      {
        return matrix_forward_iterator<T_type, I_ORDER, ORDER,
                                             STYLE>(*this);
      }
      
      template <matrix_order I_ORDER>
      inline
      const_matrix_forward_iterator <T_type, I_ORDER, ORDER, STYLE>
      begin_f () const
      {
        return const_matrix_forward_iterator <T_type, I_ORDER,
                                                   ORDER, STYLE>
          (*this);
      }

      template <matrix_order I_ORDER>
      inline 
      matrix_forward_iterator<T_type, I_ORDER, ORDER, STYLE>
      end_f ()
      {
        return (begin_f<I_ORDER>().set_end());
      }
      
      template <matrix_order I_ORDER>
      inline 
      const_matrix_forward_iterator<T_type, I_ORDER, ORDER, STYLE>
      end_f () const
      {
        return (begin_f<I_ORDER>().set_end());
      }

      /* Default Versions */
      
      inline forward_iterator begin_f ()
      {
        return forward_iterator(*this);
      }
      
      inline const_forward_iterator begin_f () const
      {
        return const_forward_iterator (*this);
      }

      inline forward_iterator end_f ()
      {
        return (begin_f().set_end());
      }
      
      inline 
      const_forward_iterator end_f () const
      {
        return (begin_f().set_end());
      }

      /* Bidirectional Iterator Factories */

      /* Generalized versions */

      template <matrix_order I_ORDER>
      inline 
      matrix_bidirectional_iterator<T_type, I_ORDER, ORDER, STYLE>
      begin_bd ()
      {
        return matrix_bidirectional_iterator<T_type, I_ORDER, ORDER,
                                             STYLE>(*this);
      }
      
      template <matrix_order I_ORDER>
      inline
      const_matrix_bidirectional_iterator<T_type, I_ORDER, ORDER, STYLE>
      begin_bd () const
      {
        return const_matrix_bidirectional_iterator<T_type, I_ORDER,
                                                   ORDER, STYLE>
          (*this);
      }

      template <matrix_order I_ORDER>
      inline 
      matrix_bidirectional_iterator<T_type, I_ORDER, ORDER, STYLE>
      end_bd ()
      {
        return (begin_bd<I_ORDER>().set_end());
      }
      
      template <matrix_order I_ORDER>
      inline 
      const_matrix_bidirectional_iterator<T_type, I_ORDER, ORDER, STYLE>
      end_bd () const
      {
        return (begin_bd<I_ORDER>.set_end());
      }

      template <matrix_order I_ORDER>
      inline std::reverse_iterator<matrix_bidirectional_iterator<T_type,
                                   I_ORDER, ORDER, STYLE> >
      rbegin_bd ()
      {
        return std::reverse_iterator<matrix_bidirectional_iterator
                                     <T_type, I_ORDER, ORDER, STYLE> > 
               (end_bd<I_ORDER>());
      }

      template <matrix_order I_ORDER>
      inline 
      std::reverse_iterator<const_matrix_bidirectional_iterator
                            <T_type, I_ORDER, ORDER, STYLE> > 
      rbegin_bd () const
      {
        return std::reverse_iterator<const_matrix_bidirectional_iterator
                                     <T_type, I_ORDER, ORDER, STYLE> > 
        (end_bd<I_ORDER>());
      }

      template <matrix_order I_ORDER>
      inline std::reverse_iterator<matrix_bidirectional_iterator
                                   <T_type, I_ORDER, ORDER, STYLE> >
      rend_bd ()
      {
        return std::reverse_iterator<matrix_bidirectional_iterator
                                     <T_type, I_ORDER, ORDER, STYLE> > 
               (begin_bd<I_ORDER>());
      }

      template <matrix_order I_ORDER>
      inline 
      std::reverse_iterator<const_matrix_bidirectional_iterator
                            <T_type, I_ORDER, ORDER, STYLE> > 
      rend_bd () const
      {
        return std::reverse_iterator<const_matrix_bidirectional_iterator
                                     <T_type, I_ORDER, ORDER, STYLE> > 
          (begin_bd<I_ORDER>());
      }

      /* Specific versions --- the generalized versions force you
       * choose the ordering explicitly.  These definitions set up
       * in-order iteration as a default */
      
      inline bidirectional_iterator begin_bd ()
      {
        return bidirectional_iterator(*this);
      }
      
      inline const_bidirectional_iterator begin_bd() const
      {
        return const_bidirectional_iterator (*this);
      }

      inline bidirectional_iterator end_bd ()
      {
        return (begin_bd().set_end());
      }
      
      inline 
      const_bidirectional_iterator end_bd () const
      {
        return (begin_bd().set_end());
      }

      inline reverse_bidirectional_iterator rbegin_bd()
      {
        return reverse_bidirectional_iterator (end_bd());
      }

      inline const_reverse_bidirectional_iterator rbegin_bd () const
      {
        return const_reverse_bidirectional_iterator (end_bd());
      }

      inline reverse_bidirectional_iterator rend_bd ()
      {
        return reverse_bidirectional_iterator (begin_bd());
      }

      inline const_reverse_iterator rend_bd () const
      {
        return const_reverse_bidirectiona_iterator (begin_bd());
      }

    protected:
      /**** INSTANCE VARIABLES ****/

      /* I know the point of C++ is to force you to write 20 times
       * more code than should be necessary but "using" inherited ivs
       * is just stupid.
       */
      using DBRef::data_;  // refer to inherited data pointer directly 
      using Base::rows_;   // " # of rows directly
      using Base::cols_;   // " # of cols directly

  }; // end class Matrix

  /**** EXTERNAL OPERATORS ****/

  /* External operators include a range of binary matrix operations
   * such as tests for equality, and arithmetic.  Style
   * (concrete/view) of the returned matrix is that of the left hand
   * side parameter by default
   *
   * There is also a question of the ordering of the returned matrix.
   * We adopt the convention of returning a matrix ordered like that
   * of the left hand side argument, by default.
   *
   * Whenever there is only one matrix argument (lhs is scalar) we use
   * its order and style as the default.
   *
   * A general template version of each operator also exists and users
   * can coerce the return type to whatever they prefer using some
   * ugly syntax; ex:
   *
   * Matrix<> A; ...  Matrix<double, Row> B = operator*<Row,Concrete>
   *                                          (A, A);
   *
   * In general, the matrix class copy constructor will quietly
   * convert whatever matrix template is returned to the type of the
   * matrix it is being copied into on return, but one might want to
   * specify the type for objects that only exist for a second (ex:
   * (operator*<Row,Concrete>(A, A)).begin()).  Also, note that the
   * fact that we return concrete matrices by default does not
   * preclude the user from taking advantage of fast view copies.  It
   * is the template type of the object being copy-constructed that
   * matters---in terms of underlying implementation all matrices are
   * views, concrete matrices just maintain a particular policy.
   *
   * TODO Consider the best type for scalar args to these functions.
   * For the most part, these will be primitives---doubles mostly.
   * Passing these by reference is probably less efficient than
   * passing by value.  But, for user-defined types pass-by-reference
   * might be the way to go and the cost in this case will be much
   * higher than the value-reference trade-off for primitives.  Right
   * now we use pass-by-reference but we might reconsider...
   */

  /**** ARITHMETIC OPERATORS ****/

  /* These macros provide templates for the basic definitions required
   * for all of the binary operators.  Each operator requires 6
   * definitions.  First, a general matrix definition must be
   * provided.  This definition can return a matrix of a different
   * style and order than its arguments but can only be called if its
   * template type is explicitly specified.  The actual logic of the
   * operator should be specified within this function.  The macros
   * provide definitions for the other 5 required templates, one
   * default matrix by matrix, general matrix by scalar, default
   * matrix by scalar, general scalar by matrix, default scalar by
   * matrix.  The default versions call the more general versions with
   * such that they will return concrete matrices with order equal to
   * the left-hand (or only) matrix passed to the default version.
   *
   */

#define SCYTHE_BINARY_OPERATOR_DMM(OP)                                \
  template <matrix_order ORDER, matrix_style L_STYLE,                 \
            matrix_order R_ORDER, matrix_style R_STYLE,               \
            typename T_type>                                          \
  inline Matrix<T_type, ORDER, Concrete>                              \
  OP (const Matrix<T_type, ORDER, L_STYLE>& lhs,                      \
      const Matrix<T_type, R_ORDER, R_STYLE>& rhs)                    \
  {                                                                   \
    return OP <T_type, ORDER, Concrete>(lhs, rhs);                    \
  }

#define SCYTHE_BINARY_OPERATOR_GMS(OP)                                \
  template <typename T_type, matrix_order ORDER, matrix_style STYLE,  \
            matrix_order L_ORDER, matrix_style L_STYLE>               \
  inline Matrix<T_type, ORDER, STYLE>                                 \
  OP (const Matrix<T_type, L_ORDER, L_STYLE>& lhs,                    \
      const typename Matrix<T_type>::ttype &rhs)                       \
  {                                                                   \
    return  (OP <T_type, ORDER, STYLE>                                \
        (lhs, Matrix<T_type, L_ORDER>(rhs)));                         \
  }

#define SCYTHE_BINARY_OPERATOR_DMS(OP)                                \
  template <matrix_order ORDER, matrix_style L_STYLE,                 \
            typename T_type>                                          \
  inline Matrix<T_type, ORDER, Concrete>                              \
  OP (const Matrix<T_type, ORDER, L_STYLE>& lhs,                      \
      const typename Matrix<T_type>::ttype &rhs)                      \
  {                                                                   \
    return (OP <T_type, ORDER, Concrete> (lhs, rhs));                 \
  }
  
#define SCYTHE_BINARY_OPERATOR_GSM(OP)                                \
  template <typename T_type, matrix_order ORDER, matrix_style STYLE,  \
            matrix_order R_ORDER, matrix_style R_STYLE>               \
  inline Matrix<T_type, ORDER, STYLE>                                 \
  OP (const typename Matrix<T_type>::ttype &lhs,                      \
      const Matrix<T_type, R_ORDER, R_STYLE>& rhs) \
  {                                                                   \
    return  (OP <T_type, ORDER, STYLE>                                \
        (Matrix<T_type, R_ORDER>(lhs), rhs));                         \
  }

#define SCYTHE_BINARY_OPERATOR_DSM(OP)                                \
  template <matrix_order ORDER, matrix_style R_STYLE,                 \
            typename T_type>                                          \
  inline Matrix<T_type, ORDER, Concrete>                              \
  OP (const typename Matrix<T_type>::ttype &lhs,                      \
      const Matrix<T_type, ORDER, R_STYLE>& rhs)                      \
  {                                                                   \
    return (OP <T_type, ORDER, Concrete> (lhs, rhs));                 \
  }

#define SCYTHE_BINARY_OPERATOR_DEFS(OP)                               \
  SCYTHE_BINARY_OPERATOR_DMM(OP)                                      \
  SCYTHE_BINARY_OPERATOR_GMS(OP)                                      \
  SCYTHE_BINARY_OPERATOR_DMS(OP)                                      \
  SCYTHE_BINARY_OPERATOR_GSM(OP)                                      \
  SCYTHE_BINARY_OPERATOR_DSM(OP)

  /* Matrix multiplication */
  
  /* General template version. Must be called with operator*<> syntax
   */
 
  /* We provide two symmetric algorithms for matrix multiplication,
   * one for col-major and the other for row-major matrices.  They are
   * designed to minimize cache misses.The decision is based on the
   * return type of the template so, when using matrices of multiple
   * orders, this can get ugly.  These optimizations only really start
   * paying dividends as matrices get big, because cache misses are
   * rare with smaller matrices.
   */

  template <typename T_type, matrix_order ORDER, matrix_style STYLE,
            matrix_order L_ORDER, matrix_style L_STYLE,
            matrix_order R_ORDER, matrix_style R_STYLE>
  inline Matrix<T_type, ORDER, STYLE>
  operator* (const Matrix<T_type, L_ORDER, L_STYLE>& lhs,
             const Matrix<T_type, R_ORDER, R_STYLE>& rhs)
  {
    if (lhs.size() == 1 || rhs.size() == 1)
      return (lhs % rhs);

    SCYTHE_CHECK_10 (lhs.cols() != rhs.rows(),
        scythe_conformation_error,
        "Matrices with dimensions (" << lhs.rows() 
        << ", " << lhs.cols()
        << ") and (" << rhs.rows() << ", " << rhs.cols()
        << ") are not multiplication-conformable");

    Matrix<T_type, ORDER, Concrete> result (lhs.rows(), rhs.cols(), false);

    T_type tmp;
    if (ORDER == Col) { // col-major optimized
     for (uint j = 0; j < rhs.cols(); ++j) {
       for (uint i = 0; i < lhs.rows(); ++i)
        result(i, j) = (T_type) 0;
       for (uint l = 0; l < lhs.cols(); ++l) {
         tmp = rhs(l, j);
         for (uint i = 0; i < lhs.rows(); ++i)
           result(i, j) += tmp * lhs(i, l);
       }
     }
    } else { // row-major optimized
     for (uint i = 0; i < lhs.rows(); ++i) {
       for (uint j = 0; j < rhs.cols(); ++j)
         result(i, j) = (T_type) 0;
       for (uint l = 0; l < rhs.rows(); ++l) {
         tmp = lhs(i, l);
         for (uint j = 0; j < rhs.cols(); ++j)
           result(i, j) += tmp * rhs(l,j);
       }
     }
    }

    SCYTHE_VIEW_RETURN(T_type, ORDER, STYLE, result)
  }

  SCYTHE_BINARY_OPERATOR_DEFS(operator*)

  
  template <typename T_type, matrix_order ORDER, matrix_style STYLE,
            matrix_order L_ORDER, matrix_style L_STYLE,
            matrix_order R_ORDER, matrix_style R_STYLE>
  inline Matrix<T_type, ORDER, STYLE>
  kronecker (const Matrix<T_type, L_ORDER, L_STYLE>& lhs,
             const Matrix<T_type, R_ORDER, R_STYLE>& rhs)
  {
    Matrix<T_type,ORDER,Concrete> res = lhs;
    res.kronecker(rhs);
    return (res);
  }

  SCYTHE_BINARY_OPERATOR_DEFS(kronecker)

  /* Macro definition for general return type templates of standard
   * binary operators (this handles, +, -, %, /, but not *)
   */
    
#define SCYTHE_GENERAL_BINARY_OPERATOR(OP,FUNCTOR)                    \
  template <typename T_type, matrix_order ORDER, matrix_style STYLE,  \
            matrix_order L_ORDER, matrix_style L_STYLE,               \
            matrix_order R_ORDER, matrix_style R_STYLE>               \
  inline Matrix<T_type, ORDER, STYLE>                                 \
  OP (const Matrix<T_type, L_ORDER, L_STYLE>& lhs,                    \
      const Matrix<T_type, R_ORDER, R_STYLE>& rhs)                    \
  {                                                                   \
    SCYTHE_CHECK_10(lhs.size() != 1 && rhs.size() != 1 &&             \
        (lhs.rows() != rhs.rows() || lhs.cols() != rhs.cols()),       \
        scythe_conformation_error,                                    \
        "Matrices with dimensions (" << lhs.rows()                    \
        << ", " << lhs.cols()                                         \
        << ") and (" << rhs.rows() << ", " << rhs.cols()              \
        << ") are not conformable");                                  \
                                                                      \
    if (lhs.size() == 1) {                                            \
      Matrix<T_type,ORDER,Concrete> res(rhs.rows(),rhs.cols(),false); \
      std::transform(rhs.begin_f(), rhs.end_f(),                      \
          res.template begin_f<R_ORDER>(),                            \
          std::bind1st(FUNCTOR <T_type>(), lhs(0)));                  \
      SCYTHE_VIEW_RETURN(T_type, ORDER, STYLE, res)                   \
    }                                                                 \
                                                                      \
    Matrix<T_type,ORDER,Concrete> res(lhs.rows(), lhs.cols(), false); \
                                                                      \
    if (rhs.size() == 1) {                                            \
      std::transform(lhs.begin_f(), lhs.end_f(),                      \
          res.template begin_f<L_ORDER> (),                           \
          std::bind2nd(FUNCTOR <T_type>(), rhs(0)));                  \
    } else {                                                          \
      std::transform(lhs.begin_f(), lhs.end_f(),                      \
          rhs.template begin_f<L_ORDER>(),                            \
          res.template begin_f<L_ORDER>(),                            \
          FUNCTOR <T_type> ());                                       \
    }                                                                 \
                                                                      \
    SCYTHE_VIEW_RETURN(T_type, ORDER, STYLE, res)                     \
  }

  /* Addition operators */

  SCYTHE_GENERAL_BINARY_OPERATOR (operator+, std::plus)
  SCYTHE_BINARY_OPERATOR_DEFS (operator+)

  /* Subtraction operators */

  SCYTHE_GENERAL_BINARY_OPERATOR (operator-, std::minus)
  SCYTHE_BINARY_OPERATOR_DEFS (operator-)

  /* Element-by-element multiplication operators */

  SCYTHE_GENERAL_BINARY_OPERATOR (operator%, std::multiplies)
  SCYTHE_BINARY_OPERATOR_DEFS(operator%)

  /* Element-by-element division */

  SCYTHE_GENERAL_BINARY_OPERATOR (operator/, std::divides)
  SCYTHE_BINARY_OPERATOR_DEFS (operator/)

  /* Element-by-element exponentiation */

  SCYTHE_GENERAL_BINARY_OPERATOR (operator^, exponentiate)
  SCYTHE_BINARY_OPERATOR_DEFS (operator^)

  /* Negation operators */

  // General return type version
  template <typename T_type, matrix_order R_ORDER, matrix_style R_STYLE,
            matrix_order ORDER, matrix_style STYLE>
  inline Matrix<T_type, R_ORDER, R_STYLE>
  operator- (const Matrix<T_type, ORDER, STYLE>& M)
  {
    Matrix<T_type, R_ORDER, Concrete> result(M.rows(), M.cols(), false);
    std::transform(M.template begin_f<ORDER>(), 
                   M.template end_f<ORDER>(), 
                   result.template begin_f<R_ORDER>(),
                   std::negate<T_type> ());
    SCYTHE_VIEW_RETURN(T_type, R_ORDER, R_STYLE, result)
  }
  
  // Default return type version
  template <matrix_order ORDER, matrix_style P_STYLE, typename T_type>
  inline Matrix<T_type, ORDER, Concrete>
  operator- (const Matrix<T_type, ORDER, P_STYLE>& M)
  {
    return operator-<T_type, ORDER, Concrete> (M);
  }

  /* Unary not operators */

  template <matrix_order R_ORDER, matrix_style R_STYLE, typename T_type,
            matrix_order ORDER, matrix_style STYLE>
  inline Matrix<bool, R_ORDER, R_STYLE>
  operator! (const Matrix<T_type, ORDER, STYLE>& M)
  {
    Matrix<bool, R_ORDER, Concrete> result(M.rows(), M.cols(), false);
    std::transform(M.template begin_f<ORDER>(), 
                   M.template end_f<ORDER>(), 
                   result.template begin_f<R_ORDER>(),
                   std::logical_not<T_type> ());
    SCYTHE_VIEW_RETURN(T_type, R_ORDER, R_STYLE, result)
  }
  
  // Default return type version
  template <typename T_type, matrix_order ORDER, matrix_style P_STYLE>
  inline Matrix<bool, ORDER, Concrete>
  operator! (const Matrix<T_type, ORDER, P_STYLE>& M)
  {
    return (operator!<ORDER, Concrete> (M));
  }
  /**** COMPARISON OPERATORS ****/

  /* These macros are analogous to those above, except they return
   * only boolean matrices and use slightly different template
   * parameter orderings.  Kind of redundant, but less confusing than
   * making omnibus macros that handle both cases.
   */
#define SCYTHE_GENERAL_BINARY_BOOL_OPERATOR(OP,FUNCTOR)               \
  template <matrix_order ORDER, matrix_style STYLE, typename T_type,  \
            matrix_order L_ORDER, matrix_style L_STYLE,               \
            matrix_order R_ORDER, matrix_style R_STYLE>               \
  inline Matrix<bool, ORDER, STYLE>                                   \
  OP (const Matrix<T_type, L_ORDER, L_STYLE>& lhs,                    \
      const Matrix<T_type, R_ORDER, R_STYLE>& rhs)                    \
  {                                                                   \
    SCYTHE_CHECK_10(lhs.size() != 1 && rhs.size() != 1 &&             \
        (lhs.rows() != rhs.rows() || lhs.cols() != rhs.cols()),       \
        scythe_conformation_error,                                    \
        "Matrices with dimensions (" << lhs.rows()                    \
        << ", " << lhs.cols()                                         \
        << ") and (" << rhs.rows() << ", " << rhs.cols()              \
        << ") are not conformable");                                  \
                                                                      \
    if (lhs.size() == 1) {                                            \
      Matrix<bool,ORDER,Concrete> res(rhs.rows(),rhs.cols(),false);   \
      std::transform(rhs.begin_f(), rhs.end_f(),                      \
          res.template begin_f<R_ORDER>(),                            \
          std::bind1st(FUNCTOR <T_type>(), lhs(0)));                  \
      SCYTHE_VIEW_RETURN(T_type, ORDER, STYLE, res)                   \
    }                                                                 \
                                                                      \
    Matrix<bool,ORDER,Concrete> res(lhs.rows(), lhs.cols(), false);   \
                                                                      \
    if (rhs.size() == 1) {                                            \
      std::transform(lhs.begin_f(), lhs.end_f(),                      \
          res.template begin_f<L_ORDER> (),                           \
          std::bind2nd(FUNCTOR <T_type>(), rhs(0)));                  \
    } else {                                                          \
      std::transform(lhs.begin_f(), lhs.end_f(),                      \
          rhs.template begin_f<L_ORDER>(),                            \
          res.template begin_f<L_ORDER>(),                            \
          FUNCTOR <T_type> ());                                       \
    }                                                                 \
                                                                      \
    SCYTHE_VIEW_RETURN(T_type, ORDER, STYLE, res)                     \
  }

#define SCYTHE_BINARY_BOOL_OPERATOR_DMM(OP)                           \
  template <typename T_type, matrix_order ORDER, matrix_style L_STYLE,\
            matrix_order R_ORDER, matrix_style R_STYLE>               \
  inline Matrix<bool, ORDER, Concrete>                                \
  OP (const Matrix<T_type, ORDER, L_STYLE>& lhs,                      \
             const Matrix<T_type, R_ORDER, R_STYLE>& rhs)             \
  {                                                                   \
    return OP <ORDER, Concrete>(lhs, rhs);                            \
  }

#define SCYTHE_BINARY_BOOL_OPERATOR_GMS(OP)                           \
  template <matrix_order ORDER, matrix_style STYLE, typename T_type,  \
            matrix_order L_ORDER, matrix_style L_STYLE>               \
  inline Matrix<bool, ORDER, STYLE>                                   \
  OP (const Matrix<T_type, L_ORDER, L_STYLE>& lhs,                    \
      const typename Matrix<T_type>::ttype &rhs)                      \
  {                                                                   \
    return  (OP <ORDER, STYLE>                                        \
        (lhs, Matrix<T_type, L_ORDER>(rhs)));                         \
  }

#define SCYTHE_BINARY_BOOL_OPERATOR_DMS(OP)                           \
  template <typename T_type, matrix_order ORDER, matrix_style L_STYLE>\
  inline Matrix<bool, ORDER, Concrete>                                \
  OP (const Matrix<T_type, ORDER, L_STYLE>& lhs,                      \
      const typename Matrix<T_type>::ttype &rhs)                      \
  {                                                                   \
    return (OP <ORDER, Concrete> (lhs, rhs));                         \
  }
  
#define SCYTHE_BINARY_BOOL_OPERATOR_GSM(OP)                           \
  template <matrix_order ORDER, matrix_style STYLE, typename T_type,  \
            matrix_order R_ORDER, matrix_style R_STYLE>               \
  inline Matrix<bool, ORDER, STYLE>                                   \
  OP (const typename Matrix<T_type>::ttype &lhs,                      \
      const Matrix<T_type, R_ORDER, R_STYLE>& rhs)                    \
  {                                                                   \
    return  (OP <ORDER, STYLE>                                        \
        (Matrix<T_type, R_ORDER>(lhs), rhs));                         \
  }

#define SCYTHE_BINARY_BOOL_OPERATOR_DSM(OP)                           \
  template <typename T_type, matrix_order ORDER, matrix_style R_STYLE>\
  inline Matrix<bool, ORDER, Concrete>                                \
  OP (const typename Matrix<T_type>::ttype &lhs,                      \
      const Matrix<T_type, ORDER, R_STYLE>& rhs)                      \
  {                                                                   \
    return (OP <ORDER, Concrete> (lhs, rhs));                         \
  }

#define SCYTHE_BINARY_BOOL_OPERATOR_DEFS(OP)                          \
  SCYTHE_BINARY_BOOL_OPERATOR_DMM(OP)                                 \
  SCYTHE_BINARY_BOOL_OPERATOR_GMS(OP)                                 \
  SCYTHE_BINARY_BOOL_OPERATOR_DMS(OP)                                 \
  SCYTHE_BINARY_BOOL_OPERATOR_GSM(OP)                                 \
  SCYTHE_BINARY_BOOL_OPERATOR_DSM(OP)

  /* Element-wise Equality operator
   * See equals() method for straight equality checks
   */

  SCYTHE_GENERAL_BINARY_BOOL_OPERATOR (operator==, std::equal_to)
  SCYTHE_BINARY_BOOL_OPERATOR_DEFS (operator==)

  SCYTHE_GENERAL_BINARY_BOOL_OPERATOR (operator!=, std::not_equal_to)
  SCYTHE_BINARY_BOOL_OPERATOR_DEFS (operator!=)

  SCYTHE_GENERAL_BINARY_BOOL_OPERATOR (operator<, std::less)
  SCYTHE_BINARY_BOOL_OPERATOR_DEFS (operator<)

  SCYTHE_GENERAL_BINARY_BOOL_OPERATOR (operator<=, std::less_equal)
  SCYTHE_BINARY_BOOL_OPERATOR_DEFS (operator<=)

  SCYTHE_GENERAL_BINARY_BOOL_OPERATOR (operator>, std::greater)
  SCYTHE_BINARY_BOOL_OPERATOR_DEFS (operator>)

  SCYTHE_GENERAL_BINARY_BOOL_OPERATOR (operator>=, std::greater_equal)
  SCYTHE_BINARY_BOOL_OPERATOR_DEFS (operator>=)

  SCYTHE_GENERAL_BINARY_BOOL_OPERATOR (operator&, std::logical_and)
  SCYTHE_BINARY_BOOL_OPERATOR_DEFS (operator&)


  SCYTHE_GENERAL_BINARY_BOOL_OPERATOR (operator|, std::logical_or)
  SCYTHE_BINARY_BOOL_OPERATOR_DEFS (operator|)

  /**** INPUT-OUTPUT ****/


  /* This function simply copies values from an input stream into a
   * matrix.  It relies on the iterators for bounds checking.
   */

  template <typename T, matrix_order O, matrix_style S>
  std::istream& operator>> (std::istream& is, Matrix<T,O,S>& M)
  {
    std::copy(std::istream_iterator<T> (is), std::istream_iterator<T>(),
         M.begin_f());

    return is;
  }

  /* Writes a matrix to an ostream in readable format.  This is
   * intended to be used to pretty-print to the terminal.
   */

  template <typename T, matrix_order O, matrix_style S>
  std::ostream& operator<< (std::ostream& os, const Matrix<T,O,S>& M)
  {
    /* This function take two passes to figure out appropriate field
     * widths.  Speed isn't really the point here.
     */

    // Store previous io settings
    std::ios_base::fmtflags preop = os.flags();

    uint mlen = os.width();
    std::ostringstream oss;
    oss.precision(os.precision());
    oss << std::setiosflags(std::ios::fixed);
    
    typename Matrix<T,O,S>::const_forward_iterator last = M.end_f();
    for (typename Matrix<T,O,S>::const_forward_iterator i = M.begin_f();
        i != last; ++i) {
      oss.str("");
      oss << (*i);
      if (oss.str().length() > mlen)
        mlen = oss.str().length();
    }


    /* Write the stream */
    // Change to a fixed with format.  Users should control precision
    os << std::setiosflags(std::ios::fixed);

    
    for (uint i = 0; i < M.rows(); ++i) {
      Matrix<T, O, View> row = M(i, _);
      //for (uint i = 0; i < row.size(); ++i)
      //  os << std::setw(mlen) << row[i] << " ";
      typename Matrix<T,O,View>::const_forward_iterator row_last 
        = row.end_f();
      for (typename 
          Matrix<T,O,View>::forward_iterator el = row.begin_f();
          el != row_last; ++el) {
        os << std::setw(mlen) << *el << " ";
      }
      os << std::endl;
    }
    
    
    // Restore pre-op flags
    os.flags(preop);

    return os;
  }

#ifdef SCYTHE_LAPACK
  /* A template specialization of operator* for col-major, concrete
   * matrices of doubles that is only visible when a LAPACK library is
   * available.  This function is an analog of the above function and
   * the above doxygen documentation serves for both.
   *
   * This needs to go below % so it can see the template definition
   * (since it isn't actually in the template itself.
   */

  template<>
  inline Matrix<>
  operator*<double,Col,Concrete,Col,Concrete>
  (const Matrix<>& lhs, const Matrix<>& rhs)
  {
    if (lhs.size() == 1 || rhs.size() == 1)
      return (lhs % rhs);

    SCYTHE_DEBUG_MSG("Using lapack/blas for matrix multiplication");
    SCYTHE_CHECK_10 (lhs.cols() != rhs.rows(),
        scythe_conformation_error,
        "Matrices with dimensions (" << lhs.rows() 
        << ", " << lhs.cols()
        << ") and (" << rhs.rows() << ", " << rhs.cols()
        << ") are not multiplication-conformable");

    Matrix<> result (lhs.rows(), rhs.cols(), false);

    // Get pointers to the internal arrays and set up some vars
    double* lhspnt = lhs.getArray();
    double* rhspnt = rhs.getArray();
    double* resultpnt = result.getArray();
    const double one(1.0);
    const double zero(0.0);
    int rows = (int) lhs.rows();
    int cols = (int) rhs.cols();
    int innerDim = (int) rhs.rows();

    // Call the lapack routine.
    lapack::dgemm_("N", "N", &rows, &cols, &innerDim, &one, lhspnt,
                   &rows, rhspnt, &innerDim, &zero, resultpnt, &rows);

    return result;
  }
#endif

} // end namespace scythe

#endif /* SCYTHE_MATRIX_H */