#ifndef MATRIX_H
#define MATRIX_H

/*
 * METRIC --- Mode expansion modeling in integrated optics / photonics
 * http://metric.computational-photonics.eu/ 
 */

/* 
 * matrix.h
 * Basic vectors, matrices, third order arrays, some linear algebra routines 
 */

#include <stdlib.h>
#include <stdio.h>
#include "complex.h"

// #define SAFE  // remove comment before compilation to enable range checks 
                 // for (.)-access to all arrays

#define INC_EIGEN  // introduce comment before compilation to  
                   // exclude the EIGEN subroutines

class Bvector;
class Ivector;
class Dvector;
class Cvector;
class Bmatrix;
class Imatrix;
class Dmatrix;
class Cmatrix;
class Barray3;
class Iarray3;
class Darray3;
class Carray3;

/* -------------------------------------------------------------------- */
/* - vectors ---------------------------------------------------------- */
/* -------------------------------------------------------------------- */

/* Bvector: bool vector */
class Bvector 
{
public:
	// number of elements
	int nel;
	// initialize
	Bvector();
	Bvector(int ne);
	// destroy
	~Bvector();
	// copy constructor 
	Bvector(const Bvector& s);
	// assignment 
	Bvector& operator=(const Bvector& s);

	// free allocated memory
	void strip();

	// input
	void read(FILE *dat);
	// output
	void write(FILE *dat) const;

	// output to MATLAB script file name.m
	void mfile(const char *name) const;

	// subscripting
#ifdef SAFE
	bool& operator() (const int& i) 
	{ 
		if(nel>i && i>=0) return m[i]; 
		else 
		{
			fprintf(stderr, "\n-- Bvector(%d): set(%d)\n", nel, i);
			exit(1);
		}
	}
	bool operator() (const int& i) const 
	{ 
		if(nel>i && i>=0) return m[i]; 
		else 
		{
			fprintf(stderr, "\n-- Bvector(%d): get(%d)\n", nel, i);
			exit(1);
		}
	}
#else
	inline bool& operator() (const int& i) { return m[i]; }
	inline bool operator() (const int& i) const { return m[i]; }
#endif

	// set all elements to d
	void init(bool b);

	// extract a subvector: n entries starting at position x
	Bvector subvector(int n, int x) const;
	// set the entries beginning at position x to the values of
	// the supplied smaller vector 
	void setsubvector(const Bvector &s, int x);

	// append an entry to the end of the vector 
	void append(bool e);
	// insert value e at position p,
	// shift the following elements accordingly 
	void insert(bool e, int p);
	// remove the element at position p, 
	// shift the following elements accordingly
	void remove(int p);

private:
	bool *m;
};

/* Ivector: int vector */
class Ivector 
{
public:
	// number of elements
	int nel;
	// initialize
	Ivector();
	Ivector(int ne);
	// destroy
	~Ivector();
	// copy constructor 
	Ivector(const Ivector& s);
	// assignment 
	Ivector& operator=(const Ivector& s);

	// free allocated memory
	void strip();

	// input
	void read(FILE *dat);
	// output
	void write(FILE *dat) const;

	// output to MATLAB script file name.m
	void mfile(const char *name) const;

	// subscripting
#ifdef SAFE
	int& operator() (const int& i) 
	{ 
		if(nel>i && i>=0) return m[i]; 
		else 
		{
			fprintf(stderr, "\n-- Ivector(%d): set(%d)\n", nel, i);
			exit(1);
		}
	}
	int operator() (const int& i) const 
	{ 
		if(nel>i && i>=0) return m[i]; 
		else 
		{
			fprintf(stderr, "\n-- Ivector(%d): get(%d)\n", nel, i);
			exit(1);
		}
	}
#else
	inline int& operator() (const int& i) { return m[i]; }
	inline int operator() (const int& i) const { return m[i]; }
#endif

	// set all elements to d
	void init(int d);

	// extract a subvector: n entries starting at position x
	Ivector subvector(int n, int x) const;
	// set the entries beginning at position x to the values of
	// the supplied smaller vector 
	void setsubvector(const Ivector &s, int x);

	// append an entry to the end of the vector 
	void append(int e);
	// insert value e at position p,
	// shift the following elements accordingly 
	void insert(int e, int p);
	// remove the element at position p, 
	// shift the following elements accordingly
	void remove(int p);

	// the Euclidian norm 
	double norm() const;
	// ... squared 
	double sqnorm() const;

	// sort the vector elements: d>=0: ascending, else descending
	void sort(int d);
	// ascending
	void sort();

	// minimum and maximum of the vector elements, corresponding indices
	int min() const;
	int min(int& mi) const;
	int max() const;
	int max(int& mi) const;

private:
	int *m;
};

/* Dvector: double vector */
class Dvector 
{
	friend class Cvector; 
	friend class Dmatrix; 
	friend class Cmatrix; 
public:
	// number of elements
	int nel;
	// initialize
	Dvector();
	Dvector(int ne);
	// destroy
	~Dvector();
	// copy constructor 
	Dvector(const Dvector& s);
	// assignment 
	Dvector& operator=(const Dvector& s);

	// free allocated memory
	void strip();

	// input
	void read(FILE *dat);
	// output
	void write(FILE *dat) const;

	// output to MATLAB script file name.m
	void mfile(const char *name) const;

	// subscripting
#ifdef SAFE
	double& operator() (const int& i) 
	{
		if(nel>i && i>=0) return m[i]; 
		else 
		{
			fprintf(stderr, "\n-- Dvector(%d): set(%d)\n", nel, i);
			exit(1);
		}
	}
	double operator() (const int& i) const 
	{ 
		if(nel>i && i>=0) return m[i]; 
		else 
		{
			fprintf(stderr, "\n-- Dvector(%d): get(%d)\n", nel, i);
			exit(1);
		}
	}
#else
	inline double& operator() (const int& i) { return m[i]; }
	inline double operator() (const int& i) const { return m[i]; }
#endif

	// set all elements to d
	void init(double d);

