/*
 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 "PETSc_factory.h"

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

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

    std::size_t nr = A.size (0), nc = A.size (1);

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

    for (std::size_t i = 0; i < nr; ++i)
      {
        std::vector <double> data_tmp;
        std::vector <std::size_t> cidx_tmp;
        A.getrow (i, cidx_tmp, data_tmp);
        std::vector <std::size_t> ridx_tmp (cidx_tmp.size (), i);
        add_to_arrays (ridx, cidx, data, ridx_tmp, cidx_tmp, data_tmp);
      }

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

    return retval;
  }

octave_value
PETSc_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;
  }

void
PETSc_factory::add_to_arrays (Array <octave_idx_type> & ridx,
                              Array <octave_idx_type> & cidx,
                              Array <double> & values,
                              std::vector <std::size_t> const& row,
                              std::vector <std::size_t> const& col,
                              std::vector <double> const& val)
  {
    octave_idx_type pos = cidx.numel ();

    dim_vector dims = cidx.dims ();
    if (dims(1) == 1)
      dims(0) += col.size ();
    else
      dims(1) += col.size ();
    ridx.resize (dims, 0);
    cidx.resize (dims, 0);
    values.resize (dims, 0);

    std::vector <std::size_t>::const_iterator rit = row.begin (),
                                                    cit = col.begin ();
    std::vector <double>::const_iterator vit = val.begin ();

    while (cit != col.end ())
      {
        ridx(pos) = *rit;
        cidx(pos) = *cit;
        values(pos) = *vit;

        ++rit;
        ++cit;
        ++vit;
        ++pos;
      }

    return;
  }