$$\phantom{******** If you see this on the webpage then the browser could not locate *********}$$
$$\phantom{******** jsMath/easy/load.js or the variable root is set wrong in this file *********}$$
$$\newcommand{\vector}[1]{\left[\begin{array}{c} #1 \end{array}\right]}$$ $$\newenvironment{matrix}{\left[\begin{array}{cccccccccc}} {\end{array}\right]}$$ $$\newcommand{\A}{{\cal A}}$$ $$\newcommand{\W}{{\cal W}}$$

/* Copyright 2008-2014 Research Foundation State University of New York   */

/* This file is part of QUB Express.                                      */

/* QUB Express is free software; you can redistribute it and/or modify    */
/* it under the terms of the GNU General Public License as published by   */
/* the Free Software Foundation, either version 3 of the License, or      */
/* (at your option) any later version.                                    */

/* QUB Express is distributed in the hope that it will be useful,         */
/* but WITHOUT ANY WARRANTY; without even the implied warranty of         */
/* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the          */
/* GNU General Public License for more details.                           */

/* You should have received a copy of the GNU General Public License,     */
/* named LICENSE.txt, in the QUB Express program directory.  If not, see  */
/* <http://www.gnu.org/licenses/>.                                        */

#ifndef UBLAS_PLUS_H
#define UBLAS_PLUS_H

/*

*/

// This header includes the vector and matrix from boost::ublas,
// and defines template functions to set up, read, print, and store them.

#include <boost/numeric/ublas/vector.hpp>
#include <boost/numeric/ublas/matrix.hpp>
#include <algorithm>
#include <vector>
#include "qubfast.h"

using namespace boost::numeric::ublas;
using std::min;
using std::max;

template<class S, class T>
 void matrix_from_arrays(matrix<S>& mat, T** dat, int m, int n)
{
  if ( mat.size1() != (size_t) m || mat.size2() != (size_t) n )
    mat.resize(m, n);
  for (int i=0; i<m; ++i) {
    S* mat_i = & mat(i, 0);
    T* dat_i = dat[i];
    for (int j=n; j; --j)
      *(mat_i++) = (S) *(dat_i++);
  }
  /* some compilers can't handle these templates:
  // typename M::iterator1; // so the compiler doesn't treat it as R-value
  typename M::iterator1 mat_i = mat.begin1();
  for (int i=0; i<m; ++i, ++mat_i) {
#ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION
    typename M::iterator2 mat_j = mat_i.begin();
#else
    typename M::iterator2 mat_j = begin(mat_i, iterator1_tag());
#endif
    T* dat_j = dat[i];
    for (int j=0; j<n; ++j)
      *(mat_j++) = *(dat_j++);
  }
  */
}

template<class S, class T>
void matrix_to_arrays(matrix<S>& mat, T** dat, int m, int n)
{
  for (int i=0; i<m; ++i) {
    S* mat_i = & mat(i, 0);
    T* dat_i = dat[i];
    for (int j=n; j; --j)
      *(dat_i++) = (T) *(mat_i++);
  }
  /* some compilers can't handle these templates:
  typename M::iterator1 mat_i = mat.begin1();
  for (int i=0; i<m; ++i, ++mat_i) {
#ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION
    typename M::iterator2 mat_j = mat_i.begin();
#else
    typename M::iterator2 mat_j = begin(mat_i, iterator1_tag());
#endif
    T* dat_j = dat[i];
    for (int j=0; j<n; ++j)
      *(dat_j++) = *(mat_j++);
  }
  */
}

template<class S, class T>
void matrix_from_array(matrix<S>& mat, T* dat, int m, int n)
{
  if ( mat.size1() != (size_t) m || mat.size2() != (size_t) n )
    mat.resize(m, n);
  T* dat_i = dat;
  for (int i=0; i<m; ++i) {
    S* mat_i = & mat(i, 0);
    for (int j=n; j; --j)
      *(mat_i++) = (S) *(dat_i++);
  }
  // typename M::iterator1; // so the compiler doesn't treat it as R-value
  /* some compilers can't handle these templates:
  typename M::iterator1 mat_i = mat.begin1();
  for (int i=0; i<m; ++i, ++mat_i) {
#ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION
    typename M::iterator2 mat_j = mat_i.begin();
#else
    typename M::iterator2 mat_j = begin(mat_i, iterator1_tag());
#endif
    for (int j=0; j<n; ++j) {
      *(mat_j++) = *(dat_i++);
    }
  }
  */
}

template<class S, class T>
  void matrix_to_array(matrix<S>& mat, T* dat, int m, int n)
{
  int M = min(m, (int) mat.size1());
  int N = min(n, (int) mat.size2());
  int i, j;
  for (i=0; i<M; ++i) {
    S *mat_i = & mat(i, 0);
    for (j=0; j<N; ++j)
      *(dat++) = (T) *(mat_i++);
    for (; j<n; ++j)
      *(dat++) = 0.0;
  }
  for (; i<m; ++i)
    for (j=0; j<m; ++j)
      *(dat++) = 0.0;
}

template<class T>
void matrix_to_partition(matrix<T>& mat, int m, int n, int *rowix, int *colix, matrix<T>& part)
{
  if ( part.size1() != (size_t) m || part.size2() != (size_t) n )
    part.resize(m, n);
  for (int i=0; i<m; ++i)
    for (int j=0; j<n; ++j)
      part(i,j) = mat(rowix[i], colix[j]);
}

template<class T>
void matrix_from_partition(matrix<T>& mat, int m, int n, int *rowix, int *colix, matrix<T>& part)
{
  if ( part.size1() != (size_t) m || part.size2() != (size_t) n )
    part.resize(m, n);
  for (int i=0; i<m; ++i)
    for (int j=0; j<n; ++j)
      mat(rowix[i], colix[j]) = part(i,j);
}

template<class T>
void fill_identity(matrix<T>& mat)
{
  int m = (int) mat.size1();
  int n = (int) mat.size2();
  for (int i=0; i<m; ++i)
    for (int j=0; j<n; ++j)
      mat(i,j) = (i==j) ? 1 : 0;
}

template<class T>
void fill_const(matrix<T>& mat, T val)
{
  for (typename matrix<T>::iterator1 ii = mat.begin1(); ii != mat.end1(); ++ii)
#ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION
    for (typename matrix<T>::iterator2 jj = ii.begin(); jj != ii.end(); ++jj)
#else
    for (typename matrix<T>::iterator2 jj = begin(ii, iterator1_tag()); jj != end(ii, iterator1_tag()); ++jj)
#endif
      *jj = val;
}

template<class T>
std::ostream& operator<<(std::ostream& ost, vector<T>& vec)
{
  ost << "[ ";
  for (typename vector<T>::iterator ii = vec.begin(); ii != vec.end(); ++ii)
    ost << *ii << ' ';
  ost << "]";
  return ost;
}

template<class T>
std::ostream& operator<<(std::ostream& ost, matrix<T>& mat)
{
  ost << "[ ";
  for (typename matrix<T>::iterator1 ii = mat.begin1(); ii != mat.end1(); ++ii) {
    ost << "[ ";
#ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION
    for (typename matrix<T>::iterator2 jj = ii.begin(); jj != ii.end(); ++jj) {
#else
    for (typename matrix<T>::iterator2 jj = begin(ii, iterator1_tag()); jj != end(ii, iterator1_tag()); ++jj) {
#endif
      ost << *jj << ' ';
    }
    ost << "]" << std::endl << "  ";
  }
  ost << "]";
  return ost;
}

  // These are to simplify working with C-interfaced code, which variously expects T* or T**
  // e.g. in max_mac_ll.cpp

template<class T> class pointy_matrix {
public:
  matrix<T> m;
  std::vector<T*> rows;
  T *p, **pp;
  int n; // == m.size1() * m.size2()
  
  pointy_matrix()
    : p(NULL), pp(NULL), n(0)
  {}
  pointy_matrix(int nr, int nc)
    : m(nr, nc), rows(nr), p(NULL), pp(NULL), n(nr*nc)
  {
    resize(nr, nc);
  }
  void resize(int nr, int nc)
  {
    n = nr*nc;
    m.resize(nr, nc);
    m.clear();
    rows.resize(nr);
    for (int i=0; i<nr; ++i)
      rows[i] = &(m(i,0));
    if ( nr && nc ) {
      p = &(m(0,0));
      pp = &(rows[0]);
    }
  }
  void clear()
  {
    m.clear();
  }
};

template<class T>
class pointy_col_vector {
public:
  matrix<T> m;
  T *p, **pp;
  int n;

  pointy_col_vector()
    : p(NULL), pp(&p), n(0)
  {}
  pointy_col_vector(int _n)
    : m(_n,1), p(NULL), pp(&p), n(_n)
  {
    resize(_n);
  }
  pointy_col_vector(int _n, T *data)
    : m(_n,1), p(NULL), pp(&p), n(_n)
  {
    resize(_n);
    memcpy(p, data, _n*sizeof(T));
  }
  void resize(int _n)
  {
    n = _n;
    m.resize(_n,1);
    m.clear();
    if ( _n )
      p = &(m(0,0));
  }
  void clear()
  {
    m.clear();
  }
};


#ifdef _WIN32
QUBFAST_TEMPLATE template class QUBFAST_API boost::numeric::ublas::vector<char>;
QUBFAST_TEMPLATE template class QUBFAST_API boost::numeric::ublas::matrix<char>;
QUBFAST_TEMPLATE template class QUBFAST_API boost::numeric::ublas::vector<short>;
QUBFAST_TEMPLATE template class QUBFAST_API boost::numeric::ublas::matrix<short>;
QUBFAST_TEMPLATE template class QUBFAST_API boost::numeric::ublas::vector<int>;
QUBFAST_TEMPLATE template class QUBFAST_API boost::numeric::ublas::matrix<int>;
QUBFAST_TEMPLATE template class QUBFAST_API boost::numeric::ublas::vector<float>;
QUBFAST_TEMPLATE template class QUBFAST_API boost::numeric::ublas::matrix<float>;
QUBFAST_TEMPLATE template class QUBFAST_API boost::numeric::ublas::vector<double>;
QUBFAST_TEMPLATE template class QUBFAST_API boost::numeric::ublas::matrix<double>;

QUBFAST_TEMPLATE template class QUBFAST_API pointy_matrix<char>;
QUBFAST_TEMPLATE template class QUBFAST_API pointy_matrix<short>;
QUBFAST_TEMPLATE template class QUBFAST_API pointy_matrix<int>;
QUBFAST_TEMPLATE template class QUBFAST_API pointy_matrix<float>;
QUBFAST_TEMPLATE template class QUBFAST_API pointy_matrix<double>;
QUBFAST_TEMPLATE template class QUBFAST_API pointy_col_vector<char>;
QUBFAST_TEMPLATE template class QUBFAST_API pointy_col_vector<short>;
QUBFAST_TEMPLATE template class QUBFAST_API pointy_col_vector<int>;
QUBFAST_TEMPLATE template class QUBFAST_API pointy_col_vector<float>;
QUBFAST_TEMPLATE template class QUBFAST_API pointy_col_vector<double>;

#endif

  
#endif