lilc_matrix_declarations.h

// -*- mode: c++ -*-
#ifndef _LILC_MATRIX_DECLARATIONS_H_
#define _LILC_MATRIX_DECLARATIONS_H_

#include <algorithm>
#include <cmath>
#include <cassert>
#include <iostream>
#include <iterator>
#include <set>

#include "swap_struct.h"

template<class el_type> 
class lilc_matrix : public lil_sparse_matrix<el_type>
{
public:
    
    //-------------- typedefs and inherited variables --------------//
    using lil_sparse_matrix<el_type>::m_idx;
    using lil_sparse_matrix<el_type>::m_x;
    using lil_sparse_matrix<el_type>::m_n_rows;
    using lil_sparse_matrix<el_type>::m_n_cols;
    using lil_sparse_matrix<el_type>::n_rows;
    using lil_sparse_matrix<el_type>::n_cols;
    using lil_sparse_matrix<el_type>::nnz;
    using lil_sparse_matrix<el_type>::nnz_count;
    using lil_sparse_matrix<el_type>::eps;


    typedef typename lil_sparse_matrix<el_type>::idx_vector_type idx_vector_type;
    typedef typename lil_sparse_matrix<el_type>::elt_vector_type elt_vector_type;
    
    typedef typename idx_vector_type::iterator idx_it;
    typedef typename elt_vector_type::iterator elt_it;

    std::vector< std::vector< int > > list; 
    
    std::vector<int> row_first; 
    std::vector<int> col_first; 
    
    block_diag_matrix<el_type> S; 
    
    //-------------- types of pivoting procedures ----------------//
    struct pivot_type {
        enum {
            BKP, 
            ROOK
        };
    };
    
public:
    
    lilc_matrix (int n_rows = 0, int n_cols = 0): 
    lil_sparse_matrix<el_type> (n_rows, n_cols) 
    {
        m_x.reserve(n_cols);
        m_idx.reserve(n_cols);
    }
    
    //----Matrix referencing/filling----//
    
    inline virtual el_type coeff(const int& i, const int& j, int offset = 0) const 
    {   
        //invariant: first elem in each col of a is the diagonal elem if it exists.
        if (i == j) {
            if (m_idx[j].size() == 0) return 0;
            return (m_idx[j][0] == i ? m_x[j][0] : 0);
        }
        
        for (int k = offset, end = m_idx[j].size(); k < end; k++) {
            if (m_idx[j][k] == i) return m_x[j][k];
        }
        
        return 0;
    }
    
    inline bool coeffRef(const int& i, const int& j, std::pair<idx_it, elt_it>& its)
    {   
        for (unsigned int k = 0; k < m_idx[j].size(); k++) {
            if (m_idx[j][k] == i) {
                its = make_pair(m_idx[j].begin() + k, m_x[j].begin() + k);
                return true;
            }
        }
        
        its = make_pair(m_idx[j].end(), m_x[j].end());
        return false;
    }
    
    void resize(int n_rows, int n_cols)
    {
        m_n_rows = n_rows;
        m_n_cols = n_cols;
    
        m_x.clear();
        m_idx.clear();
        row_first.clear();
        col_first.clear();
        list.clear();
    
        m_x.resize(n_cols);
        m_idx.resize(n_cols);
        
        row_first.resize(n_cols, 1);
        col_first.resize(n_cols, 1);
        list.resize(n_cols);
        
        S.resize(n_cols, 1);
        for (int i = 0; i < n_cols; i++) {
            m_x[i].clear();
            m_idx[i].clear();
            list[i].clear();
        }
    }
    
    //-----Reorderings/Rescalings------//
    void find_root(int& s);
    
    inline bool find_level_set(vector<int>& lvl_set, vector<bool>& visited);
    
    void sym_rcm(vector<int>& perm);
    
    inline void sym_amd(vector<int>& perm);
        
    void sym_perm(vector<int>& perm);
    
    void sym_equil();
    
    //----Factorizations----//
    void ildl(lilc_matrix<el_type>& L, block_diag_matrix<el_type>& D, idx_vector_type& perm, const double& fill_factor, const double& tol, const double& pp_tol, int piv_type = pivot_type::BKP);
    
