/*******************************************************************************
+
+  LEDA 3.5
+
+  matrix.h
+
+  This file is part of the LEDA research version (LEDA-R) that can be 
+  used free of charge in academic research and teaching. Any commercial
+  use of this software requires a license which is distributed by the
+  LEDA Software GmbH, Postfach 151101, 66041 Saarbruecken, FRG
+  (fax +49 681 31104).
+
+  Copyright (c) 1991-1997  by  Max-Planck-Institut fuer Informatik
+  Im Stadtwald, 66123 Saarbruecken, Germany     
+  All rights reserved.
+ 
*******************************************************************************/

#ifndef LEDA_MATRIX_H
#define LEDA_MATRIX_H

#if !defined(LEDA_ROOT_INCL_ID)
#define LEDA_ROOT_INCL_ID 350073
#include <LEDA/REDEFINE_NAMES.h>
#endif


//------------------------------------------------------------------------------
//  matrices
//------------------------------------------------------------------------------

#include <LEDA/basic.h>
#include <LEDA/vector.h>


/*{\Manpage {matrix} {}  {Real-Valued Matrices} }*/

class matrix {

/*{\Mdefinition
An instance of the data type |matrix| is a matrix of variables of type |double|.  
}*/

  vector** v;
  int  d1;
  int  d2;

  void     flip_rows(int,int);
  void     check_dimensions(const matrix&) const; 
  double&  elem(int i, int j) const { return v[i]->v[j]; }
  double** triang(const matrix&, int&) const;
    
public:

/*{\Mcreation M }*/

 matrix(int n=0, int m=0);

/*{\Mcreate creates an instance $M$ of type $matrix$, $M$ is initialized to 
           the $n \times m$ - zero matrix. }*/

  matrix(const matrix&);
  matrix(const vector&);
  matrix(int,int,double**);

  matrix& operator=(const matrix&);

 ~matrix();

  LEDA_MEMORY(matrix)


/*{\Moperations 2 4 }*/

int     dim1()  const  {  return d1; }
/*{\Mop  returns $n$, the number of rows of $M$. }*/

int     dim2()  const  {  return d2; }
/*{\Mop  returns $m$, the number of columns of $M$. }*/


vector& row(int i) const;
/*{\Mop  returns the $i$-th row of $M$ (an $m$-vector).\\
	 \precond  $0 \le i \le n-1$. }*/

vector  col(int i) const;
/*{\Mop  returns the $i$-th column of $M$ (an $n$-vector).\\
	 \precond  $0 \le i \le m-1$. }*/

matrix  trans() const;
/*{\Mop  returns  $M^T$ ($m\times n$ - matrix). }*/


matrix  inv()   const;
/*{\Mop  returns the inverse matrix of $M$.\\
	 \precond  $M$.det() $\neq$ 0. }*/

double  det()   const;
/*{\Mop  returns the determinant of $M$.\\
	 \precond  $M$ is quadratic. }*/

/*{\Moptions nextwarning=no}*/
matrix solve(const matrix&) const;
vector solve(const vector& b) const { return vector(solve(matrix(b))); }
/*{\Mop  returns vector $x$ with $M\cdot x = b$.\\
	 \precond $M$.dim1() = $M$.dim2() = $b$.dim()
	 and $M$.det() $\neq$ 0. }*/


operator vector() const; 

vector& operator[](int i)    const { return row(i); }

double& operator()(int i, int j);
/*{\Mfunop returns $M_{i,j}$. \\
	  \precond $0\le i\le n-1$ and $0\le j\le m-1$.}*/

double  operator()(int,int) const;

int     operator==(const matrix&)    const;
int     operator!=(const matrix& x)  const { return !(*this == x); }

matrix operator+(const matrix& M1);
/*{\Mbinop Addition. \\
 	   \precond $M$.dim1() = $M1$.dim1() and $M$.dim2() = $M1$.dim2().}*/

matrix operator-(const matrix& M1);
/*{\Mbinop Subtraction. \\
	   \precond $M$.dim1() = $M1$.dim1() and
	   $M$.dim2() = $M1$.dim2().}*/

matrix operator-(); // unary

matrix& operator-=(const matrix&);
matrix& operator+=(const matrix&);

matrix operator*(const matrix& M1);
/*{\Mbinop Multiplication. \\
	   \precond $M$.dim2() = $M1$.dim1().}*/

vector operator*(const vector& vec) { return vector(*this * matrix(vec)); }
/*{\Mbinop  Multiplication with vector.\\
	    \precond $M$.dim2() = $vec$.dim(). }*/

matrix operator*(double x);
/*{\Mbinop Multiplication with double x.}*/

void print(ostream& O);
/*{\Mop  prints $M$ row by row to ostream $O$. }*/

void print() { print(cout); }
/*{\Mop  prints $M$ $cout$. }*/

void read(istream& I);
/*{\Mop  reads $|M.dim1()| \times |M.dim2()|$ numbers from input stream
         $I$ and writes them row by row into matrix $M$. }*/

void read() { read(cin); }
/*{\Mop  prints $M$ from $cin$. }*/

friend ostream& operator<<(ostream& O, const matrix& M);
/*{\Mbinopfunc  writes matrix $M$ row by row to the output stream $O$. }*/

friend istream& operator>>(istream& I, matrix& M);
/*{\Mbinopfunc  reads a matrix row by row from the input stream $I$
                and assigns it to $M$. }*/

/*{\Mimplementation
Data type $matrix$ is implemented by two-dimensional arrays of double numbers. 
Operations det, solve, and inv take time $O(n^3)$, dim1, dim2, row, and col 
take constant time, all other operations take time $O(nm)$.  
The space requirement is $O(nm)$.

Be aware that the operations on vectors and matrices incur rounding error and
hence are not completely reliable. For example, if $M$ is a matrix, $b$ is a
vector, and $x$ is computed by |x = M.solve(b)| it is not necessarly true that
the test |b == M * b| evaluates to true. The types |integer_vector| and |integer_matrix| provide exact linear algebra.

}*/

};


inline char* leda_tname(const matrix*) { return "matrix"; }


#if LEDA_ROOT_INCL_ID == 350073
#undef LEDA_ROOT_INCL_ID
#include <LEDA/UNDEFINE_NAMES.h>
#endif

#endif
