/*
 Copyright (C) 2014 Eugenio Gianniti <eugenio.gianniti@mail.polimi.it>

 This program 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.

 This program 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 along with
 this program; if not, see <http://www.gnu.org/licenses/>.
*/

#include "uBLAS_factory.h"

femfenics_base_factory const &
uBLAS_factory::instance (void)
{
  static uBLAS_factory const theinstance;
  return theinstance;
}

octave_value
uBLAS_factory::matrix (dolfin::Matrix const & A) const
{
  octave_value retval;

  // Get capacity of the dolfin sparse matrix
  boost::tuples::tuple <const std::size_t *,
                        const std::size_t *,
                        const double *, int>
    aa = A.data ();

  int nnz = aa.get<3> ();
  std::size_t nr = A.size (0), nc = A.size (1);
  std::vector<double> data_tmp;
  std::vector<std::size_t> cidx_tmp;

  dim_vector dims (nnz, 1);
  octave_idx_type nz = 0, ii = 0;
  Array<octave_idx_type> ridx (dims, 0), cidx (dims, 0);
  Array<double> data (dims, 0);

  octave_idx_type * orow = ridx.fortran_vec ();
  octave_idx_type * oc = cidx.fortran_vec ();
  double * ov = data.fortran_vec ();

  for (std::size_t i = 0; i < nr; ++i)
    {
      A.getrow (i, cidx_tmp, data_tmp);
      nz += cidx_tmp.size ();

      for (octave_idx_type j = 0;
           j < cidx_tmp.size (); ++j)
        {
          orow[ii + j] = i;
          oc[ii + j] = cidx_tmp[j];
          ov[ii + j] = data_tmp[j];
        }

      ii = nz;
    }

  dims(0) = ii;
  ridx.resize (dims);
  cidx.resize (dims);
  data.resize (dims);

  SparseMatrix sm (data, ridx, cidx, nr, nc);
  retval = sm;

  return retval;
}

octave_value
uBLAS_factory::vector (dolfin::Vector const & b) const
{
  octave_value retval;

  dim_vector dims;
  dims.resize (2);
  dims(0) = b.size ();
  dims(1) = 1;
  Array<double> myb (dims);

  for (std::size_t i = 0; i < b.size (); ++i)
    { myb.xelem (i) = b[i]; }

  retval = myb;

  return retval;
}