	// extract a subvector: n entries starting at position x
	Dvector subvector(int n, int x) const;
	// set the entries beginning at position x to the values of
	// the supplied smaller vector 
	void setsubvector(const Dvector &s, int x);

	// append an entry to the end of the vector 
	void append(double e);
	// insert value e at position p,
	// shift the following elements accordingly 
	void insert(double e, int p);
	// remove the element at position p, 
	// shift the following elements accordingly
	void remove(int p);

	// the Euclidian norm 
	double norm() const;
	// ... squared 
	double sqnorm() const;

	// sort the vector elements: d>=0: ascending, else descending
	void sort(int d);
	// ascending
	void sort();

	// minimum and maximum of the vector elements, corresponding indices
	double min() const;
	double min(int& mi) const;
	double max() const;
	double max(int& mi) const;

	// addition & assignment
	void addeq(const Dvector& b);
	// subtraction & assignment
	void subeq(const Dvector& b);
	// multiplication & assignment
	void multeq(double c);

	/* ----------------------------------------------------- */

	/* multiplications */
	friend Cmatrix mult(const Cmatrix& a, const Cmatrix& b);
	friend Cmatrix mult(const Cmatrix& a, const Dmatrix& b);
	friend Cmatrix mult(const Dmatrix& a, const Cmatrix& b);
	friend Dmatrix mult(const Dmatrix& a, const Dmatrix& b);
	friend Cvector mult(const Cmatrix& a, const Cvector& b);
	friend Cvector mult(const Dmatrix& a, const Cvector& b);
	friend Cvector mult(const Cmatrix& a, const Dvector& b);
	friend Dvector mult(const Dmatrix& a, const Dvector& b);
	friend Cmatrix mult(const Cmatrix& a, const Complex& b);
	friend Cmatrix mult(const Cmatrix& a, double b);
	friend Cmatrix mult(const Dmatrix& a, const Complex& b);
	friend Dmatrix mult(const Dmatrix& a, double b);
	friend Cmatrix mult(const Complex& a, const Cmatrix& b);
	friend Cmatrix mult(const Complex& a, const Dmatrix& b);
	friend Cmatrix mult(double a,  const Cmatrix& b);
	friend Dmatrix mult(double a,  const Dmatrix& b);
	friend Cvector mult(const Cvector& a, const Complex& b);
	friend Cvector mult(const Cvector& a, double b);
	friend Cvector mult(const Dvector& a, const Complex& b);
	friend Dvector mult(const Dvector& a, double b);
	friend Cvector mult(const Complex& a, const Cvector& b);
	friend Cvector mult(const Complex& a, const Dvector& b);
	friend Cvector mult(double a,  const Cvector& b);
	friend Dvector mult(double a,  const Dvector& b);

	/* additions */
	friend Cmatrix add(const Cmatrix& a, const Cmatrix& b);
	friend Cmatrix add(const Cmatrix& a, const Dmatrix& b);
	friend Cmatrix add(const Dmatrix& a, const Cmatrix& b);
	friend Dmatrix add(const Dmatrix& a, const Dmatrix& b);
	friend Cvector add(const Cvector& a, const Cvector& b);
	friend Cvector add(const Dvector& a, const Cvector& b);
	friend Cvector add(const Cvector& a, const Dvector& b);
	friend Dvector add(const Dvector& a, const Dvector& b);

	/* subtraction */
	friend Cmatrix sub(const Cmatrix& a, const Cmatrix& b);
	friend Cmatrix sub(const Cmatrix& a, const Dmatrix& b);
	friend Cmatrix sub(const Dmatrix& a, const Cmatrix& b);
	friend Dmatrix sub(const Dmatrix& a, const Dmatrix& b);
	friend Cvector sub(const Cvector& a, const Cvector& b);
	friend Cvector sub(const Dvector& a, const Cvector& b);
	friend Cvector sub(const Cvector& a, const Dvector& b);
	friend Dvector sub(const Dvector& a, const Dvector& b);

	/* multiplications, 
	   the argument d represents the elements of a square 
	   diagonal matrix */
	friend Dmatrix diagmult(const Dvector& d, const Dmatrix& m);
	friend Dmatrix diagmult(const Dmatrix& m, const Dvector& d);
	friend Dvector diagmult(const Dvector& d, const Dvector& a);
	friend Cmatrix diagmult(const Cvector& d, const Cmatrix& m);
	friend Cmatrix diagmult(const Cmatrix& m, const Cvector& d);
	friend Cvector diagmult(const Cvector& d, const Cvector& a);

	/* scalar products */
	friend Complex scalp(const Cvector& a, const Cvector& b);
	friend Complex scalp(const Cvector& a, const Dvector& b);
	friend Complex scalp(const Dvector& a, const Cvector& b);
	friend double  scalp(const Dvector& a, const Dvector& b);

	/* ----------------------------------------------------- */
	
	/* sort a, b in parallel, such that a is ascending */
	friend void sort(Dvector& a, Dvector &b);

private:
	double *m;
};

/* Cvector: complex vector */
class Cvector 
{
	friend class Dvector; 
	friend class Dmatrix; 
	friend class Cmatrix; 
public:
	// number of elements
	int nel;
	// initialize
	Cvector();
	Cvector(int ne);
	Cvector(const Dvector& rp);
	Cvector(const Dvector& rp, const Dvector& ip);
	// destroy
	~Cvector();
	// copy constructor 
	Cvector(const Cvector& s);
	// assignment 
	Cvector& operator=(const Cvector& s);

	// free allocated memory
	void strip();

	// input
	void read(FILE *dat);
	// output
	void write(FILE *dat) const;

	// output to MATLAB script file name.m
	void mfile(const char *name) const;

