tnt_sparse_matrix.h

#ifndef TNT_SPARSE_MATRIX_H
#define TNT_SPARSE_MATRIX_H
 
/*
*
* Template Numerical Toolkit (TNT)
*
* Mathematical and Computational Sciences Division
* National Institute of Technology,
* Gaithersburg, MD USA
*
*
* This software was developed at the National Institute of Standards and
* Technology (NIST) by employees of the Federal Government in the course
* of their official duties. Pursuant to title 17 Section 105 of the
* United States Code, this software is not subject to copyright protection
* and is in the public domain. NIST assumes no responsibility whatsoever for
* its use by other parties, and makes no guarantees, expressed or implied,
* about its quality, reliability, or any other characteristic.
*
*/
 
 
 
#include <vector>
#include "tnt_vector.h"
#include "tnt_sparse_vector.h"
 
namespace TNT
{
 
//namespace Linear_Algebra
//{
 
 
#if 0
template <class T, class Integer>
class Sparse_Matrix_Coordinate_Element
{
        private:
 
                T val_;
                Integer row_index_;
                Integer col_index_;
 
        public:
 
                Sparse_Matrix_Coordinate_Element(
                        const T& a, const Integer &i, const Integer &j) : 
                        val_(a), row_index_(i), col_index(i) {}
                
 
                const T& value() const { return val_; } 
                Integer row_index() const { return row_index_; } 
                Integer col_index() const { return col_index_; } 
 
                T& value() { return val_; }
                Integer& row_index() { return row_index_; }
                Integer& col_index() { return col_index_; }
 
};
#endif
 
template <class T>
class Sparse_Matrix
{
  private:
 
        // compressed row storage M rows of <value, col_index> pairs
        //
        std::vector< Sparse_Vector<T> > S_;    
 
    int num_rows_;        // number of rows
    int num_cols_;        // number of cols
        int num_nonzeros_;      // number of nonzeros
 
                                                    // Used only in multi-step constructions.  This
                                                    // allows one to build a sparse matrix with 
                                                    // multiple calls to insert().
                                                    //
        int internal_state_;    // 0 if closed (no more inserts) , 1 if open;
 
 
 
public:
 
 
    Sparse_Matrix(Subscript M, Subscript N):
                                S_(M),  
                                num_rows_(M), num_cols_(N), num_nonzeros_(0),
                                internal_state_(1) {};
 
            
 
    Sparse_Matrix(Subscript M, Subscript N, Subscript nz, const T* val,
                const Subscript *r, const Subscript *c):
                    S_(M),
                    num_rows_(M), 
                    num_cols_(N), 
                    num_nonzeros_(0),
                    internal_state_(1)
                    {
 
                        insert(nz, val, r, c);  
                        close();
                    };
   
 
    int is_closed() { return internal_state_; }
 
    void   insert(const T& val, Subscript i, Subscript j)
    {
          if (internal_state_ == 0) return;
 
 
            S_[i].insert(val, j);
            num_nonzeros_++;
    }
 
    void   insert(Subscript nz, const T*  val, const Subscript *i, 
                    const Subscript *j)
    {
          if (internal_state_ == 0) return;
 
            for (int count=0; count<nz; count++)
            {
                insert(val[count], i[count], j[count]); 
            }
    }
 
    void   insert_base_one(const T& val, Subscript i, Subscript j)
    {
            insert_one_base(val, i-1, j-1);
    }
 
    void   insert_base_one(Subscript nz, const T*  val, const Subscript *i, 
                    const Subscript *j)
    {
            for (int count=0; count<nz; count++)
            {
                insert(val[count], i[count]-1, j[count]-1); 
            }
    }
 
    
    void        close()
    {
            /* 
                    After this, there are no more inserts. Now one could optimize 
               storage layout.
                 
            */
 
            internal_state_ = 0;
    }
 
  inline   int    num_rows() const {return num_rows_;}
  inline   int    num_cols() const {return num_cols_;}
  inline   int    num_columns() const {return num_cols_;}
  int          num_nonzeros() const {return num_nonzeros_;}
 
 
    Vector<T> diag() const
    {
        int minMN = num_rows() < num_columns() ? num_rows() : num_columns();  
        Vector<T> diag_(minMN, T(0));
 
        for (int i=0; i<minMN; i++)
        {
            for (typename Sparse_Vector<T>::const_iterator p = S_[i].begin();
                        p < S_[i].end(); p++ )
            {
                    if (p->index() == i)
                    {
                        diag_[i] += p->value();
                        break;
                    }
            }
        }
        return diag_;
    }
 
    Vector<T> mult(const Vector<T> &x) const
    {
        int M = num_rows();
        Vector<T> y(M);
        for (int i=0; i<M; i++)
        {
            y[i] = dot_product(S_[i], x);
        }
        return y;
    }
 
 
    inline double norm() const
    {
        T sum(0.0);
        for (int i=0; i<num_rows_; i++)
        {
            for (typename Sparse_Vector<T>::const_iterator p = S_[i].begin();
                        p < S_[i].end(); p++ )
            {
                sum += p->value() * p->value();
            }
        }
        return sqrt(sum);
    }
 
    std::ostream & print(std::ostream &s) const
    {
        for (int i=0; i<num_rows_; i++)
        {
            for (typename Sparse_Vector<T>::const_iterator p = S_[i].begin();
                        p < S_[i].end(); p++ )
            {
                 s << "( " << p->value() << " , " << i << ", " << p->index() << " )\n";
            }
        }
 
    return s;
    }
 
    std::ostream & print_base_one(std::ostream &s) const
    {
        for (int i=0; i<num_rows_; i++)
        {
            for (typename Sparse_Vector<T>::const_iterator p = S_[i].begin();
                        p < S_[i].end(); p++ )
            {
                 s <<"( "<<p->value()<<" , "<<i+1<<", "<< p->index()+1 << " )\n";
            }
        }
 
    return s;
    }
            
            
};
 
 
template <class T>
inline Vector<T> operator*(const Sparse_Matrix<T> &S, const Vector<T> &x) 
{
    return S.mult(x);
}
 
template <class T>
inline double norm(const Sparse_Matrix<T> &S)
{
    return S.norm();
}
 
template <class T>
inline  std::ostream& operator<<(std::ostream &s, const Sparse_Matrix<T> &A)
{
    return A.print(s);
}
 
 
 
//}  /* namspace TNT::Linear_Algebra */
}    /* namspace TNT */
 
#endif