// nUtil - An utility Library for gnurbs
// Copyright (C) 2008-2026 Eric Bechet
//
// See the LICENSE file for contributions and license information.
// Please report all bugs and problems to <bechet@cadxfem.org>.
//

//---------------------------------------------------------------------------
#ifndef SPARSE_MATRIX_H
#define SPARSE_MATRIX_H
//---------------------------------------------------------------------------

#include "sparse_iterator.h"
#include <ostream>
#include <iostream>
#include <vector>
#include <map>
#include <set>

class sparse_matrix
{
public :
  typedef sparse_iterator iterator;

  sparse_matrix(int i, int j) :endval(0.0)
  {
    alloc(i,j);
  }

  sparse_matrix() :endval(0.0) {}

  sparse_matrix(sparse_matrix &m) :endval(* (m.endval))
  {
    alloc(m.lines.size(),m.columns.size());
    for (unsigned i=0;i<m.lines.size();++i)
      for (iterator it=m.beginl(i); it!=m.endl(i);it=it.nextl())
        setval(it.line(),it.column(),*it);
  }


  ~sparse_matrix()
  {
    clean();
  }

  iterator& beginl(int i)
  {
    return * ((* (lines[i].begin())).second);
  }
  iterator& beginc(int j)
  {
    return * ((* (columns[j].begin())).second);
  }
  iterator& endl(int i)
  {
    return endval;
  }
  iterator& endc(int j)
  {
    return endval;
  }

  int nbl()
  {
    return lines.size();
  }
  int nbc()
  {
    return columns.size();
  }

  void max_abs(int&imax,int&jmax,double& max,std::set<int> &map);

  void llinearcombination(std::map<int,double> &coeffs,int dest);

  template <class T> void diagonalize(std::vector<std::pair<int,int> > &map,
                                      std::vector<int>& indep,
                                      T linearcombination);

  void addval(int i,int j,double val);
  void setval(int i,int j,double val);
  double delval(int i,int j);


  iterator& operator()(int i,int j)
  {
    sparse_container::iterator it=lines[i].find(j);
    if (it!=lines[i].end())
      return * ((*it).second) ;
    else
      return endval;
  }


  void alloc(int i, int j);
  void clean();
  void display(std::ostream &out);
  void mult(sparse_matrix &a, sparse_matrix &b);
  unsigned nnzero(void);
  void displaycompressed(std::ostream &out);

private :
  std::vector<sparse_container> lines;
  std::vector<sparse_container> columns;
  sparse_iterator endval;
};







template <class T> void sparse_matrix::diagonalize(std::vector<std::pair<int,int> >& mappingv,
    std::vector<int>& indepv,
    T linearcombination)
{
  int imax,jmax;
  std::set<int> maps;
  double max;
  max_abs(imax,jmax,max,maps);

  while (max!=0.0)
  {
    mappingv.push_back(std::pair<int,int> (imax,jmax));
    maps.insert(imax);
    //        std::cout << "pivot i=" << imax << " j=" << jmax << " val=" << max << std::endl;
    iterator it=beginc(jmax);
    while (it!=endc(jmax))
    {
      iterator it2=it.nextc();
      if (maps.find(it.line()) ==maps.end())
      {
        double val=*it;
        int line=it.line();
        double ratio=val/max;
        std::map<int,double> group;
        group[line]=1;
        group[imax]=-ratio;
        //  std::cout << "LIN " << line << " " << imax << " 1.0 " << -ratio << std::endl;
        llinearcombination(group,line);
        linearcombination(group,line);
        delval(line,jmax);
      }
      it=it2;
    }
    std::map<int,double> group;
    group[imax]=1/max;
    llinearcombination(group,imax);
    linearcombination(group,imax);
    //    display(std::cout);
    max_abs(imax,jmax,max,maps);
  }
  //    std::cout << "remonte" << std::endl;
  maps.clear();
  std::vector<std::pair<int,int> >::reverse_iterator itr;
  std::vector<std::pair<int,int> >::reverse_iterator itre=mappingv.rend();
  for (itr=mappingv.rbegin();itr!=itre;++itr)
  {
    imax=itr->first;
    jmax=itr->second;
    max=* (*this)(imax,jmax);
    maps.insert(imax);
    //        std::cout << "pivot i=" << imax << " j=" << jmax << " val=" << max << std::endl;
    iterator it=beginc(jmax);
    while (it!=endc(jmax))
    {
      iterator it2=it.nextc();
      if (maps.find(it.line()) ==maps.end())
      {
        double val=*it;
        int line=it.line();
        double ratio=val/max;
        std::map<int,double> group;
        group[line]=1;
        group[imax]=-ratio;
        // std::cout << "LIN " << line << " " << imax << " 1.0 " << -ratio << std::endl;
        llinearcombination(group,line);
        linearcombination(group,line);
        delval(line,jmax);
      }
      std::map<int,double> group;
      group[imax]=1/max;
      llinearcombination(group,imax);
      linearcombination(group,imax);
      it=it2;
    }
    // display(std::cout);
  }
}




#endif //SPARSE_MATRIX_H