	// subscripting
#ifdef SAFE
	Complex& operator() (const int& i) 
	{ 
		if(nel>i && i>=0) return m[i]; 
		else 
		{
			fprintf(stderr, "\n-- Cvector(%d): set(%d)\n", nel, i);
			exit(1);
		}
	}
	Complex operator() (const int& i) const 
	{ 
		if(nel>i && i>=0) return m[i]; 
		else 
		{
			fprintf(stderr, "\n-- Cvector(%d): get(%d)\n", nel, i);
			exit(1);
		}
	}
#else
	inline Complex& operator() (const int& i) { return m[i]; }
	inline Complex operator() (const int& i) const { return m[i]; }
#endif

	// set all elements to c
	void init(Complex c);

	// complex conjugate
	void conj();

	// real and imaginary parts
	Dvector re() const;
	Dvector im() const;
	// vectors of absolute values and absolute squares
	Dvector abs() const;
	Dvector sqabs() const;
	// vector containing the argument (phase) of the original entries 
	Dvector arg() const;

	// extract a subvector: n entries starting at position x
	Cvector subvector(int n, int x) const;
	// set the entries beginning at position x to the values of
	// the supplied smaller vector 
	void setsubvector(const Cvector &s, int x);

	// append an entry to the end of the vector 
	void append(Complex e);
	// insert value e at position p,
	// shift the following elements accordingly 
	void insert(Complex e, int p);
	// remove the element at position p, 
	// shift the following elements accordingly
	void remove(int p);

	// the Euclidian norm 
	double norm() const;
	// ... squared 
	double sqnorm() const;

	// maximum absolute value of the vector elements, corresponding index
	double max() const;
	double max(int& mi) const;

	// addition & assignment 
	void addeq(const Cvector& b);
	void addeq(const Dvector& b);
	// subtraction & assignment
	void subeq(const Cvector& b);
	void subeq(const Dvector& b);
	// multiplication & assignment
	void multeq(const Complex& c);
	void multeq(double c);

	/* ----------------------------------------------------- */

	/* multiplications */
	friend Cmatrix mult(const Cmatrix& a, const Cmatrix& b);
	friend Cmatrix mult(const Cmatrix& a, const Dmatrix& b);
	friend Cmatrix mult(const Dmatrix& a, const Cmatrix& b);
	friend Dmatrix mult(const Dmatrix& a, const Dmatrix& b);
	friend Cvector mult(const Cmatrix& a, const Cvector& b);
	friend Cvector mult(const Dmatrix& a, const Cvector& b);
	friend Cvector mult(const Cmatrix& a, const Dvector& b);
	friend Dvector mult(const Dmatrix& a, const Dvector& b);
	friend Cmatrix mult(const Cmatrix& a, const Complex& b);
	friend Cmatrix mult(const Cmatrix& a, double b);
	friend Cmatrix mult(const Dmatrix& a, const Complex& b);
	friend Dmatrix mult(const Dmatrix& a, double b);
	friend Cmatrix mult(const Complex& a, const Cmatrix& b);
	friend Cmatrix mult(const Complex& a, const Dmatrix& b);
	friend Cmatrix mult(double a,  const Cmatrix& b);
	friend Dmatrix mult(double a,  const Dmatrix& b);
	friend Cvector mult(const Cvector& a, const Complex& b);
	friend Cvector mult(const Cvector& a, double b);
	friend Cvector mult(const Dvector& a, const Complex& b);
	friend Dvector mult(const Dvector& a, double b);
	friend Cvector mult(const Complex& a, const Cvector& b);
	friend Cvector mult(const Complex& a, const Dvector& b);
	friend Cvector mult(double a,  const Cvector& b);
	friend Dvector mult(double a,  const Dvector& b);

	/* additions */
	friend Cmatrix add(const Cmatrix& a, const Cmatrix& b);
	friend Cmatrix add(const Cmatrix& a, const Dmatrix& b);
	friend Cmatrix add(const Dmatrix& a, const Cmatrix& b);
	friend Dmatrix add(const Dmatrix& a, const Dmatrix& b);
	friend Cvector add(const Cvector& a, const Cvector& b);
	friend Cvector add(const Dvector& a, const Cvector& b);
	friend Cvector add(const Cvector& a, const Dvector& b);
	friend Dvector add(const Dvector& a, const Dvector& b);

	/* subtraction */
	friend Cmatrix sub(const Cmatrix& a, const Cmatrix& b);
	friend Cmatrix sub(const Cmatrix& a, const Dmatrix& b);
	friend Cmatrix sub(const Dmatrix& a, const Cmatrix& b);
	friend Dmatrix sub(const Dmatrix& a, const Dmatrix& b);
	friend Cvector sub(const Cvector& a, const Cvector& b);
	friend Cvector sub(const Dvector& a, const Cvector& b);
	friend Cvector sub(const Cvector& a, const Dvector& b);
	friend Dvector sub(const Dvector& a, const Dvector& b);

	/* multiplications, 
	   the argument d represents the elements of a square 
	   diagonal matrix */
	friend Dmatrix diagmult(const Dvector& d, const Dmatrix& m);
	friend Dmatrix diagmult(const Dmatrix& m, const Dvector& d);
	friend Dvector diagmult(const Dvector& d, const Dvector& a);
	friend Cmatrix diagmult(const Cvector& d, const Cmatrix& m);
	friend Cmatrix diagmult(const Cmatrix& m, const Cvector& d);
	friend Cvector diagmult(const Cvector& d, const Cvector& a);

	/* scalar products */
	friend Complex scalp(const Cvector& a, const Cvector& b);
	friend Complex scalp(const Cvector& a, const Dvector& b);
	friend Complex scalp(const Dvector& a, const Cvector& b);
	friend double  scalp(const Dvector& a, const Dvector& b);

private:
	Complex *m;
};

/* -------------------------------------------------------------------- */
/* - matrices --------------------------------------------------------- */
/* -------------------------------------------------------------------- */

/* Bmatrix: twodimensional bool array */
class Bmatrix 
{
public:
	// number of rows
	int r;
	// number of columns
	int c;
	// initialize
	Bmatrix();
	Bmatrix(int rows, int columns);
	// destroy
	~Bmatrix();
	// copy constructor 
	Bmatrix(const Bmatrix& s);
	// assignment 
	Bmatrix& operator=(const Bmatrix& s);

