/*
 *  Copyright 2008-2009 NVIDIA Corporation
 *
 *  Licensed under the Apache License, Version 2.0 (the "License");
 *  you may not use this file except in compliance with the License.
 *  You may obtain a copy of the License at
 *
 *      http://www.apache.org/licenses/LICENSE-2.0
 *
 *  Unless required by applicable law or agreed to in writing, software
 *  distributed under the License is distributed on an "AS IS" BASIS,
 *  WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
 *  See the License for the specific language governing permissions and
 *  limitations under the License.
 */

/*
 * Modifications to this file:
 * Copyright (c) 2014-2015, The University of Queensland
 * Licensed under the Apache License, Version 2.0.
 *
 */



#pragma once

#include <cusp/cmath.h>
#include <cusp/ell_matrix.h>
#include <cusp/exception.h>

#include <cusp/detail/host/conversion_utils.h>

#include <thrust/fill.h>
#include <thrust/extrema.h>
#include <thrust/count.h>

namespace cusp
{
namespace detail
{
namespace host
{

/////////////////////
// CDS Conversions //
/////////////////////
    
template <typename Matrix1, typename Matrix2>
void cds_to_csr(const Matrix1& src, Matrix2& dst)
{
    //typedef typename Matrix2::index_type IndexType;
    typedef typename Matrix2::value_type ValueType;
    
    size_t num_entries = 0;
    size_t num_diagonals = src.diagonal_offsets.size();

    // count nonzero entries
    for (size_t i = 0; i < num_diagonals; i++)
    {
        num_entries += src.block_size*src.block_size*(src.num_rows-cusp::abs(src.diagonal_offsets[i])*src.block_size)/src.block_size;
    }
    
    dst.resize(src.num_rows, src.num_rows, num_entries);

    num_entries = 0;
    dst.row_offsets[0] = 0;

    // copy nonzero entries to CSR structure
    for(size_t i = 0; i < src.num_rows; i++)
    {
        for(size_t n = 0; n < num_diagonals; n++)
        {
            if (src.diagonal_offsets[n] < 0 && i < cusp::abs(src.diagonal_offsets[n])*src.block_size)
                continue;
            if (src.diagonal_offsets[n] > 0 && i >= src.num_rows-src.diagonal_offsets[n]*src.block_size)
                continue;

            for (size_t j = 0; j < src.block_size; j++)
            {
                const ValueType value = src.values(i, n*src.block_size+j);

                dst.column_indices[num_entries] =
                    (src.diagonal_offsets[n]+i/src.block_size)*src.block_size+j;
                dst.values[num_entries] = value;
                num_entries++;
            }
        }

        dst.row_offsets[i + 1] = num_entries;
    }
}

template <typename Matrix1, typename Matrix2>
void cds_to_dia(const Matrix1& src, Matrix2& dst)
{
    typedef typename Matrix2::index_type IndexType;
    typedef typename Matrix2::value_type ValueType;
    
    const size_t num_diagonalblocks = src.diagonal_offsets.size();

    // number of diagonals for DIA format
    const size_t num_diagonals = num_diagonalblocks*(2*src.block_size-1);
    cusp::array1d<IndexType,cusp::host_memory> diagonal_offsets(num_diagonals, 0);

    size_t num_entries = 0;

    // count entries (with zero fill-in) for DIA format
    // and determine diagonal offsets
    for (size_t i = 0, n = 0; i < num_diagonalblocks; i++)
    {
        for (IndexType j = 1-(IndexType)src.block_size; j < (IndexType)src.block_size; j++)
        {
            const IndexType offset = src.diagonal_offsets[i] * src.block_size + j;
            diagonal_offsets[n++] = offset;
            num_entries += src.num_rows - cusp::abs(offset);
        }
    }

    // prepare destination matrix
    dst.resize(src.num_rows, src.num_rows, num_entries, num_diagonals);
    cusp::copy(diagonal_offsets, dst.diagonal_offsets);
    thrust::fill(dst.values.values.begin(), dst.values.values.end(), ValueType(0));

    for (size_t i = 0; i < num_diagonalblocks; i++)
    {
        const IndexType k = src.diagonal_offsets[i]*src.block_size;
        const IndexType first_row = std::max(0, -k);
        const IndexType num_rows = src.num_rows - cusp::abs(k);
        // row index is identical for CDS and DIA
        for (IndexType row = first_row; row < first_row+num_rows; row++)
        {
            for (IndexType j = 0; j < src.block_size; j++)
            {
                // index of "column" in DIA structure
                const IndexType col = i*(2*src.block_size-1)+(src.block_size-1-row%src.block_size)+j;
                dst.values(row,col) = src.values(row, i*src.block_size+j);
            }
        }
    }
}

/////////////////////
// COO Conversions //
/////////////////////
    
template <typename Matrix1, typename Matrix2>
void coo_to_csr(const Matrix1& src, Matrix2& dst)
{
    typedef typename Matrix2::index_type IndexType;
    //typedef typename Matrix2::value_type ValueType;

    dst.resize(src.num_rows, src.num_cols, src.num_entries);
    
    // compute number of non-zero entries per row of A 
    thrust::fill(dst.row_offsets.begin(), dst.row_offsets.end(), IndexType(0));

    for (size_t n = 0; n < src.num_entries; n++)
        dst.row_offsets[src.row_indices[n]]++;

    // cumsum the num_entries per row to get dst.row_offsets[]
    IndexType cumsum = 0;
    for(size_t i = 0; i < src.num_rows; i++)
    {
        IndexType temp = dst.row_offsets[i];
        dst.row_offsets[i] = cumsum;
        cumsum += temp;
    }
    dst.row_offsets[src.num_rows] = cumsum;

    // write Aj,Ax into dst.column_indices,dst.values
    for(size_t n = 0; n < src.num_entries; n++)
    {
        IndexType row  = src.row_indices[n];
        IndexType dest = dst.row_offsets[row];

        dst.column_indices[dest] = src.column_indices[n];
        dst.values[dest]         = src.values[n];

        dst.row_offsets[row]++;
    }

    IndexType last = 0; 
    for(size_t i = 0; i <= src.num_rows; i++)
    {
        IndexType temp = dst.row_offsets[i];
        dst.row_offsets[i]  = last;
        last   = temp;
    }

    //csr may contain duplicates
}

template <typename Matrix1, typename Matrix2>
void coo_to_array2d(const Matrix1& src, Matrix2& dst)
{
    //typedef typename Matrix2::index_type IndexType;
    typedef typename Matrix2::value_type ValueType;

    dst.resize(src.num_rows, src.num_cols);

    thrust::fill(dst.values.begin(), dst.values.end(), ValueType(0));

    for(size_t n = 0; n < src.num_entries; n++)
        dst(src.row_indices[n], src.column_indices[n]) += src.values[n]; //sum duplicates
}

/////////////////////
// CSR Conversions //
/////////////////////

template <typename Matrix1, typename Matrix2>
void csr_to_cds(const Matrix1& src, Matrix2& dst,
                const size_t alignment = 32)
{
    typedef typename Matrix2::index_type IndexType;
    typedef typename Matrix2::value_type ValueType;

    // compute number of occupied diagonals and enumerate them
    size_t num_diagonals = 0;

    // block size is always 1 for this conversion
    size_t block_size = 1;

    cusp::array1d<IndexType,cusp::host_memory> diag_map(src.num_rows + src.num_cols, 0);

    for(size_t i = 0; i < src.num_rows; i++)
    {
        for(IndexType jj = src.row_offsets[i]; jj < src.row_offsets[i+1]; jj++)
        {
            size_t j         = src.column_indices[jj];
            size_t map_index = (src.num_rows - i) + j; //offset shifted by + num_rows

            if(diag_map[map_index] == 0)
            {
                diag_map[map_index] = 1;
                num_diagonals++;
            }
        }
    }
   

    // allocate CDS structure
    dst.resize(src.num_rows, src.num_entries, num_diagonals, block_size, alignment);

    // fill in diagonal_offsets array
    for(size_t n = 0, diag = 0; n < src.num_rows + src.num_cols; n++)
    {
        if(diag_map[n] == 1)
        {
            diag_map[n] = diag;
            dst.diagonal_offsets[diag] = (IndexType) n - (IndexType) src.num_rows;
            diag++;
        }
    }

    // fill in values array
    thrust::fill(dst.values.values.begin(), dst.values.values.end(), ValueType(0));

    for(size_t i = 0; i < src.num_rows; i++)
    {
        for(IndexType jj = src.row_offsets[i]; jj < src.row_offsets[i+1]; jj++)
        {
            size_t j = src.column_indices[jj];
            size_t map_index = (src.num_rows - i) + j; //offset shifted by + num_rows
            size_t diag = diag_map[map_index];
        
            dst.values(i, diag) = src.values[jj];
        }
    }
}

template <typename Matrix1, typename Matrix2>
void csr_to_coo(const Matrix1& src, Matrix2& dst)
{
    typedef typename Matrix2::index_type IndexType;
    //typedef typename Matrix2::value_type ValueType;

    dst.resize(src.num_rows, src.num_cols, src.num_entries);
   
    // TODO replace with offsets_to_indices
    for(size_t i = 0; i < src.num_rows; i++)
        for(IndexType jj = src.row_offsets[i]; jj < src.row_offsets[i + 1]; jj++)
            dst.row_indices[jj] = i;

    cusp::copy(src.column_indices, dst.column_indices);
    cusp::copy(src.values,         dst.values);
}

template <typename Matrix1, typename Matrix2>
void csr_to_dia(const Matrix1& src, Matrix2& dst,
                const size_t alignment = 32)
{
    typedef typename Matrix2::index_type IndexType;
    typedef typename Matrix2::value_type ValueType;

    // compute number of occupied diagonals and enumerate them
    size_t num_diagonals = 0;

    cusp::array1d<IndexType,cusp::host_memory> diag_map(src.num_rows + src.num_cols, 0);

    for(size_t i = 0; i < src.num_rows; i++)
    {
        for(IndexType jj = src.row_offsets[i]; jj < src.row_offsets[i+1]; jj++)
        {
            size_t j         = src.column_indices[jj];
            size_t map_index = (src.num_rows - i) + j; //offset shifted by + num_rows

            if(diag_map[map_index] == 0)
            {
                diag_map[map_index] = 1;
                num_diagonals++;
            }
        }
    }
   

    // allocate DIA structure
    dst.resize(src.num_rows, src.num_cols, src.num_entries, num_diagonals, alignment);

    // fill in diagonal_offsets array
    for(size_t n = 0, diag = 0; n < src.num_rows + src.num_cols; n++)
    {
        if(diag_map[n] == 1)
        {
            diag_map[n] = diag;
            dst.diagonal_offsets[diag] = (IndexType) n - (IndexType) src.num_rows;
            diag++;
        }
    }

    // fill in values array
    thrust::fill(dst.values.values.begin(), dst.values.values.end(), ValueType(0));

    for(size_t i = 0; i < src.num_rows; i++)
    {
        for(IndexType jj = src.row_offsets[i]; jj < src.row_offsets[i+1]; jj++)
        {
            size_t j = src.column_indices[jj];
            size_t map_index = (src.num_rows - i) + j; //offset shifted by + num_rows
            size_t diag = diag_map[map_index];
        
            dst.values(i, diag) = src.values[jj];
        }
    }
}

template <typename Matrix1, typename Matrix2>
void csr_to_hyb(const Matrix1& src, Matrix2& dst,
                const size_t num_entries_per_row,
                const size_t alignment = 32)
{
    typedef typename Matrix2::index_type IndexType;
    typedef typename Matrix2::value_type ValueType;

    // The ELL portion of the HYB matrix will have 'num_entries_per_row' columns.
    // Nonzero values that do not fit within the ELL structure are placed in the 
    // COO format portion of the HYB matrix.
    
    // compute number of nonzeros in the ELL and COO portions
    size_t num_ell_entries = 0;
    for(size_t i = 0; i < src.num_rows; i++)
        num_ell_entries += thrust::min<size_t>(num_entries_per_row, src.row_offsets[i+1] - src.row_offsets[i]); 

    IndexType num_coo_entries = src.num_entries - num_ell_entries;

    dst.resize(src.num_rows, src.num_cols, 
               num_ell_entries, num_coo_entries, 
               num_entries_per_row, alignment);

    const IndexType invalid_index = cusp::ell_matrix<IndexType, ValueType, cusp::host_memory>::invalid_index;

    // pad out ELL format with zeros
    thrust::fill(dst.ell.column_indices.values.begin(), dst.ell.column_indices.values.end(), invalid_index);
    thrust::fill(dst.ell.values.values.begin(),         dst.ell.values.values.end(),         ValueType(0));

    for(size_t i = 0, coo_nnz = 0; i < src.num_rows; i++)
    {
        size_t n = 0;
        IndexType jj = src.row_offsets[i];

        // copy up to num_cols_per_row values of row i into the ELL
        while(jj < src.row_offsets[i+1] && n < num_entries_per_row)
        {
            dst.ell.column_indices(i,n) = src.column_indices[jj];
            dst.ell.values(i,n)         = src.values[jj];
            jj++, n++;
        }

        // copy any remaining values in row i into the COO
        while(jj < src.row_offsets[i+1])
        {
            dst.coo.row_indices[coo_nnz]    = i;
            dst.coo.column_indices[coo_nnz] = src.column_indices[jj];
            dst.coo.values[coo_nnz]         = src.values[jj];
            jj++; coo_nnz++;
        }
    }
}


template <typename Matrix1, typename Matrix2>
void csr_to_ell(const Matrix1& src, Matrix2& dst,
                const size_t num_entries_per_row, const size_t alignment = 32)
{
    typedef typename Matrix2::index_type IndexType;
    typedef typename Matrix2::value_type ValueType;

    // compute number of nonzeros
    size_t num_entries = 0;
#pragma omp parallel for reduction( +: num_entries)
    for(size_t i = 0; i < src.num_rows; i++)
        num_entries += thrust::min<size_t>(num_entries_per_row, src.row_offsets[i+1] - src.row_offsets[i]); 

    dst.resize(src.num_rows, src.num_cols, num_entries, num_entries_per_row, alignment);

    const IndexType invalid_index = cusp::ell_matrix<IndexType, ValueType, cusp::host_memory>::invalid_index;

    // pad out ELL format with zeros
    thrust::fill(dst.column_indices.values.begin(), dst.column_indices.values.end(), invalid_index);
    thrust::fill(dst.values.values.begin(),         dst.values.values.end(),         ValueType(0));

#pragma omp parallel for
    for(size_t i = 0; i < src.num_rows; i++)
    {
        size_t n = 0;
        IndexType jj = src.row_offsets[i];

        // copy up to num_cols_per_row values of row i into the ELL
        while(jj < src.row_offsets[i+1] && n < num_entries_per_row)
        {
            dst.column_indices(i,n) = src.column_indices[jj];
            dst.values(i,n)         = src.values[jj];
            jj++, n++;
        }
    }
}

    
template <typename Matrix1, typename Matrix2>
void csr_to_array2d(const Matrix1& src, Matrix2& dst)
{
    typedef typename Matrix2::index_type IndexType;
    typedef typename Matrix2::value_type ValueType;

    dst.resize(src.num_rows, src.num_cols);

    thrust::fill(dst.values.begin(), dst.values.end(), ValueType(0));

    for(size_t i = 0; i < src.num_rows; i++)
        for(IndexType jj = src.row_offsets[i]; jj < src.row_offsets[i+1]; jj++)
            dst(i, src.column_indices[jj]) += src.values[jj]; //sum duplicates
}


/////////////////////
// DIA Conversions //
/////////////////////

template <typename Matrix1, typename Matrix2>
void dia_to_cds(const Matrix1& src, Matrix2& dst)
{
    //typedef typename Matrix2::index_type IndexType;
    //typedef typename Matrix2::value_type ValueType;

    // allocate CDS structure
    dst.resize(src.num_rows, src.num_entries, src.diagonal_offsets.size(), 1, src.values.pitch);

    cusp::copy(src.diagonal_offsets, dst.diagonal_offsets);
    cusp::copy(src.values,           dst.values);
}

template <typename Matrix1, typename Matrix2>
void dia_to_csr(const Matrix1& src, Matrix2& dst)
{
    typedef typename Matrix2::index_type IndexType;
    typedef typename Matrix2::value_type ValueType;
    
    size_t num_entries = 0;
    size_t num_diagonals = src.diagonal_offsets.size();
    
    // count nonzero entries
    for(size_t i = 0; i < src.num_rows; i++)
    {
        for(size_t n = 0; n < num_diagonals; n++)
        {
            const IndexType j = i + src.diagonal_offsets[n];

            if(j >= 0 && static_cast<size_t>(j) < src.num_cols && src.values(i,n) != ValueType(0))
                num_entries++;
        }
    }

    dst.resize(src.num_rows, src.num_cols, num_entries);

    num_entries = 0;
    dst.row_offsets[0] = 0;

    // copy nonzero entries to CSR structure
    for(size_t i = 0; i < src.num_rows; i++)
    {
        for(size_t n = 0; n < num_diagonals; n++)
        {
            const IndexType j = i + src.diagonal_offsets[n];

            if(j >= 0 && static_cast<size_t>(j) < src.num_cols)
            {
                const ValueType value = src.values(i, n);

                if (value != ValueType(0))
                {
                    dst.column_indices[num_entries] = j;
                    dst.values[num_entries] = value;
                    num_entries++;
                }
            }
        }

        dst.row_offsets[i + 1] = num_entries;
    }
}

/////////////////////
// ELL Conversions //
/////////////////////

template <typename Matrix1, typename Matrix2>
void ell_to_coo(const Matrix1& src, Matrix2& dst)
{
    typedef typename Matrix2::index_type IndexType;
    typedef typename Matrix2::value_type ValueType;
    
    const IndexType invalid_index = cusp::ell_matrix<IndexType, ValueType, cusp::host_memory>::invalid_index;
    
    dst.resize(src.num_rows, src.num_cols, src.num_entries);

    size_t num_entries = 0;

    const size_t num_entries_per_row = src.column_indices.num_cols;

    for(size_t i = 0; i < src.num_rows; i++)
    {
        for(size_t n = 0; n < num_entries_per_row; n++)
        {
            const IndexType j = src.column_indices(i,n);
            const ValueType v = src.values(i,n);

            if(j != invalid_index)
            {
                dst.row_indices[num_entries]    = i;
                dst.column_indices[num_entries] = j;
                dst.values[num_entries]         = v;
                num_entries++;
            }
        }
    }
}

template <typename Matrix1, typename Matrix2>
void ell_to_csr(const Matrix1& src, Matrix2& dst)
{
    typedef typename Matrix2::index_type IndexType;
    typedef typename Matrix2::value_type ValueType;
    
    const IndexType invalid_index = cusp::ell_matrix<IndexType, ValueType, cusp::host_memory>::invalid_index;

    dst.resize(src.num_rows, src.num_cols, src.num_entries);

    size_t num_entries = 0;
    dst.row_offsets[0] = 0;

    const size_t num_entries_per_row = src.column_indices.num_cols;

    for(size_t i = 0; i < src.num_rows; i++)
    {
        for(size_t n = 0; n < num_entries_per_row; n++)
        {
            const IndexType j = src.column_indices(i,n);
            const ValueType v = src.values(i,n);

            if(j != invalid_index)
            {
                dst.column_indices[num_entries] = j;
                dst.values[num_entries]         = v;
                num_entries++;
            }
        }

        dst.row_offsets[i + 1] = num_entries;
    }
}

/////////////////////
// HYB Conversions //
/////////////////////

template <typename Matrix1, typename Matrix2>
void hyb_to_coo(const Matrix1& src, Matrix2& dst)
{
    typedef typename Matrix2::index_type IndexType;
    typedef typename Matrix2::value_type ValueType;
    
    dst.resize(src.num_rows, src.num_cols, src.num_entries);

    const IndexType invalid_index = cusp::ell_matrix<IndexType, ValueType, cusp::host_memory>::invalid_index;

    const size_t num_entries_per_row = src.ell.column_indices.num_cols;

    size_t num_entries  = 0;
    size_t coo_progress = 0;
    
    // merge each row of the ELL and COO parts into a single COO row
    for(size_t i = 0; i < src.num_rows; i++)
    {
        // append the i-th row from the ELL part
        for(size_t n = 0; n < num_entries_per_row; n++)
        {
            const IndexType j = src.ell.column_indices(i,n);
            const ValueType v = src.ell.values(i,n);

            if(j != invalid_index)
            {
                dst.row_indices[num_entries]    = i;
                dst.column_indices[num_entries] = j;
                dst.values[num_entries]         = v;
                num_entries++;
            }
        }

        // append the i-th row from the COO part
        while (coo_progress < src.coo.num_entries && static_cast<size_t>(src.coo.row_indices[coo_progress]) == i)
        {
            dst.row_indices[num_entries]    = i;
            dst.column_indices[num_entries] = src.coo.column_indices[coo_progress];
            dst.values[num_entries]         = src.coo.values[coo_progress];
            num_entries++;
            coo_progress++;
        }
    }
}

template <typename Matrix1, typename Matrix2>
void hyb_to_csr(const Matrix1& src, Matrix2& dst)
{
    typedef typename Matrix2::index_type IndexType;
    typedef typename Matrix2::value_type ValueType;
    
    dst.resize(src.num_rows, src.num_cols, src.num_entries);

    const IndexType invalid_index = cusp::ell_matrix<IndexType, ValueType, cusp::host_memory>::invalid_index;

    const size_t num_entries_per_row = src.ell.column_indices.num_cols;

    size_t num_entries = 0;
    dst.row_offsets[0] = 0;
    
    size_t coo_progress = 0;
    
    // merge each row of the ELL and COO parts into a single CSR row
    for(size_t i = 0; i < src.num_rows; i++)
    {
        // append the i-th row from the ELL part
        for(size_t n = 0; n < num_entries_per_row; n++)
        {
            const IndexType j = src.ell.column_indices(i,n);
            const ValueType v = src.ell.values(i,n);

            if(j != invalid_index)
            {
                dst.column_indices[num_entries] = j;
                dst.values[num_entries]         = v;
                num_entries++;
            }
        }

        // append the i-th row from the COO part
        while (coo_progress < src.coo.num_entries && static_cast<size_t>(src.coo.row_indices[coo_progress]) == i)
        {
            dst.column_indices[num_entries] = src.coo.column_indices[coo_progress];
            dst.values[num_entries]         = src.coo.values[coo_progress];
            num_entries++;
            coo_progress++;
        }

        dst.row_offsets[i + 1] = num_entries;
    }
}


/////////////////////////
// Array1d Conversions //
/////////////////////////
template <typename Matrix1, typename Matrix2>
void array2d_to_array1d(const Matrix1& src, Matrix2& dst)
{
  if (src.num_rows == 0 && src.num_cols == 0)
  {
    dst.resize(0);
  }
  else if (src.num_cols == 1)
  {
    dst.resize(src.num_rows);

    for (size_t i = 0; i < src.num_rows; i++)
      dst[i] = src(i,0);
  }
  else if (src.num_rows == 1)
  {
    dst.resize(src.num_cols);

    for (size_t j = 0; j < src.num_cols; j++)
      dst[j] = src(0,j);
  }
  else
  {
    throw cusp::format_conversion_exception("array2d to array1d conversion is only defined for row or column vectors");
  }
}

/////////////////////////
// Array2d Conversions //
/////////////////////////
template <typename Matrix1, typename Matrix2>
void array1d_to_array2d(const Matrix1& src, Matrix2& dst)
{
  dst.resize(src.size(),1);

  for (size_t i = 0; i < src.size(); i++)
    dst(i,0) = src[i];
}

template <typename Matrix1, typename Matrix2>
void array2d_to_coo(const Matrix1& src, Matrix2& dst)
{
  //typedef typename Matrix2::index_type IndexType;
  typedef typename Matrix2::value_type ValueType;

  // count number of nonzero entries in array
  size_t nnz = 0;
  
  for(size_t i = 0; i < src.num_rows; i++)
  {
    for(size_t j = 0; j < src.num_cols; j++)
    {
      if (src(i,j) != ValueType(0))
        nnz++;
    }
  }

  dst.resize(src.num_rows, src.num_cols, nnz);

  nnz = 0;

  for(size_t i = 0; i < src.num_rows; i++)
  {
    for(size_t j = 0; j < src.num_cols; j++)
    {
      if (src(i,j) != ValueType(0))
      {
        dst.row_indices[nnz]    = i;
        dst.column_indices[nnz] = j;
        dst.values[nnz]         = src(i,j);
        nnz++;
      }
    }
  }
}

template <typename Matrix1, typename Matrix2>
void array2d_to_csr(const Matrix1& src, Matrix2& dst)
{
  typedef typename Matrix2::index_type IndexType;
  typedef typename Matrix2::value_type ValueType;
  
  IndexType nnz = src.num_entries - thrust::count(src.values.begin(), src.values.end(), ValueType(0));

  dst.resize(src.num_rows, src.num_cols, nnz);

  IndexType num_entries = 0;

  for(size_t i = 0; i < src.num_rows; i++)
  {
    dst.row_offsets[i] = num_entries;

    for(size_t j = 0; j < src.num_cols; j++)
    {
      if (src(i,j) != ValueType(0))
      {
        dst.column_indices[num_entries] = j;
        dst.values[num_entries]         = src(i,j);
        num_entries++;
      }
    }
  }

  dst.row_offsets[src.num_rows] = num_entries;
}

} // end namespace host
} // end namespace detail
} // end namespace cusp