// $Id: dsmatrix.h,v 1.1.1.1 2000/11/02 08:47:16 fritzi Exp $

// This file is part of the DEAL Library
// DEAL is Copyright(1995) by
// Roland Becker, Guido Kanschat, Franz-Theo Suttmeier


#ifndef __dsmatrix_h
#define __dsmatrix_h

#include <math.h>
#ifndef __dvector_h
#include <deal/dvector.h>
#endif
#include <iostream.h>
/*
CLASS
   dSMatrixStruct
   */
class dSMatrixStruct
{
  friend class dSMatrix;
  friend class fSMatrix;
  int max_dim;
  int rows, cols;
  int vec_len, max_vec_len;
  int compressed;

public:
				     /**
				      * Print the sparsity of the matrix
				      * in a format that #gnuplot# understands
				      * and which can be used to plot the
				      * sparsity pattern in a graphical
				      * way. The format consists of pairs
				      * #i j# of nonzero elements, each
				      * representing one entry of this
				      * matrix, one per line of the output
				      * file. Indices are counted from
				      * zero on, as usual. Since sparsity
				      * patterns are printed in the same
				      * way as matrices are displayed, we
				      * print the negative of the column
				      * index, which means that the
				      * #(0,0)# element is in the top left
				      * rather than in the bottom left
				      * corner.
				      *
				      * Print the sparsity pattern in
				      * gnuplot by setting the data style
				      * to dots or points and use the
				      * #plot# command.
				      */
  void print_gnuplot (ostream &out) const;
  int* rowstart;
  int* colnums;
  //////////
  void reinit(int m, int n, int max_per_row);
  //////////
  dSMatrixStruct(int m, int n, int max_per_row) 
    : max_dim(0), max_vec_len(0), rowstart(0), colnums(0)
      {
	reinit(m,n,max_per_row);
      }
  //////////
  dSMatrixStruct(int n, int max_per_row)
    : max_dim(0), max_vec_len(0), rowstart(0), colnums(0)
      {
	reinit(n,n,max_per_row);
      }
  //////////
  ~dSMatrixStruct()
    {
      delete[] rowstart;
      delete[] colnums;
    }
  //////////
  void compress();
  //////////
  int operator () (int i, int j);
  //////////
  void add(int i, int j);
  //////////
  void add_matrix(int n, int* rowcols);
};

/*
CLASS
   dSMatrix
   */
class dSMatrix
{
 protected:
  dSMatrixStruct& cols;
  double* val;
  int max_len;
 public:
  int* rowstart() { return cols.rowstart; }
  int* colnums()  { return cols.colnums; }
 public:
  //////////
  void reinit();
  //////////
  dSMatrix(dSMatrixStruct& c)
      : max_len(0), cols(c), val(0)
    {
      reinit();
    }
  //////////
  ~dSMatrix()
    {
      delete[] val;
    }

  //////////
  int m() const { return cols.rows; }
  //////////
  int n() const { return cols.cols; }

  double get(int i,int j) { return  val[cols(i,j)]; }
  //////////
  void set(int i,int j,double value) { val[cols(i,j)] = value; }
  //////////
  void add(int i,int j,double value) { val[cols(i,j)]+= value; }
  
  //////////
  void vmult (dVector& dst,const dVector& src,int adding=0);
  //////////
  void Tvmult(dVector& dst,const dVector& src,int adding=0);
  
  //////////
  double residual (dVector& dst,const dVector& x,const dVector& b);
  //////////
  void SSOR_precond(dVector& dst,const dVector& src,double om = 1.);
  //////////
  void SOR_precond(dVector& dst,const dVector& src,double om = 1.);
  //////////
  void SSOR(dVector& dst,double om = 1.);
  //////////
  void SOR(dVector& dst,double om = 1.);
  //////////
  void precondition(dVector& dst,const dVector& src) { dst=src; }
};
#endif