	// free allocated memory
	void strip();

	// input
	void read(FILE *dat);
	// output
	void write(FILE *dat) const;

	// output to MATLAB script file name.m
	void mfile(const char *name) const;

	// subscripting
#ifdef SAFE
	bool& operator() (const int& x, const int& y)
	{ 
		if(r>x && c>y && x>=0 && y>=0) return m[x+y*r]; 
		else 
		{
			fprintf(stderr, "\n-- Bmatrix(%d, %d): set(%d, %d)\n", r, c, x, y);
			exit(1);
		}
	}
	bool operator() (const int& x, const int& y) const
	{ 
		if(r>x && c>y && x>=0 && y>=0) return m[x+y*r]; 
		else 
		{
			fprintf(stderr, "\n-- Bmatrix(%d, %d): get(%d, %d)\n", r, c, x, y);
			exit(1);
		}
	}
#else
	inline bool& operator() (const int& x, const int& y) { return m[x+y*r]; }
	inline bool operator() (const int& x, const int& y) const { return m[x+y*r]; }
#endif

	// set all elements to i
	void init(bool b);

	// extract a row or a column
	Bvector row(int x) const;
	Bvector col(int y) const;
	// set a row or a column
	void row(Bvector v, int x);
	void col(Bvector v, int y);

private:
	int dim;
	bool *m;
};

/* Imatrix: twodimensional int array */
class Imatrix 
{
public:
	// number of rows
	int r;
	// number of columns
	int c;
	// initialize
	Imatrix();
	Imatrix(int rows, int columns);
	// destroy
	~Imatrix();
	// copy constructor 
	Imatrix(const Imatrix& s);
	// assignment 
	Imatrix& operator=(const Imatrix& s);

	// free allocated memory
	void strip();

	// input
	void read(FILE *dat);
	// output
	void write(FILE *dat) const;

	// output to MATLAB script file name.m
	void mfile(const char *name) const;

	// subscripting
#ifdef SAFE
	int& operator() (const int& x, const int& y)
	{ 
		if(r>x && c>y && x>=0 && y>=0) return m[x+y*r]; 
		else 
		{
			fprintf(stderr, "\n-- Imatrix(%d, %d): set(%d, %d)\n", r, c, x, y);
			exit(1);
		}
	}
	int operator() (const int& x, const int& y) const
	{ 
		if(r>x && c>y && x>=0 && y>=0) return m[x+y*r]; 
		else 
		{
			fprintf(stderr, "\n-- Imatrix(%d, %d): get(%d, %d)\n", r, c, x, y);
			exit(1);
		}
	}
#else
	inline int& operator() (const int& x, const int& y) { return m[x+y*r]; }
	inline int operator() (const int& x, const int& y) const { return m[x+y*r]; }
#endif

	// set all elements to i
	void init(int i);
	// the unit matrix of dimension n
	void unit(int n);

	// extract a row or a column
	Ivector row(int x) const;
	Ivector col(int y) const;
	// set a row or a column
	void row(Ivector v, int x);
	void col(Ivector v, int y);

	// minimum and maximum of the matrix elements, corresponding indices
	int min() const;
	int min(int& x, int& y) const;
	int max() const;
	int max(int& x, int& y) const;

	// transpose 
	void transpose();

private:
	int dim;
	int *m;
};

/* Dmatrix: twodimensional double array */
class Dmatrix 
{
	friend class Dvector; 
	friend class Cmatrix; 
	friend class Cvector; 
public:
	// number of rows
	int r;
	// number of columns
	int c;
	// initialize
	Dmatrix();
	Dmatrix(int rows, int columns);
	// destroy
	~Dmatrix();
	// copy constructor 
	Dmatrix(const Dmatrix& s);
	// assignment 
	Dmatrix& operator=(const Dmatrix& s);

	// free allocated memory
	void strip();

	// input
	void read(FILE *dat);
	// output
	void write(FILE *dat) const;

	// output to MATLAB script file name.m
	void mfile(const char *name) const;

	// subscripting
#ifdef SAFE
	double& operator() (const int& x, const int& y)
	{ 
		if(r>x && c>y && x>=0 && y>=0) return m[x+y*r]; 
		else 
		{
			fprintf(stderr, "\n-- Dmatrix(%d, %d): set(%d, %d)\n", r, c, x, y);
			exit(1);
		}
	}
	double operator() (const int& x, const int& y) const
	{ 
		if(r>x && c>y && x>=0 && y>=0) return m[x+y*r]; 
		else 
		{
			fprintf(stderr, "\n-- Dmatrix(%d, %d): get(%d, %d)\n", r, c, x, y);
			exit(1);
		}
	}
#else
	inline double& operator() (const int& x, const int& y) { return m[x+y*r]; }
	inline double operator() (const int& x, const int& y) const { return m[x+y*r]; }
#endif

	// set all elements to d
	void init(double d);
	// the unit matrix of dimension n
	void unit(int n);
	// the diagonal matrix with diagonal elements d
	void diag(Dvector d);

	// transpose 
	void transpose();

	// extract a row or a column
	Dvector row(int x) const;
	Dvector col(int y) const;
	// set a row or a column
	void row(Dvector v, int x);
	void col(Dvector v, int y);

	// extract a submatrix: nr rows and nc columns 
	//                      starting at position x, y
	Dmatrix submatrix(int nr, int nc, int x, int y) const;

	// set the entries beginning at position x, y to the values of
	// the supplied smaller matrix
	void setsubmatrix(const Dmatrix &s, int x, int y);

	// minimum and maximum of the matrix elements, corresponding indices
	double min() const;
	double min(int& x, int& y) const;
	double max() const;
	double max(int& x, int& y) const;

	// the inverse matrix, only for square matrices ... 
	void inverse();

	// singular value decomposition 
	// based on: "svdcmp.c"
	//           (C) Copr. 1986-92 Numerical Recipes Software 5.)13. 
	// "this"=A is replaced by a matrix U, with a matrix v and a 
	// vector w set, such that A = U * diag(w) * v^T 
	void svd(Dvector& w, Dmatrix& v); 