    void ildl_inplace(block_diag_matrix<el_type>& D, idx_vector_type& perm, const double& fill_factor, const double& tol, const double& pp_tol, int piv_type = pivot_type::BKP);
    
    //------Helpers------//
    void backsolve(const elt_vector_type& b, elt_vector_type& x) {
        assert(b.size() == x.size());
        x = b;
        // simple forward substitution
        for (int i = 0; i < m_n_cols; i++) {
            x[i] /= m_x[i][0];
            for (int k = 1; k < m_idx[i].size(); k++) {
                x[m_idx[i][k]] -= x[i]*m_x[i][k];
            }
        }
    }
    
    void forwardsolve(const elt_vector_type& b, elt_vector_type& x) {
        assert(b.size() == x.size());
        // simple back substitution
        for (int i = m_n_cols-1; i >= 0; i--) {
            x[i] = b[i]/m_x[i][0];
            for (int k = 1; k < m_idx[i].size(); k++) {
                x[i] -= x[m_idx[i][k]]*m_x[i][k]/m_x[i][0];
            }
        }
    }
    
    void multiply(const elt_vector_type& x, elt_vector_type& y, bool full_mult = true) {
        y.clear(); y.resize(x.size(), 0);
        for (int i = 0; i < m_n_cols; i++) {
            for (int k = 0; k < m_idx[i].size(); k++) {
                y[m_idx[i][k]] += x[i]*m_x[i][k];
                if (full_mult && i != m_idx[i][k]) {
                    y[i] += x[m_idx[i][k]]*m_x[i][k];
                }
            }
        }
    }
    
    inline void pivot(swap_struct<el_type>& s, vector<bool>& in_set, lilc_matrix<el_type>& L, const int& k, const int& r);
    
    inline void pivotA(swap_struct<el_type>& s, vector<bool>& in_set, const int& k, const int& r);
    
    template <class Container>
    inline void ensure_invariant(const int& j, const int& k, Container& con, bool update_list = false) {
        int offset;
        if (update_list) offset = row_first[j];
        else offset = col_first[j];
        
        if ((offset >= (int) con.size()) || con.empty() || con[offset] == k) return;
        
        int i, min(offset);
        for (i = offset; i < (int) con.size(); i++) {
            if (con[i] == k) {
                min = i; 
                break;
            } else if ( con[i] < con[min] ) {
                min = i;
            }
        }
        
        if (update_list)
            std::swap(con[offset], con[min]);
        else {
            std::swap(con[offset], con[min]);
            std::swap(m_x[j][offset], m_x[j][min]);
        }
    }
    
    inline void advance_first(const int& k) {
        for (idx_it it = list[k].begin(); it != list[k].end(); it++) {  
            ensure_invariant(*it, k, m_idx[*it]); //make sure next element is good before we increment.
            col_first[*it]++; //should have ensured invariant now
        }
    }
    
    inline void advance_list(const int& k) {
        for (idx_it it = m_idx[k].begin(); it != m_idx[k].end(); it++) {
            if (*it == k) continue;
            ensure_invariant(*it, k, list[*it], true); //make sure next element is good.
            row_first[*it]++; //invariant ensured.
        }           
    }
    
    //----IO Functions----//
    
    std::string to_string () const;
    
    bool load(std::string filename);

    bool load(const std::vector<int>& ptr, const std::vector<int>& row, const std::vector<el_type>& val);

    bool load(const int* ptr, const int* row, const el_type* val, int dim);
    
    bool save(std::string filename, bool sym = false);

};

//------------------ include files for class functions -------------------//

#include "lilc_matrix_find_level_set.h"
#include "lilc_matrix_find_root.h"
#include "lilc_matrix_sym_rcm.h"
#include "lilc_matrix_sym_amd.h"
#include "lilc_matrix_sym_perm.h"
#include "lilc_matrix_sym_equil.h"
#include "lilc_matrix_ildl_helpers.h"
#include "lilc_matrix_ildl.h"
#include "lilc_matrix_ildl_inplace.h"
#include "lilc_matrix_pivot.h"
#include "lilc_matrix_load.h"
#include "lilc_matrix_save.h"
#include "lilc_matrix_to_string.h"

#endif