SparseMatrixOperator.h

Go to the documentation of this file.

/*
 * 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.
 * 
 * Written (W) 2013 Soumyajit De
 */

#ifndef SPARSE_MATRIX_OPERATOR_H_
#define SPARSE_MATRIX_OPERATOR_H_

#include <shogun/lib/config.h>
#include <shogun/lib/SGSparseVector.h>
#include <shogun/lib/SGSparseMatrix.h>
#include <shogun/mathematics/linalg/linop/MatrixOperator.h>

namespace shogun
{
template<class T> class SGVector;
template<class T> class SGSparseMatrix;

struct SparsityStructure
{
    SparsityStructure() : m_num_rows(0), m_ptr(NULL) {}

    SparsityStructure(index_t* row_offsets, index_t* column_indices,
        index_t num_rows)
    : m_num_rows(num_rows),
        m_ptr(new int32_t*[num_rows]())
    {
        for (index_t i=0; i<m_num_rows; ++i)
        {
            index_t current_index=row_offsets[i];
            index_t new_index=row_offsets[i+1];
            index_t length_row=(new_index-current_index);

            m_ptr[i]=new int32_t[length_row+1]();
            m_ptr[i][0]=length_row;
    
            for (index_t j=1; j<=length_row; ++j)
                m_ptr[i][j]=column_indices[current_index++];
        }
    }

    ~SparsityStructure()
    {
        for (index_t i=0; i<m_num_rows; ++i)
            delete [] m_ptr[i];
        delete [] m_ptr;
    }

    void display_sparsity_structure()
    {
        for (index_t i=0; i<m_num_rows; ++i)
        {
            index_t nnzs=m_ptr[i][0];
            SG_SPRINT("Row number %d. Number of Non-zeros %d. Colums ", i, nnzs);
            for(index_t j=1; j<=nnzs; ++j)
            {
                SG_SPRINT("%d", m_ptr[i][j]);
                if (j<nnzs)
                    SG_SPRINT(", ");
            }
            SG_SPRINT("\n");
        }
    }

    index_t m_num_rows;

    int32_t **m_ptr;
};


template<class T> class CSparseMatrixOperator : public CMatrixOperator<T>
{
typedef bool supports_complex128_t;

public:
    CSparseMatrixOperator();

    explicit CSparseMatrixOperator(SGSparseMatrix<T> op);

    CSparseMatrixOperator(const CSparseMatrixOperator<T>& orig);

    ~CSparseMatrixOperator();

    virtual SGVector<T> apply(SGVector<T> b) const;

    virtual void set_diagonal(SGVector<T> diag);

    virtual SGVector<T> get_diagonal() const;

    SGSparseMatrix<T> get_matrix_operator() const;

    SparsityStructure* get_sparsity_structure(int64_t power=1) const;

    template<class Scalar>
    inline operator CSparseMatrixOperator<Scalar>*() const
    {
        REQUIRE(m_operator.sparse_matrix, "Matrix is not initialized!\n");
        typedef SGSparseVector<Scalar> vector;
        typedef SGSparseVectorEntry<Scalar> entry;

        vector* rows=SG_MALLOC(vector, m_operator.num_vectors);

        for (index_t i=0; i<m_operator.num_vectors; ++i)
        {
            entry* features=SG_MALLOC(entry, m_operator[i].num_feat_entries);
            for (index_t j=0; j<m_operator[i].num_feat_entries; ++j)
            {
                features[j].feat_index=m_operator[i].features[j].feat_index;
                features[j].entry=static_cast<Scalar>(m_operator[i].features[j].entry);
            }
            rows[i].features=features;
            rows[i].num_feat_entries=m_operator[i].num_feat_entries;
        }

        SGSparseMatrix<Scalar> casted_m;
        casted_m.sparse_matrix=rows;
        casted_m.num_vectors=m_operator.num_vectors;
        casted_m.num_features= m_operator.num_features;

        return new CSparseMatrixOperator<Scalar>(casted_m);
    }

    virtual const char* get_name() const
    {
        return "SparseMatrixOperator";
    }

private:
    SGSparseMatrix<T> m_operator;

    void init();

};

}
#endif // SPARSE_MATRIX_OPERATOR_H_