	// Cholesky-decomposition of a symmetrical, positive matrix
	// based on: "choldc.c"
	//           (C) Copr. 1986-92 Numerical Recipes Software 5.)13. 
	// return value: 0 - okay, 1 - matrix is numerically not positive
	int choldc();

	// backward substitution following Cholesky-decomposition 
	// calculate x as a solution of this * x = rhs, 
	// r x r - matrix this, left-triangular,
	// r - dim vector x, returned,
	// r - dim vector rhs 
	Dvector bwdsubst(const Dvector& rhs);

	// backward substitution following Cholesky-decomposition 
	// calculate xmat as a solution of this * xmat = rmat, 
	// r x r - matrix this, left-triangular,
	// r x c - matrix xmat, returned,
	// r x c - matrix rmat 
	Dmatrix bwdsubst(const Dmatrix& rmat);

	// backward substitution following Cholesky-decomposition 
	// calculate vector x as a solution of this^t * x = rhs,
	// r x r - matrix this, left-triangular,
	// r - dim vector x, returned,
	// r - dim vector rhs 
	Dvector bwdsubst_t(const Dvector& rhs);

	// solve a linear system with positive matrix by Cholesky decomp.,
	// calculate vector x as a solution of this * x = rhs,
        // r x r - matrix this, positive, destroyed, 
	//                      replaced by a factor of the Chol. dec.
        // r - dim vector x, returned,
        // r - dim vector rhs
        Dvector CDsolve(const Dvector& rhs);

	// calculate eigenvalues / corresp. eigenvectors 
	// of a real symmetric matrix,
	// based on: "tred2.c", "tqli.c", "dpythag.c"
	//           (C) Copr. 1986-92 Numerical Recipes Software 5.)13. 
	// Only the upper (right) part has to be supplied, 
	// the eigenvalues are returned, the matrix is replaced by 
	// its eigenvectors, returned(j) corresponds to the 
	// normalized eigenvector (*this).col(j)
	Dvector eigenv();

	// addition & assignment 
	void addeq(const Dmatrix& b);
	// subtraction & assignment
	void subeq(const Dmatrix& b);
	// multiplication & assignment
	void multeq(double c);

	/* ----------------------------------------------------- */

	/* multiplications */
	friend Cmatrix mult(const Cmatrix& a, const Cmatrix& b);
	friend Cmatrix mult(const Cmatrix& a, const Dmatrix& b);
	friend Cmatrix mult(const Dmatrix& a, const Cmatrix& b);
	friend Dmatrix mult(const Dmatrix& a, const Dmatrix& b);
	friend Cvector mult(const Cmatrix& a, const Cvector& b);
	friend Cvector mult(const Dmatrix& a, const Cvector& b);
	friend Cvector mult(const Cmatrix& a, const Dvector& b);
	friend Dvector mult(const Dmatrix& a, const Dvector& b);
	friend Cmatrix mult(const Cmatrix& a, const Complex& b);
	friend Cmatrix mult(const Cmatrix& a, double b);
	friend Cmatrix mult(const Dmatrix& a, const Complex& b);
	friend Dmatrix mult(const Dmatrix& a, double b);
	friend Cmatrix mult(const Complex& a, const Cmatrix& b);
	friend Cmatrix mult(const Complex& a, const Dmatrix& b);
	friend Cmatrix mult(double a,  const Cmatrix& b);
	friend Dmatrix mult(double a,  const Dmatrix& b);
	friend Cvector mult(const Cvector& a, const Complex& b);
	friend Cvector mult(const Cvector& a, double b);
	friend Cvector mult(const Dvector& a, const Complex& b);
	friend Dvector mult(const Dvector& a, double b);
	friend Cvector mult(const Complex& a, const Cvector& b);
	friend Cvector mult(const Complex& a, const Dvector& b);
	friend Cvector mult(double a,  const Cvector& b);
	friend Dvector mult(double a,  const Dvector& b);

	/* additions */
	friend Cmatrix add(const Cmatrix& a, const Cmatrix& b);
	friend Cmatrix add(const Cmatrix& a, const Dmatrix& b);
	friend Cmatrix add(const Dmatrix& a, const Cmatrix& b);
	friend Dmatrix add(const Dmatrix& a, const Dmatrix& b);
	friend Cvector add(const Cvector& a, const Cvector& b);
	friend Cvector add(const Dvector& a, const Cvector& b);
	friend Cvector add(const Cvector& a, const Dvector& b);
	friend Dvector add(const Dvector& a, const Dvector& b);

	/* subtraction */
	friend Cmatrix sub(const Cmatrix& a, const Cmatrix& b);
	friend Cmatrix sub(const Cmatrix& a, const Dmatrix& b);
	friend Cmatrix sub(const Dmatrix& a, const Cmatrix& b);
	friend Dmatrix sub(const Dmatrix& a, const Dmatrix& b);
	friend Cvector sub(const Cvector& a, const Cvector& b);
	friend Cvector sub(const Dvector& a, const Cvector& b);
	friend Cvector sub(const Cvector& a, const Dvector& b);
	friend Dvector sub(const Dvector& a, const Dvector& b);

	/* multiplications, 
	   the argument d represents the elements of a square 
	   diagonal matrix */
	friend Dmatrix diagmult(const Dvector& d, const Dmatrix& m);
	friend Dmatrix diagmult(const Dmatrix& m, const Dvector& d);
	friend Dvector diagmult(const Dvector& d, const Dvector& a);
	friend Cmatrix diagmult(const Cvector& d, const Cmatrix& m);
	friend Cmatrix diagmult(const Cmatrix& m, const Cvector& d);
	friend Cvector diagmult(const Cvector& d, const Cvector& a);

	/* scalar products */
	friend Complex scalp(const Cvector& a, const Cvector& b);
	friend Complex scalp(const Cvector& a, const Dvector& b);
	friend Complex scalp(const Dvector& a, const Cvector& b);
	friend double  scalp(const Dvector& a, const Dvector& b);

private:
	int dim;
	double *m;
};

/* Cmatrix: twodimensional complex array */
class Cmatrix 
{
	friend class Cvector; 
	friend class Dmatrix; 
	friend class Dvector; 
public:
	// number of rows
	int r;
	// number of columns
	int c;
	// initialize
	Cmatrix();
	Cmatrix(int rows, int columns);
	Cmatrix(const Dmatrix& rp, const Dmatrix& ip);
	// destroy
	~Cmatrix();
	// copy constructor 
	Cmatrix(const Cmatrix& s);
	// assignment 
	Cmatrix& operator=(const Cmatrix& s);

	// free allocated memory
	void strip();

	// input
	void read(FILE *dat);
	// output
	void write(FILE *dat) const;

	// output to MATLAB script file name.m
	void mfile(const char *name) const;

	// subscripting
#ifdef SAFE
	Complex& operator() (const int& x, const int& y)
	{ 
		if(r>x && c>y && x>=0 && y>=0) return m[x+y*r]; 
		else 
		{
			fprintf(stderr, "\n-- Cmatrix(%d, %d): set(%d, %d)\n", r, c, x, y);
			exit(1);
		}
	}
	Complex operator() (const int& x, const int& y) const
	{ 
		if(r>x && c>y && x>=0 && y>=0) return m[x+y*r]; 
		else 
		{
			fprintf(stderr, "\n-- Cmatrix(%d, %d): get(%d, %d)\n", r, c, x, y);
			exit(1);
		}
	}
#else
	inline Complex& operator() (const int& x, const int& y) { return m[x+y*r]; }
	inline Complex operator() (const int& x, const int& y) const { return m[x+y*r]; }
#endif

	// set all elements to c
	void init(Complex c);
	// complex conjugate
	void conj();
	// the unit matrix of dimension n
	void unit(int n);
	// the diagonal matrix with diagonal elements d
	void diag(Cvector d);

	// extract a row or a column
	Cvector row(int x) const;
	Cvector col(int y) const;
	// set a row or a column
	void row(Cvector v, int x);
	void col(Cvector v, int y);

	// real and imaginary parts 
	Dmatrix re() const;
	Dmatrix im() const;
	// matrices of absolute values (squared) of the original entries  
	Dmatrix abs() const;
	Dmatrix sqabs() const;
	// matrix containing the argument (phase) of the original entries 
	Dmatrix arg() const;

	// maximum absolute value of the matrix elements, corresponding indices
	double max() const;
	double max(int& x, int& y) const;

	// transpose 
	void transpose();

	// adjoint 
	void adjoint();

	// the inverse matrix, only for square matrices ... 
	void inverse();

	// Cholesky-decomposition of a Hermitian, positive matrix
	// based on: "choldc.c"
	//           (C) Copr. 1986-92 Numerical Recipes Software 5.)13. 
	// return value: 0 - okay, 1 - matrix is numerically not positive
	int choldc();

	// backward substitution following Cholesky-decomposition 
	// calculate x as a solution of this * x = rhs, 
	// r x r - matrix this, left-triangular,
	// r - dim vector x, returned,
	// r - dim vector rhs 
	Cvector bwdsubst(const Cvector& rhs);

	// backward substitution following Cholesky-decomposition 
	// calculate vector x as a solution of this^t * x = rhs,
	// r x r - matrix this, left-triangular,
	// r - dim vector x, returned,
	// r - dim vector rhs 
	Cvector bwdsubst_t(const Cvector& rhs);

	// solve a linear system with positive matrix by Cholesky decomp.,
	// calculate vector x as a solution of this * x = rhs,
        // r x r - matrix this, positive, destroyed, 
	//                      replaced by a factor of the Chol. dec.
        // r - dim vector x, returned,
        // r - dim vector rhs
        Cvector CDsolve(const Cvector& rhs);

	// calculate eigenvalues / corresp. eigenvectors
	// of a complex Hermitian matrix,
	// the eigenvalues are returned, the matrix is replaced by
	// its eigenvectors, returned(j) corresponds to the
	// normalized eigenvector (*this).col(j)
	// based on converting the Hermitian eigensystem to a real symmetric 
	// system with twice the dimension, 
	// see W.H. Press et. al., "Numerical Recipes in C", chapter 11.4 
	Dvector eigenv();

	// calculate eigenvalues / corresp. eigenvectors
	// of a general complex matrix,
	// the eigenvalues are returned, the matrix is replaced by
	// its eigenvectors, returned(j) corresponds to the
	// normalized eigenvector (*this).col(j)
	// based on the ComplexEigenSolver of the Eigen C++ template library;
	// info: 0: computation successful,
	//       1: erroneous prerequisites,
	//       2: no convergence 
	Cvector geneigenv(int &info);

	// extract a submatrix: nr rows and nc columns 
	//                      starting at position x, y
	Cmatrix submatrix(int nr, int nc, int x, int y) const;

	// set the entries beginning at position x, y to the values of
	// the supplied smaller matrix
	void setsubmatrix(const Cmatrix &s, int x, int y);

	// addition & assignment 
	void addeq(const Cmatrix& b);
	void addeq(const Dmatrix& b);
	// subtraction & assignment 
	void subeq(const Cmatrix& b);
	void subeq(const Dmatrix& b);
	// multiplication & assignment
	void multeq(const Complex& c);
	void multeq(double c);

	/* ----------------------------------------------------- */

	/* multiplications */
	friend Cmatrix mult(const Cmatrix& a, const Cmatrix& b);
	friend Cmatrix mult(const Cmatrix& a, const Dmatrix& b);
	friend Cmatrix mult(const Dmatrix& a, const Cmatrix& b);
	friend Dmatrix mult(const Dmatrix& a, const Dmatrix& b);
	friend Cvector mult(const Cmatrix& a, const Cvector& b);
	friend Cvector mult(const Dmatrix& a, const Cvector& b);
	friend Cvector mult(const Cmatrix& a, const Dvector& b);
	friend Dvector mult(const Dmatrix& a, const Dvector& b);
	friend Cmatrix mult(const Cmatrix& a, const Complex& b);
	friend Cmatrix mult(const Cmatrix& a, double b);
	friend Cmatrix mult(const Dmatrix& a, const Complex& b);
	friend Dmatrix mult(const Dmatrix& a, double b);
	friend Cmatrix mult(const Complex& a, const Cmatrix& b);
	friend Cmatrix mult(const Complex& a, const Dmatrix& b);
	friend Cmatrix mult(double a,  const Cmatrix& b);
	friend Dmatrix mult(double a,  const Dmatrix& b);
	friend Cvector mult(const Cvector& a, const Complex& b);
	friend Cvector mult(const Cvector& a, double b);
	friend Cvector mult(const Dvector& a, const Complex& b);
	friend Dvector mult(const Dvector& a, double b);
	friend Cvector mult(const Complex& a, const Cvector& b);
	friend Cvector mult(const Complex& a, const Dvector& b);
	friend Cvector mult(double a,  const Cvector& b);
	friend Dvector mult(double a,  const Dvector& b);

	/* additions */
	friend Cmatrix add(const Cmatrix& a, const Cmatrix& b);
	friend Cmatrix add(const Cmatrix& a, const Dmatrix& b);
	friend Cmatrix add(const Dmatrix& a, const Cmatrix& b);
	friend Dmatrix add(const Dmatrix& a, const Dmatrix& b);
	friend Cvector add(const Cvector& a, const Cvector& b);
	friend Cvector add(const Dvector& a, const Cvector& b);
	friend Cvector add(const Cvector& a, const Dvector& b);
	friend Dvector add(const Dvector& a, const Dvector& b);

	/* subtraction */
	friend Cmatrix sub(const Cmatrix& a, const Cmatrix& b);
	friend Cmatrix sub(const Cmatrix& a, const Dmatrix& b);
	friend Cmatrix sub(const Dmatrix& a, const Cmatrix& b);
	friend Dmatrix sub(const Dmatrix& a, const Dmatrix& b);
	friend Cvector sub(const Cvector& a, const Cvector& b);
	friend Cvector sub(const Dvector& a, const Cvector& b);
	friend Cvector sub(const Cvector& a, const Dvector& b);
	friend Dvector sub(const Dvector& a, const Dvector& b);

	/* multiplications, 
	   the argument d represents the elements of a square 
	   diagonal matrix */
	friend Dmatrix diagmult(const Dvector& d, const Dmatrix& m);
	friend Dmatrix diagmult(const Dmatrix& m, const Dvector& d);
	friend Dvector diagmult(const Dvector& d, const Dvector& a);

	friend Cmatrix diagmult(const Cvector& d, const Cmatrix& m);
	friend Cmatrix diagmult(const Cmatrix& m, const Cvector& d);
	friend Cvector diagmult(const Cvector& d, const Cvector& a);

	/* scalar products */
	friend Complex scalp(const Cvector& a, const Cvector& b);
	friend Complex scalp(const Cvector& a, const Dvector& b);
	friend Complex scalp(const Dvector& a, const Cvector& b);
	friend double  scalp(const Dvector& a, const Dvector& b);

private:
	int dim;
	Complex *m;
};

/* -------------------------------------------------------------------- */
/* - some linear algebra routines ------------------------------------- */
/* -------------------------------------------------------------------- */

/* solve the real eigenvalue problem 
   a*x = lambda*x
   for all eigenvalues lambda and eigenvector x (columns of x), 
   ordered, lambda(i) decreases with i */
void solveEVproblem(Dmatrix& a, Dvector& lambda, Dmatrix& x);

/* solve the generalized real eigenvalue problem 
   a*x = lambda*b*x  with positive b
   for all eigenvalues lambda and eigenvectors x (columns of x) */ 
void solvegeneralEVproblem(Dmatrix& a, Dmatrix& b, Dvector& lambda, Dmatrix& x);

/* -------------------------------------------------------------------- */

/* write pairs (x, y) to file name.xyf */
void toxyf(const char *name, Dvector x, Dvector y);
/* write triples (x, y.re, y.im) to file name.xyf, for complex y */
void toxyf(const char *name, Dvector x, Cvector y);

/* -------------------------------------------------------------------- */

/* sort a, b in parallel, such that a is ascending */
void sort(Dvector& a, Dvector &b);

/* -------------------------------------------------------------------- */
/* - third order arrays ------------------------------------------------ */
/* -------------------------------------------------------------------- */

/* Barray3: threedimensional bool array */
class Barray3 
{
public:
	// first dimension 
	int d1;
	// second dimension 
	int d2;
	// third dimension 
	int d3;
	// initialize
	Barray3();
	Barray3(int dim1, int dim2, int dim3);
	// destroy
	~Barray3();
	// copy constructor 
	Barray3(const Barray3& s);
	// assignment 
	Barray3& operator=(const Barray3& s);

	// free allocated memory
	void strip();

	// input
	void read(FILE *dat);
	// output
	void write(FILE *dat) const;

	// subscripting
#ifdef SAFE
	bool& operator() (const int& i, const int& j, const int& k) 
	{ 
		if(d1>i && d2>j && d3>k && i>=0 && j>=0 && k>=0) 
			return m[i+j*d1+k*d1d2];
		else 
		{
			fprintf(stderr, 
			        "\n-- Barray3(%d, %d, %d): set(%d, %d, %d)\n", 
			        d1, d2, d3, i, j, k);
			exit(1);
		}
	}
	bool operator() (const int& i, const int& j, const int& k) const 
	{ 
		if(d1>i && d2>j && d3>k && i>=0 && j>=0 && k>=0) 
			return m[i+j*d1+k*d1d2];
		else 
		{
			fprintf(stderr, 
			        "\n-- Barray3(%d, %d, %d): get(%d, %d, %d)\n", 
			        d1, d2, d3, i, j, k);
			exit(1);
		}
	}
#else
	inline bool& operator() (const int& i, const int& j, const int& k) 
		{ return m[i+j*d1+k*d1d2]; }
	inline bool operator() (const int& i, const int& j, const int& k) const 
		{ return m[i+j*d1+k*d1d2]; }
#endif

	// set all elements to i
	void init(bool b);

private:
	int d1d2;
	int dim;
	bool *m;
};


/* Iarray3: threedimensional integer array */
class Iarray3 
{
public:
	// first dimension 
	int d1;
	// second dimension 
	int d2;
	// third dimension 
	int d3;
	// initialize
	Iarray3();
	Iarray3(int dim1, int dim2, int dim3);
	// destroy
	~Iarray3();
	// copy constructor 
	Iarray3(const Iarray3& s);
	// assignment 
	Iarray3& operator=(const Iarray3& s);

	// free allocated memory
	void strip();

	// input
	void read(FILE *dat);
	// output
	void write(FILE *dat) const;

	// subscripting
#ifdef SAFE
	int& operator() (const int& i, const int& j, const int& k) 
	{ 
		if(d1>i && d2>j && d3>k && i>=0 && j>=0 && k>=0) 
			return m[i+j*d1+k*d1d2];
		else 
		{
			fprintf(stderr, 
			        "\n-- Iarray3(%d, %d, %d): set(%d, %d, %d)\n", 
			        d1, d2, d3, i, j, k);
			exit(1);
		}
	}
	int operator() (const int& i, const int& j, const int& k) const 
	{ 
		if(d1>i && d2>j && d3>k && i>=0 && j>=0 && k>=0) 
			return m[i+j*d1+k*d1d2];
		else 
		{
			fprintf(stderr, 
			        "\n-- Iarray3(%d, %d, %d): get(%d, %d, %d)\n", 
			        d1, d2, d3, i, j, k);
			exit(1);
		}
	}
#else
	inline int& operator() (const int& i, const int& j, const int& k) 
		{ return m[i+j*d1+k*d1d2]; }
	inline int operator() (const int& i, const int& j, const int& k) const 
		{ return m[i+j*d1+k*d1d2]; }
#endif

	// set all elements to i
	void init(int b);

private:
	int d1d2;
	int dim;
	int *m;
};


/* Darray3: threedimensional double array */
class Darray3 
{
public:
	// first dimension 
	int d1;
	// second dimension 
	int d2;
	// third dimension 
	int d3;
	// initialize
	Darray3();
	Darray3(int dim1, int dim2, int dim3);
	// destroy
	~Darray3();
	// copy constructor 
	Darray3(const Darray3& s);
	// assignment 
	Darray3& operator=(const Darray3& s);

	// free allocated memory
	void strip();

	// input
	void read(FILE *dat);
	// output
	void write(FILE *dat) const;

	// subscripting
#ifdef SAFE
	double& operator() (const int& i, const int& j, const int& k) 
	{ 
		if(d1>i && d2>j && d3>k && i>=0 && j>=0 && k>=0) 
			return m[i+j*d1+k*d1d2];
		else 
		{
			fprintf(stderr, 
			        "\n-- Darray3(%d, %d, %d): set(%d, %d, %d)\n", 
			        d1, d2, d3, i, j, k);
			exit(1);
		}
	}
	double operator() (const int& i, const int& j, const int& k) const 
	{ 
		if(d1>i && d2>j && d3>k && i>=0 && j>=0 && k>=0) 
			return m[i+j*d1+k*d1d2];
		else 
		{
			fprintf(stderr, 
			        "\n-- Darray3(%d, %d, %d): get(%d, %d, %d)\n", 
			        d1, d2, d3, i, j, k);
			exit(1);
		}
	}
#else
	inline double& operator() (const int& i, const int& j, const int& k) 
		{ return m[i+j*d1+k*d1d2]; }
	inline double operator() (const int& i, const int& j, const int& k) const 
		{ return m[i+j*d1+k*d1d2]; }
#endif

	// set all elements to i
	void init(double b);

private:
	int d1d2;
	int dim;
	double *m;
};


/* Carray3: threedimensional Complex array */
class Carray3 
{
public:
	// first dimension 
	int d1;
	// second dimension 
	int d2;
	// third dimension 
	int d3;
	// initialize
	Carray3();
	Carray3(int dim1, int dim2, int dim3);
	// destroy
	~Carray3();
	// copy constructor 
	Carray3(const Carray3& s);
	// assignment 
	Carray3& operator=(const Carray3& s);

	// free allocated memory
	void strip();

	// input
	void read(FILE *dat);
	// output
	void write(FILE *dat) const;

	// subscripting
#ifdef SAFE
	Complex& operator() (const int& i, const int& j, const int& k) 
	{ 
		if(d1>i && d2>j && d3>k && i>=0 && j>=0 && k>=0) 
			return m[i+j*d1+k*d1d2];
		else 
		{
			fprintf(stderr, 
			        "\n-- Carray3(%d, %d, %d): set(%d, %d, %d)\n", 
			        d1, d2, d3, i, j, k);
			exit(1);
		}
	}
	Complex operator() (const int& i, const int& j, const int& k) const 
	{ 
		if(d1>i && d2>j && d3>k && i>=0 && j>=0 && k>=0) 
			return m[i+j*d1+k*d1d2];
		else 
		{
			fprintf(stderr, 
			        "\n-- Carray3(%d, %d, %d): get(%d, %d, %d)\n", 
			        d1, d2, d3, i, j, k);
			exit(1);
		}
	}
#else
	inline Complex& operator() (const int& i, const int& j, const int& k) 
		{ return m[i+j*d1+k*d1d2]; }
	inline Complex operator() (const int& i, const int& j, const int& k) const 
		{ return m[i+j*d1+k*d1d2]; }
#endif

	// set all elements to i
	void init(Complex b);

private:
	int d1d2;
	int dim;
	Complex *m;
};
#undef SAFE

#endif  // MATRIX_H