/*
author:
	seth from http://www.bierdatenbank.de
tab stops:
	2 (otherwise hardly readable!)
abr.:
	cbr - call-by-reference
*/
#ifndef __seth_matrix_h
#define __seth_matrix_h

#include "vektor.h"

class matrix{
	private:
		int dim_zeilen;								// anzahl(zeilen)
		int dim_spalten;							// anzahl(spalten)
		double **elem;								// komponenten
		int swap(double&, double&);		// tausche zwei double-zahlen
		double max(double, double);		// liefer max(a,b) a,b double
		int max(int, int);						// liefer max(a,b) a,b int
		double min(double, double);		// liefer min(a,b) a,b double
		int min(int, int);						// liefer min(a,b) a,b int
		int sgn(double, int=0);				// liefer sgn(a, sgn(0)) a double (vorzeichenfunktion)
		int sgn(int, int=0);					// liefer sgn(a, sgn(0)) a int (vorzeichenfunktion)
		int pow(int, int=2);					// liefer x^y x,y int, y>=0
	public:
		matrix(int=1, int=1);								// constructor
		matrix(const matrix&);							// copy-constructor
		matrix(double**, int, int);					// constructor fuer import aus array
		matrix(char*, int=0);								// constructor fuer import aus datei
		void export_to_file(char*, int=0);	// export in datei, wahlweise mit oder ohne dimensionen in den ersten beiden zeilen
		matrix(char*, vektor&, int=0);			// constructor fuer import aus datei mit zusaetzlichem vektor
		~matrix();													// destructor
		void deallocate();									// deallokierer; speicher freigeben
		matrix(const vektor&);							// vektor2matrix
		int isIndex(int, int);							// ueberpruefe, ob dimensionen eingehalten
		int set(int, int, double);					// setze komponente
		double get(int, int);								// liefer komponente
		matrix get_submatrix(int, int, int, int);// liefer untermatrix (x1,y1,x2,y2)
		vektor get_spalte(int);							// liefer spalte
		vektor get_zeile(int);							// liefer zeile
		vektor get_diag();									// erzeuge vektor mit diagonalelementen der matrix M
		int get_dim_spalten();							// liefer anzahl(zeilen)
		int get_dim_zeilen();								// liefer anzahl(spalten)
		double max();												// liefer max(M) M matrix
		double min();												// liefer min(M) M matrix
		int count(double=0, double=0);			// liefer anzahl der vorkommen einer bel. zahl; zweiter parameter definiert toleranz
		int nnz();													// liefer anzahl der nicht-null-stellen; also wie dim_zeilen*dim_spalten-count(0,0), nur schneller 
		void identity();										// aendert matrix zu identitaet
		void zero();												// aendert matrix zu nullmatrix
		void diag(double);									// erzeuge diagonalmatrix mit konstantem diag-elem.
		void tridiag(double, double, double);	// erzeuge tridiagonalmatrix mit 3 konstanten diag-elem.
		void hilbert();											// erzeuge hilbertmatrix
		matrix transpose();									// liefer maxtrixtransponierte
		int is_symmetric();									// pruefe, ob matrix symmetrisch
		double norm(int=1);									// matrix-norm
		double pnorm(int=1);								// matrix-p-norm
		double ssn();												// spaltensummennorm (kurz)
		double spaltensummennorm();					// spaltensummennorm
		double zsn();												// zeilensummennorm (kurz)
		double zeilensummennorm();					// zeilensummennorm
		matrix inverse();										// liefer matrixinverse
		vektor solve(const vektor&);				// loese lgs Ax=b
		vektor solve_linksunten(const vektor&); // loese lgs Lx=b
		vektor solve_rechtsoben(const vektor&); // loese lgs Rx=b
		int lr(matrix&, matrix&);						// LR-zerlegung (kurz)
		int LR(matrix&, matrix&);						// LR-zerlegung
		int lu(matrix&, matrix&);						// LU-zerlegung (kurz)
		int LU(matrix&, matrix&);						// LU-zerlegung
		int qr_givens(matrix&, matrix&);		// QR-zerlegung mit givens-rotation liefert 0/1, cbr: Q, R
		int QR_givens(matrix&, matrix&);		// QR-zerlegung mit givens-rotation liefert 0/1, cbr: Q, R
		void householder_corner_subtraction(matrix&, const matrix&);	// wird nur fuer QR_householder benoetigt. evtl. sogar besser private...
		int qr_householder(matrix&, matrix&);					// QR-zerlegung mit householder-spiegelungen liefert 0/1, cbr: Q, R
		int QR_householder(matrix&, matrix&);					// QR-zerlegung mit householder-spiegelungen liefert 0/1, cbr: Q, R
		int qr_givens_iter(matrix&, int=50, int=0);		// QR-verfahren mit givens-rotation, liefert anzahl(iterationen), cbr: A_iteriert; benoetigt: hessenberg-matrix
		int QR_givens_iter(matrix&, int=50, int=0);		// QR-verfahren mit givens-rotation, liefert anzahl(iterationen), cbr: A_iteriert; benoetigt: hessenberg-matrix
		int QR_givens_iter2(matrix&, int=50, int=0);	// QR-verfahren mit givens-rotation, liefert anzahl(iterationen), cbr: A_iteriert; benoetigt: hessenberg-matrix
		int hessenberg_householder(matrix&, int);			// bringt matrix auf hessenberg-gestalt mit householder-spiegelungen liefert aufwand, cbr: hessenbergmatrix
		int hessenberg_householder2(matrix&, int);		// bringt matrix auf hessenberg-gestalt mit householder-spiegelungen liefert aufwand, cbr: hessenbergmatrix
		vektor tikhonov(const vektor&, double, int=0);				// tikhonov-regularisierung (Ax=b)
		matrix tikhonov(const vektor&, const vektor&, int=0);	// tikhonov-regularisierung (Ax=b) fuer mehrere alphas
		double m2s();												// matrix2skalar (kurz)
		double m2skalar();									// matrix2skalar
		vektor m2v();												// matrix2vektor (kurz)
		vektor m2vektor();									// matrix2vektor
		matrix tensorprodukt(const vektor&, const vektor&, int=-1, int=-1, int=-1, int=-1);	// tensorprodukt v*v=m
		matrix mult_transp();								// schnellerer algorithmus fuer M*M.transpose();
		matrix transp_mult();								// schnellerer algorithmus fuer M.transpose()*M;
		matrix operator-();									// unitaerer matrix operator -m
		matrix operator*(const matrix&);		// m*m = m
		matrix operator*(double);						// m*s = m (komponentenweise)
		matrix operator*(const vektor&);		// m*v = m (muss u.u. manuell gecastet werden mit m2v)
		vektor multiply_sub(const vektor&, int, int, int, int, int, char='r');							// berechne teilweise v*m oder m*v
		friend matrix operator*(double, const matrix&);		// s*m = m
		friend matrix operator*(vektor&, const matrix&);	// v*m = m (muss manuell gecastet werden)
		friend matrix operator*(vektor&, vektor&);				// v*v = m (muss u.u. manuell gecastet werden mit m2s) (tensor-/skalarprodukt)
		matrix& operator*=(const matrix&);	// m*=m = m
		matrix& operator*=(double);					// m*=s = m
		matrix  operator+(const matrix&);		// m+m = m
		matrix  operator+(double);					// m+s*Id = m; addiere nur auf der diagonalen!
		matrix& operator+=(const matrix&);	// m+=m = m
		matrix& operator+=(double d);				// m+=s := m+=s*Id = m; addiere nur auf der diagonalen!
		matrix  operator-(const matrix&);		// m-m = m
		matrix  operator-(double);					// m-s*Id = m; subtrahiere nur auf der diagonalen!
		matrix& operator-=(const matrix&);	// m-=m = m
		matrix& operator-=(double d);				// m-=s := m-=s*Id = m; subtrahiere nur auf der diagonalen!
		void ausgabe(int=0);								// gib komponenten auf bildschirm aus; parameter bestimmt ausgabe-art
		ostream& operator<<(ostream&);			// ausgabe-kram
};

inline ostream& operator<<(ostream& strm, matrix& A){
	A.ausgabe();
	return strm;
};

int matrix::swap(double &a, double &b){
	double h;
	h=a;
	a=b;
	b=h;
	return 0;
}

double matrix::max(double a, double b){
	return (a>b)?a:b;
}

int matrix::max(int a, int b){
	return (a>b)?a:b;
}

double matrix::min(double a, double b){
	return (a<b)?a:b;
}

int matrix::min(int a, int b){
	return (a<b)?a:b;
}

int matrix::sgn(double a, int zero){
	int res=(a>0)?1:(a==0)?zero:-1;
	return res;
}

int matrix::sgn(int a, int zero){
	int res=(a>0)?1:(a==0)?zero:-1;
	return res;
}

int matrix::pow(int x, int p){
	if(p<0) cout <<"error: exponent bei int pow(int, int) muss ganzzahlig und >=0 sein"<<endl;
	int res=1;
	for(int i=0;i<p;++i) res*=x;
	return res;
}

// constructor
matrix::matrix(int zeilen, int spalten){
	if(zeilen>0 && spalten>0){
		dim_zeilen=zeilen;
		dim_spalten=spalten;
		elem=new double*[dim_zeilen];
		for(int i=0;i<dim_zeilen;++i){
			elem[i]=new double[dim_spalten];
//			for(int j=0;j<dim_spalten;++j) elem[i][j]=0;
		}
	}else cout <<"error: dimensionen muessen positiv sein";
}

// copy-constructor
matrix::matrix(const matrix& M){
	dim_zeilen=M.dim_zeilen;
	dim_spalten=M.dim_spalten;
	elem=new double*[dim_zeilen];
	for(int i=0;i<dim_zeilen;++i){
		elem[i]=new double[dim_spalten];
		for(int j=0;j<dim_spalten;++j) elem[i][j]=M.elem[i][j];
	}
}

// constructor fuer import aus array
matrix::matrix(double** array, int zeilen, int spalten){
	if(zeilen>0 && spalten>0){
		dim_zeilen=zeilen;
		dim_spalten=spalten;
		elem=new double*[dim_zeilen];
		for(int i=0;i<dim_zeilen;++i){
			elem[i]=new double[dim_spalten];
			for(int j=0;j<dim_spalten;++j) elem[i][j]=array[i][j];
		}
	}else cout <<"error: dimensionen muessen positiv sein";
}

// constructor fuer import aus datei
matrix::matrix(char* filename, int verbose){
	int error=1;
	std::ifstream infile; //nur einlesen
	double komp;
	infile.open(filename,std::ios::in);
	if(infile >>dim_zeilen) // lese dimensionen ein
		if(infile >>dim_spalten)
			if(dim_zeilen>0 && dim_spalten>0) error=0;
	if(!error){
		elem=new double*[dim_zeilen]; // erstelle entsprechende matrix
		for(int i=0;i<dim_zeilen;++i) elem[i]=new double[dim_spalten];
		int row=0, col=0;
		while(infile >>komp){ // einlesen der datei in matrix
			elem[row][col]=komp;
			if(++col==dim_spalten){
				col=0;
				if(verbose) cout <<"\r"<<(dim_zeilen-row);
				if(++row==dim_zeilen) break;
			}
		}
		if(verbose) cout <<"\r "<<endl;
	}else cout <<"datei konnte nicht geoeffnet werden oder enthaelt keine zulaessigen dimensionen in den ersten beiden zeilen";
	infile.close();
}

// export in datei, wahlweise mit oder ohne dimensionen in den ersten beiden zeilen
void matrix::export_to_file(char* filename, int save_dim){
	std::ofstream outfile; //nur schreiben
	outfile.open(filename,std::ios::out);
	if(save_dim){
		outfile <<dim_zeilen<<"\n"; // schreibe dimensionen
		outfile <<dim_spalten<<"\n";
	}
	outfile.setf(ios::scientific,  ios::floatfield);	// scientific-notation
	outfile.setf(ios::adjustfield, ios::right);			// rechtsbuendig
	outfile.setf(ios::showpoint);										// dezimalpunkt und abschließende nullen ausgeben
	outfile.precision(4);														// genauigkeit der gleitkommawerte=2
	for(int row=0;row<dim_zeilen;++row){
		for(int col=0;col<dim_spalten;++col)
			outfile <<setw(13)<<elem[row][col];
		outfile <<"\n";
	}
	outfile.close();
}

// constructor fuer import aus datei mit zusaetzlichem vektor
matrix::matrix(char* filename, vektor& v, int verbose){
	int error=1;
	std::ifstream infile; //nur einlesen
	double komp;
	infile.open(filename,ios::in);
	if(infile >>dim_zeilen) // lese dimensionen ein
		if(infile >>dim_spalten)
			if(dim_zeilen>0 && dim_spalten>0) error=0;
	if(!error){
		elem=new double*[dim_zeilen]; // erstelle entsprechende matrix
		for(int i=0;i<dim_zeilen;++i) elem[i]=new double[dim_spalten];
		int row=0, col=0;
		while(infile >>komp){ // einlesen der datei in matrix
			elem[row][col]=komp;
			if(++col==dim_spalten){
				col=0;
				if(verbose) cout <<"\r"<<(dim_zeilen-row)<<" ";
				if(++row==dim_zeilen) break;
			}
		}
		if(verbose) cout <<"\r "<<endl;
		v.dim=dim_zeilen;
		delete[] v.elem;
		v.elem=new double[v.dim];
		while(infile >>komp && col<dim_zeilen) v.elem[col++]=komp; // einlesen der datei in vektor
	}else cout <<"datei konnte nicht geoeffnet werden oder enthaelt keine zulaessigen dimensionen in den ersten beiden zeilen";
	infile.close();
}

// destructor
matrix::~matrix(){
//	for(int i=0;i<dim_zeilen;++i) delete[] elem[i];
//	delete[] elem;
}

// deallokierer; speicher freigeben
void matrix::deallocate(){
	for(int i=0;i<dim_zeilen;++i) delete[] elem[i];
	delete[] elem;
}

// vektor2matrix-copy constructor
matrix::matrix(const vektor& v){
	if(v.transp==0){
		dim_spalten=1;
		dim_zeilen=v.dim;
		for(int i=0;i<v.dim;++i) elem[i][0]=v.elem[i];
	}else{
		dim_spalten=v.dim;
		dim_zeilen=1;
		for(int i=0;i<v.dim;++i) elem[0][i]=v.elem[i];
	}
}

// pruefe, ob zugriff auf komponente innerhalb der matrix
int matrix::isIndex(int x, int y){
	return ((0<=x) && (x<dim_zeilen) && (0<=y) && (y<dim_spalten))?1:0;
}

// setze komponente auf gegebenen wert
int matrix::set(int x,int y, double d){
	int res=1;
  if(isIndex(x,y)==1) elem[x][y]=d;
	else{
		cout <<"error: komponentenindex > dimension!"<<endl;
		res=0;
	}
	return res;
}

// liefer wert einer gegebenen komponente
double matrix::get(int x, int y){
	if(isIndex(x,y)==1) return elem[x][y];
	else{
		cout <<"error: komponentenindex > dimension!"<<endl;
		return 0;
	}
}

// liefer untermatrix (zeilenbeginn, spaltenbegin, zeilenende, spaltenende)
matrix matrix::get_submatrix(int x1, int y1, int x2, int y2){
	matrix S(abs(x2-x1)+1, abs(y2-y1)+1);
	if(isIndex(x1,y1) && isIndex(x2,y2) && x1<=x2 && y1<=y2){
		for(int i=0;i<S.dim_zeilen;++i)
			for(int j=0;j<S.dim_spalten;++j)
				S.elem[i][j]=elem[i+x1][j+y1];
	}
	else{
		S.zero();
		cout <<"error in get_submatrix(): untermatrix laesst sich nicht erstellen!"<<endl;
		if(!isIndex(x1,y1)) cout <<"!isIndex(x1,y1)"<<endl;
		if(!isIndex(x2,y2)) cout <<"!isIndex(x2,y2)"<<endl;
		if(x1>x2) cout <<"x1>x2"<<endl;
		if(y1>y2) cout <<"y1>y2"<<endl;
	}
	return S;
}

// liefer spalte
vektor matrix::get_spalte(int col){
	vektor v(dim_zeilen);
	for(int i=0;i<dim_zeilen;++i) v.elem[i]=elem[i][col];
	return v;
}

// liefer zeile
vektor matrix::get_zeile(int row){
	vektor v(dim_spalten);
	for(int i=0;i<dim_spalten;++i) v.elem[i]=elem[row][i];
	v.transpose();
	return v;
}

// erzeuge vektor mit diagonalelementen der matrix M
vektor matrix::get_diag(){
	int dim=min(dim_zeilen, dim_spalten);
	vektor v(dim);
	for(int i=0;i<dim;++i) v.elem[i]=elem[i][i];
	return v;
}

// liefer anzahl der spalten
int matrix::get_dim_spalten(){
	return dim_spalten;
}

// liefer anzahl der zeilen
int matrix::get_dim_zeilen(){
	return dim_zeilen;
}

// liefer betragsgroesste komponente
double matrix::max(){
	double res=0;
	for(int i=0;i<dim_zeilen;++i)
		for(int j=0;j<dim_spalten;++j)
			res=max(res, fabs(elem[i][j]));
	return res;
}

// liefer betragskleinste komponente
double matrix::min(){
	double res=DBL_MAX;
	for(int i=0;i<dim_zeilen;++i)
		for(int j=0;j<dim_spalten;++j)
			res=min(res, fabs(elem[i][j]));
	return res;
}

// liefer anzahl der vorkommen einer bel. zahl; zweiter parameter definiert toleranz
int matrix::count(double x, double tol){
	int res=0;
	for(int i=0;i<dim_zeilen;++i)
		for(int j=0;j<dim_spalten;++j)
			res+=(fabs(elem[i][j]-x)<=tol);
	return res;
}

// liefer anzahl der nicht-null-stellen; also wie dim_zeilen*dim_spalten-count(0,0), nur schneller 
int matrix::nnz(){
	int res=0;
	for(int i=0;i<dim_zeilen;++i)
		for(int j=0;j<dim_spalten;++j)
			res+=(elem[i][j]!=0);
	return res;
}

// aendere matrix zu identitaet ab
void matrix::identity(){
	for(int i=0;i<dim_zeilen;++i)
		for(int j=0;j<dim_spalten;++j)
			elem[i][j]=(i==j);
}

// aendere matrix zu nullmatrix ab
void matrix::zero(){
	for(int i=0;i<dim_zeilen;++i)
		for(int j=0;j<dim_spalten;++j)
			elem[i][j]=0;
}

// aendere matrix zu diagonalmatrix ab: diag(d)
void matrix::diag(double d){
	for(int i=0;i<dim_zeilen;++i){
	  if(isIndex(i,i)==1) elem[i][i]=d;
	}
}

// aendere matrix zu tridiagonalmatrix ab: tridiag(x,y,z)
void matrix::tridiag(double x, double y, double z){
	for(int i=0;i<dim_zeilen;++i){
	  if(isIndex(i,i-1)==1) elem[i][i-1]=x;
	  if(isIndex(i,i)==1) elem[i][i]=y;
	  if(isIndex(i,i+1)==1) elem[i][i+1]=z;
	}
}

// aendere matrix zu hilbertmatrix ab: hilbert()
void matrix::hilbert(){
	for(int i=0;i<dim_zeilen;++i)
		for(int j=0;j<dim_spalten;++j)
			elem[i][j]=1.0/(i+j+1); // eigentlich -1 aber wir fangen ja bei 0 an zu zaehlen
}

// liefer transponierte matrix
matrix matrix::transpose(){
	matrix H(dim_spalten, dim_zeilen);
	for(int i=0;i<H.dim_zeilen;++i)
		for(int j=0;j<H.dim_spalten;++j)
			H.elem[i][j]=elem[j][i];
	return H;
}

// pruefe, ob matrix symmetrisch
int matrix::is_symmetric(){
	int res=1;
	for(int i=0;i<dim_zeilen && res==1;++i)
		for(int j=0;j<dim_spalten && res==1;++j)
			if(elem[i][j]!=elem[j][i]) res=0;
	return res;
}

// berechne p-norm (geht momentan nur fuer p=1)
double matrix::norm(int p){
	return (*this).pnorm(p);
}

// berechne p-norm (geht momentan nur fuer p=1)
double matrix::pnorm(int p){
	double res=0;
	if(p==1) res=(*this).spaltensummennorm();
	else cout <<"error: noch nicht implementiert"<<endl;
	return res;
}

// berechne spaltensummennorm
double matrix::ssn(){
	return (*this).spaltensummennorm();
}

// berechne spaltensummennorm
double matrix::spaltensummennorm(){
	double res=0, sum;
	for(int i=0;i<dim_spalten;++i){
		sum=0;
		for(int j=0;j<dim_zeilen;++j) sum+=fabs(elem[j][i]);
		res=max(res,sum);
	}
	return res;
}

// berechne zeilensummennorm
double matrix::zsn(){
	return (*this).zeilensummennorm();
}

// berechne zeilensummennorm
double matrix::zeilensummennorm(){
	double res=0, sum;
	for(int i=0;i<dim_zeilen;++i){
		sum=0;
		for(int j=0;j<dim_spalten;++j) sum+=fabs(elem[i][j]);
		res=max(res,sum);
	}
	return res;
}

// inverts nxn-matrix with Gauss-Jordan-algorithm with column pivoting
matrix matrix::inverse(){
	matrix H(*this);
	if(H.dim_zeilen==H.dim_spalten){
		int dim=H.dim_zeilen;
		int col,row,el;
		int *v;
		v=new int[dim];
		double b,max;
		for(col=0;col<dim;++col){
			v[col]=col;
			max=0;
			for(row=col;row<dim;++row){	// suche max:=betragsgroesstes elem. in spalte col, ab col-tem element
				if(fabs(H.elem[row][col])>fabs(max)){
					max=H.elem[row][col];
					el=row;									// max = el-te komponente in col-ter spalte
				}
			}
			if(max==0){									// (nullspalte =>) matrix singulaer
				cout <<"error: matrix nicht invbar"<<endl;
				return (*this);
			}
			if(el!=col){								// wenn max nicht ganz oben (in col-ter colxcol-hauptuntermatrix) steht vertausche zeilen, so dass max ganz oben steht also =a[col][col]
				for(int j=0;j<dim;++j) swap(H.elem[col][j], H.elem[el][j]);
				v[col]=el;								// speichere nummer der gewechselten zeile
			}
			H.elem[col][col]=1;					// setze max auf 1
			for(int j=0;j<dim;++j) H.elem[col][j]/=max;// normalisiere zeile col (mit 1/max)
			for(row=0;row<dim;++row){
				if(row!=col){							// in allen anderen zeilen
					b=H.elem[row][col];
					H.elem[row][col]=0;			// anulliere spalte col
					for(int j=0;j<dim;++j) H.elem[row][j]-=b*H.elem[col][j];
				}
			}
		}
		for(col=dim-1;col>=0;--col){
			if(v[col]!=col) for(int j=0;j<dim;++j) swap(H.elem[j][col], H.elem[j][v[col]]);
		}
	}else cout <<"error: matrix nicht quadratisch => nicht invbar"<<endl;
	return H;
}

// loese Ax=b
vektor matrix::solve(const vektor &b){
	vektor x(b.dim);
	x=((*this).inverse()*b).m2v();
	return x;
}

// loese Ax=b, mit A=L
vektor matrix::solve_linksunten(const vektor &b){
	matrix H(*this);
	vektor x(b);
	for(int j=0;j<dim_spalten;++j){
		if(elem[j][j]==0){
			cout <<"error: division by zero (vektor matrix::solve_linksunten(const vektor &b))"<<endl;
			x.elem[j]=DBL_MAX; // ui, ui, jetzt wuerde alles falsch werden... ;-)
		}
		else x.elem[j]/=elem[j][j];
		H.elem[j][j]=1;
		for(int i=j+1;i<dim_zeilen;++i){
			x.elem[i]-=elem[i][j]*x.elem[j];
			H.elem[i][j]=0;
		}
	}
	H.deallocate();
	return x;
}

// loese Ax=b, mit A=R
vektor matrix::solve_rechtsoben(const vektor &b){
	matrix H(*this);
	vektor x(b);
	for(int j=dim_spalten-1;j>=0;--j){
		if(elem[j][j]==0){
			cout <<"error: division by zero (vektor matrix::solve_rechtsoben(const vektor &b))"<<endl;
			x.elem[j]=DBL_MAX; // ui, ui, jetzt wuerde alles falsch werden... ;-)
		}
		else x.elem[j]/=elem[j][j];
		H.elem[j][j]=1;
		for(int i=j-1;i>=0;--i){
			x.elem[i]-=elem[i][j]*x.elem[j];
			H.elem[i][j]=0;
		}
	}
	H.deallocate();
	return x;
}

// LR-zerlegung
int matrix::lr(matrix &L, matrix &R){
	return (*this).LR(L, R);
}

// LR-zerlegung
int matrix::LR(matrix &L, matrix &R){
	int res=0;
	matrix H(*this);
	if(H.dim_zeilen==H.dim_spalten){
		int dim=H.dim_zeilen;
		for(int i=1;i<dim;++i) H.elem[i][0]/=H.elem[0][0];
		for(int i=1;i<dim-1;++i){
			// bestimme i-te zeile der matrix R
			for(int j=i;j<dim;++j)
			  for(int k=0;k<=i-1;++k) H.elem[i][j]-=H.elem[i][k]*H.elem[k][j];
			// bestimme i-te spalte der matrix L
			for(int j=i+1;j<dim;++j){
				for(int k=0;k<=i-1;++k) H.elem[j][i]-=H.elem[j][k]*H.elem[k][i];
				H.elem[j][i]/=H.elem[i][i];
			}
		}
		// bestimmung des wertes rnn
  	for(int k=0;k<dim-1;++k) H.elem[dim-1][dim-1]-=H.elem[dim-1][k]*H.elem[k][dim-1];
  	res=1;
	}else cout <<"error: matrix nicht quadratisch => nicht LR-zerlegbar"<<endl;
	for(int i=0;i<H.dim_zeilen;++i){
		for(int j=0;j<i;++j){
			L.elem[i][j]=H.elem[i][j];
			R.elem[i][j]=0;
		}
		for(int j=i;j<L.dim_spalten;++j){
			R.elem[i][j]=H.elem[i][j];
			L.elem[i][j]=0;
		}
		L.elem[i][i]=1;
	}
	H.deallocate();
	return res;
}

// LU-zerlegung
int matrix::lu(matrix &L, matrix &U){
	return (*this).LU(L, U);
}

// LU-zerlegung
int matrix::LU(matrix &L, matrix &U){
	int res=0;
	double sum;
	if(dim_zeilen==dim_spalten){
		int dim=dim_zeilen;
		for(int i=0;i<dim;++i){
			for(int j=0;j<=i;++j){
				sum=0;
				for(int k=0;k<j;++k) sum+=L.elem[i][k]*U.elem[k][j];
				L.elem[i][j]=elem[i][j]-sum;
			}
			for(int j=i;j<dim;++j){
				sum=0;
				for(int k=0;k<i;++k) sum+=L.elem[i][k]*U.elem[k][j];
				U.elem[i][j]=(elem[i][j]-sum)/L.elem[i][i];
			}
		}
  	res=1;
	}else cout <<"error: matrix nicht quadratisch => nicht LU-zerlegbar"<<endl;
	return res;
}

// QR-zerlegung mit givens-rotationen
int matrix::qr_givens(matrix &Q, matrix &R){
	return (*this).QR_givens(Q, R);
}

// QR-zerlegung mit givens-rotationen
int matrix::QR_givens(matrix &Q, matrix &R){
	matrix R_k(dim_zeilen, dim_spalten);
	double x, y, xydist, negsine, cosine;
	// initialize the matrices
	Q.identity();
	R=(*this);
	// find each of the Givens rotation matrices
	for(int col=0;col<dim_spalten-1;++col){
		for(int row=col+1;row<dim_zeilen;++row){
			// only compute a rotation matrix if the corresponding entry in the matrix is nonzero
			y=R.elem[row][col];
			if(y!=0.0){
				// find the values of -sin and cos
				x=R.elem[col][col];
				xydist=sqrt(x*x + y*y);
				negsine=y/xydist;
				cosine=x/xydist;
				// calculate the R2 and R^t matrices
				R_k.identity();
				R_k.elem[col][col]=cosine;
				R_k.elem[col][row]=negsine;
				R_k.elem[row][col]=-negsine;
				R_k.elem[row][row]=cosine;
				// "add" R2 to the Q and R matrices
				Q*=R_k.transpose();
				R=R_k*R;
			}
		}
	}
	R_k.deallocate();
	return 1;
}

// QR-zerlegung mit householder-spiegelungen
int matrix::qr_householder(matrix &Q, matrix &R){
	return (*this).QR_householder(Q, R);
}

void matrix::householder_corner_subtraction(matrix &left_large, const matrix &right_small){
	/* calculate start offsets for the large matrix */
	int row_offset=left_large.dim_zeilen-right_small.dim_zeilen;
	int col_offset=left_large.dim_spalten-right_small.dim_spalten;
	/* subtract the elements individually */
	for(int row=0;row<right_small.dim_zeilen;++row)
		for(int col=0;col<right_small.dim_spalten;++col)
			left_large.elem[row_offset+row][col_offset+col]-= right_small.elem[row][col];
}

// QR-zerlegung mit householder-spiegelungen
int matrix::QR_householder(matrix &Q, matrix &R){
	matrix H(*this);
	// initialize the matrices
	Q.identity();
	R=(*this);
	// find each of the Householder reflection matrices
	for(int col=0;col<dim_spalten-1;++col){
		// get the x vector
		matrix x(dim_spalten - col, 1);
		x=R.get_submatrix(col, col, dim_spalten-1, col);
		// make x into u~ (it only involves changing one element, so there's no point in allocating memory for a brand new vector)
		x.elem[0][0]+=sgn(x.elem[0][0], 1)*x.m2v().norm();
		// create H_temp (the 2uu^t part of H)
		matrix H_temp(x.dim_spalten, x.dim_spalten);
		H_temp=x*x.transpose()*(2.0/((x.transpose()*x).m2s()));
		// create H (I - 2uu^t)
		H.identity();
		householder_corner_subtraction(H, H_temp);
		// "add" h to the q and r matrices
		Q*=H;
		R=H*R;
		H_temp.deallocate();
		x.deallocate();
	}
	H.deallocate();
	return 1;
}

// QR-verfahren mit givens einer hessenbergmatrix
int matrix::qr_givens_iter(matrix &A_k, int maxIT, int verbose){
	return (*this).QR_givens_iter(A_k, maxIT, verbose);
}

// QR-verfahren mit givens einer hessenbergmatrix liefert #iterationsschritte
int matrix::QR_givens_iter(matrix &A_k, int maxIT, int verbose){
	int k=0, n=A_k.dim_zeilen-1;
	double x, y, xydist, s, c, b_i, b_i2;
	double epsilon=1.e-16; // werte fuer abbruch-kriterium
	double A_1, A_2, pi=3.1415926535897932384626433832795; // nur fuer ausgabe
	while(k<(maxIT+1) && (min(fabs(A_k.elem[n][n-1]), fabs(A_k.elem[n-1][n-2])) >= epsilon*(fabs(A_k.elem[n][n])+fabs(A_k.elem[n-1][n-1])))){
		for(int i=0;i<A_k.dim_spalten-1;++i){
			// ermittle drehkaestchen
			x=A_k.elem[i][i];
			y=A_k.elem[i+1][i];
			xydist=sqrt(x*x+y*y);
			c=x/xydist;
			s=y/xydist;
			for(int j=0;j<A_k.dim_zeilen;++j){	// multipliziere G (=givens-rot) von links
				b_i=A_k.elem[i][j];
				b_i2=A_k.elem[i+1][j];
				A_k.elem[i][j]=c*b_i+s*b_i2;
				A_k.elem[i+1][j]=-s*b_i+c*b_i2;
			}
			for(int j=0;j<A_k.dim_spalten;++j){	// multipliziere G (=givens-rot) von rechts
				b_i=A_k.elem[j][i];
				b_i2=A_k.elem[j][i+1];
				A_k.elem[j][i]=c*b_i+s*b_i2;
				A_k.elem[j][i+1]=-s*b_i+c*b_i2;
			}
		}
		if(verbose && k<=10){
			cout <<((k<10)?"0":"")<<k<<".   ";
			printf("%5.4f", fabs(A_k.elem[1][0]));
			cout <<"   ";
			A_1=2*(1-cos(10.0/11*pi));
			A_2=2*(1-cos(9.0/11*pi));
			printf("%5.4f", std::pow(fabs(A_2/A_1), k));
			cout <<endl;
		}
		++k;
	}
	return k;
}

// QR-verfahren mit givens einer hessenbergmatrix liefert #iterationsschritte
int matrix::QR_givens_iter2(matrix &A_k, int maxIT, int verbose){
	if(A_k.dim_zeilen!=A_k.dim_spalten) cout <<"error in matrix::QR_givens_iter(...): matrix nicht quadratisch"<<endl;
	int it=0;
	double x, y, xydist, s, c, b_i, b_i2;
	double tol=1.e-8; // wert fuer abbruch-kriterium
	double *d, *d1, *d2;
	int anzahl_kappa=4, kappa_max=15;
	d=new double[anzahl_kappa];
	d1=new double[anzahl_kappa];
	d2=new double[anzahl_kappa];
	matrix kappa(kappa_max, anzahl_kappa); kappa.zero();
	for(int i=A_k.dim_zeilen;i>1;--i){
		while(it<=maxIT && fabs(A_k.elem[i-1][i-2])>=tol){
			for(int j=0;j<anzahl_kappa;++j){ //eoc
				d2[j]=d1[j];
				d1[j]=d[j];
				d[j]=(A_k.isIndex(j+1,j))?A_k.elem[j+1][j]:0;
			}
			for(int j=0;j<A_k.dim_spalten-1;++j){
				// ermittle drehkaestchen
				x=A_k.elem[j][j];
				y=A_k.elem[j+1][j];
				xydist=sqrt(x*x+y*y);
				c=x/xydist;
				s=y/xydist;
				for(int k=0;k<i;++k){	// multipliziere G (=givens-rot) von links
					b_i=A_k.elem[j][k];
					b_i2=A_k.elem[j+1][k];
					A_k.elem[j][k]=c*b_i+s*b_i2;
					A_k.elem[j+1][k]=-s*b_i+c*b_i2;
				}
				for(int k=0;k<i;++k){	// multipliziere G (=givens-rot) von rechts
					b_i=A_k.elem[k][j];
					b_i2=A_k.elem[k][j+1];
					A_k.elem[k][j]=c*b_i+s*b_i2;
					A_k.elem[k][j+1]=-s*b_i+c*b_i2;
				}
			}
			for(int j=0;j<anzahl_kappa;++j) d[j]=(A_k.isIndex(j+1,j))?fabs(A_k.elem[j+1][j]-d[j]):0;
			if(it>1){
				if(verbose) cout <<it<<". schritt: d[0]="<<d[0]<<endl;
				for(int j=0;j<anzahl_kappa;++j){
					if(d2[j]==0 || d1[j]==0 || d1[j]==d2[j]){}// cout <<"error: division by zero (berechnung von kappa)"<<endl;
					else kappa.elem[it][j]=log(d[j]/d1[j])/log(d1[j]/d2[j]);
				}
			}
			++it;
		}
	}
	delete[] d, d1, d2;
	if(verbose) cout <<"kappa = "<<endl;
	kappa.ausgabe(1);
	kappa.deallocate();
	return it;
}

// bringt matrix auf hessenberg-gestalt mit householder-spiegelungen liefert aufwand, cbr: hessenbergmatrix
int matrix::hessenberg_householder(matrix &HB, int verbose){
	if(dim_zeilen!=dim_spalten) cout <<"error in int matrix::hessenberg_householder(matrix), dim_zeilen!=dim_spalten"<<endl;
	int dim=dim_zeilen;
	HB=(*this);
	int aufwand=0;
	vektor v(dim), w(dim);
	double v_faktor;
	for(int i=0;i<dim-2;++i){
		if(verbose) cout <<(i+1)<<". schritt"<<endl;
		for(int j=i+1;j<dim;++j) v.elem[j]=HB.elem[j][i]; // wichtig von v sind nur die komponenten i+1..dim-1, rest ist abfall, wird i.f. ignoriert
		v_faktor=0; // erst mal norm, spaeter norm^2, noch spaeter 2/norm^2
		for(int j=i+2;j<dim;++j) v_faktor+=v.elem[j]*v.elem[j]; // |v|
		aufwand+=dim-i-1;
		v.elem[i+1]+=sgn(v.elem[i+1], 1)*sqrt(v_faktor+v.elem[i+1]*v.elem[i+1]);
		aufwand+=3; // wurzel wird mit aufwand 1 gewertet, obwohl eigentlich >1...
		if(verbose) cout <<"v="<<v;
		v_faktor=2.0/(v_faktor+v.elem[i+1]*v.elem[i+1]); // 2/(|v|^2)
		aufwand+=3;

		// QA
		w=HB.multiply_sub(v, i+1, i+1, i+1, dim-1, dim-1, 'l'); // berechnet teil von v*A.
		aufwand+=(dim-i-1)*(dim-i-1);
		for(int j=i+1;j<dim;++j) w.elem[i]+=v.elem[j]*HB.elem[j][i];
		aufwand+=(dim-i-1);
		for(int j=i;j<dim;++j) w.elem[j]*=v_faktor;
		aufwand+=(dim-i);
		HB-=tensorprodukt(v, w, i+1, i, dim-1, dim-1);
		if(verbose){
			cout <<"tensorprodukt="<<endl;
			(tensorprodukt(v, w, i+1, i, dim-1, dim-1)).ausgabe();
		}
		aufwand+=(dim-i-1)*(dim-i);
		if(verbose){
			cout <<"QA="<<endl;
			HB.ausgabe();
		}

		// (QA)*Qt
		w=HB.multiply_sub(v, i+1, i+1, i+1, dim-1, dim-1, 'r'); // berechnet teil von A*v.
		aufwand+=(dim-i-1)*(dim-i-1);
		for(int j=i+1;j<dim;++j) w.elem[i]+=HB.elem[i][j]*v.elem[j];
		aufwand+=(dim-i-1);
		for(int j=i;j<dim;++j) w.elem[j]*=v_faktor;
		aufwand+=(dim-i);
		HB-=tensorprodukt(w, v, i, i+1, dim-1, dim-1);
		if(verbose){
			cout <<"tensorprodukt="<<endl;
			(tensorprodukt(w, v, i, i+1, dim-1, dim-1)).ausgabe();
		}
		aufwand+=(dim-i-1)*(dim-i);
		if(verbose){
			cout <<"QAQt="<<endl;
			HB.ausgabe();
		}

		verbose=min((i<3), verbose);
	}
	v.deallocate();
	w.deallocate();
	return aufwand;
}

// bringt matrix auf hessenberg-gestalt mit householder-spiegelungen liefert aufwand, cbr: hessenbergmatrix
int matrix::hessenberg_householder2(matrix &HB, int verbose){
	if(dim_zeilen!=dim_spalten) cout <<"error in int matrix::hessenberg_householder(matrix), dim_zeilen!=dim_spalten"<<endl;
	int dim=dim_zeilen;
	HB=(*this);
	int aufwand=0;
	vektor v(dim);
	double v_faktor, sum;
	matrix Q(dim, dim); Q.zero();
	for(int i=0;i<dim-2;++i){
		if(verbose) cout <<(i+1)<<". schritt"<<endl;
		for(int j=i+1;j<dim;++j) v.elem[j]=HB.elem[j][i];
		v_faktor=0; // erst mal norm, spaeter norm^2, noch spaeter 2/norm^2
		for(int j=i+2;j<dim;++j) v_faktor+=v.elem[j]*v.elem[j]; // |v|
		aufwand+=dim-i-3;
		v.elem[i+1]+=sgn(v.elem[i+1], 1)*sqrt(v_faktor+v.elem[i+1]*v.elem[i+1]);
		aufwand+=2;
		if(verbose) cout <<"v="<<v;
		v_faktor+=v.elem[i+1]*v.elem[i+1]; // |v|^2
		++aufwand;
		v_faktor=2.0/v_faktor;
		for(int j=0;j<dim;++j){
			Q.elem[j][i]=(j==i);
			Q.elem[i][j]=(j==i);
		}
		for(int j=i+1;j<dim;++j)
			for(int k=i+1;k<dim;++k)
				Q.elem[j][k]=(j==k)-v.elem[j]*v.elem[k]*v_faktor;
		aufwand+=(dim-i-2)*(dim-i-2)*2;
		if(verbose){
			cout <<"Q="<<endl;
			Q.ausgabe();
		}
		matrix A(HB.get_submatrix(i,i,dim-1,dim-1));
		for(int j=i;j<dim;++j){
			for(int k=i;k<dim;++k){
				HB.elem[j][k]=0;
				for(int l=i;l<dim;++l){
					sum=0;
					for(int m=i;m<dim;++m)
						sum+=Q.elem[j][m]*A.elem[m-i][l-i];
					HB.elem[j][k]+=sum*Q.elem[k][l];
				}
			}
		}
		aufwand+=pow(dim-i-1,4)+pow(dim-i-1,3);
		A.deallocate();
		if(verbose){
			cout <<"QAQt="<<endl;
			HB.ausgabe();
		}
		verbose=min((i<3), verbose);
	}
	v.deallocate();
	Q.deallocate();
	return aufwand;
}

// tikhonov-regularisierung (Ax=b)
vektor matrix::tikhonov(const vektor& b, double alpha, int verbose){
	clock_t t1, t2;
	vektor result(dim_spalten);
	if(verbose){
		cout <<"tikhonov:"<<endl;	cout <<"At*A+alpha: ";
		t1=clock();
		matrix B(((*this).transp_mult())+alpha);
		t2=clock();
		cout <<(((t2-t1)*1000)/CLOCKS_PER_SEC)<<" msec"<<endl;
		cout <<"B=inverse: ";
		t1=clock();
		B=B.inverse();
		t2=clock();
		cout <<(((t2-t1)*1000)/CLOCKS_PER_SEC)<<" msec"<<endl;
		cout <<"B*(At*b): ";
		t1=clock();
		result=(B*((*this).transpose()*b)).m2v();
		t2=clock();
		cout <<(((t2-t1)*1000)/CLOCKS_PER_SEC)<<" msec"<<endl;
		B.deallocate();
	}else{
		result=((((*this).transp_mult())+alpha).inverse()*((*this).transpose()*b)).m2v();
	}
	return result;
}

// tikhonov-regularisierung (Ax=b) fuer mehrere alphas
matrix matrix::tikhonov(const vektor& b, const vektor& alpha, int verbose){
	clock_t t1, t2;
	if(verbose){cout <<"tikhonov:"<<endl;	cout <<"B=At*A: "; t1=clock();}
	matrix B(((*this).transp_mult()));
	if(verbose){t2=clock(); cout <<(((t2-t1)*1000)/CLOCKS_PER_SEC)<<" msec"<<endl;}
	if(verbose){cout <<"D=At*b: "; t1=clock();}
	matrix D((*this).transpose()*b); // D ist eigentlich vektor (casten kostet aber nur zeit...)
	if(verbose){t2=clock(); cout <<(((t2-t1)*1000)/CLOCKS_PER_SEC)<<" msec"<<endl;}
	matrix C(B.dim_zeilen, B.dim_spalten);
	if(verbose) cout <<"loese LGS fuer alpha[0.."<<(alpha.dim-1)<<"]"<<endl;
	matrix result(alpha.dim, dim_spalten);
	for(int i=0;i<alpha.dim;++i){
		if(verbose){cout <<"C=inv(B+alpha["<<i<<"]), mit alpha["<<i<<"]="<<alpha.elem[i]<<": "; t1=clock();}
		C=(B+alpha.elem[i]).inverse();
		if(verbose){t2=clock(); cout <<(((t2-t1)*1000)/CLOCKS_PER_SEC)<<" msec"<<endl;}
		C*=D;
		for(int j=0;j<dim_spalten;++j) result.elem[i][j]=C.elem[j][0];
	}
	B.deallocate();
	C.deallocate();
	D.deallocate();
	return result;
}

// matrix2skalar (extrahiere erste komponente als double skalar)
double matrix::m2s(){
	return (*this).m2skalar();
}

// matrix2skalar (extrahiere erste komponente als double skalar)
double matrix::m2skalar(){
	double s=elem[0][0];
	if(dim_spalten!=1 || dim_zeilen!=1) cout <<"error: bei matrix2skalar wird mist gebaut: matrix nicht 1x1"<<endl;
	return s;
}

// matrix2vektor (extrahiere erste spalte/zeile als vektor)
vektor matrix::m2v(){
	return (*this).m2vektor();
}

// matrix2vektor (extrahiere erste spalte/zeile als vektor)
vektor matrix::m2vektor(){
	vektor v(max(dim_spalten, dim_zeilen));
	v.transp=(dim_spalten>dim_zeilen);
	for(int i=0;i<v.dim;++i) v.elem[i]=(v.transp)?elem[0][i]:elem[i][0];
	if(min(dim_spalten, dim_zeilen)!=1) cout <<"error: bei matrix2vektor wird mist gebaut: matrix nicht eindimensional"<<endl;
	return v;
}

// tensorprodukt: v*v=m
matrix matrix::tensorprodukt(const vektor& v_left, const vektor& v_right, int v_l_begin, int v_r_begin, int v_l_end, int v_r_end){
	matrix H(v_left.dim, v_right.dim);
	if(v_l_begin==-1 && v_r_begin==-1 && v_l_end==-1 && v_r_end==-1){
		v_l_begin=0;
		v_r_begin=0;
		v_l_end=v_left.dim-1;
		v_r_end=v_right.dim-1;
	}
	int dist_row=v_l_end-v_l_begin;
	int dist_col=v_r_end-v_r_begin;
	if(dist_row<0 || dist_col<0 || v_r_end>=v_right.dim || v_l_end>=v_left.dim)
		cout <<"error in matrix::tensorprodukt(const vektor&, const vektor&, int, int, int, int): zugriff ausserhalb der dimensionen oder so"<<endl;
	else{
		for(int i=v_l_begin;i<=v_l_end;++i)
			for(int j=v_r_begin;j<=v_r_end;++j)
				H.elem[i][j]=v_left.elem[i]*v_right.elem[j];
	}
	return H;
}

// schnellerer algorithmus fuer M*M.transpose();
matrix matrix::mult_transp(){
	matrix H(dim_zeilen, dim_zeilen);
	for(int i=0;i<dim_zeilen;++i)
		for(int j=i;j<dim_zeilen;++j){
			H.elem[i][j]=0;
			for(int k=0;k<dim_spalten;++k) H.elem[i][j]+=elem[i][k]*elem[j][k];
			if(i!=j) H.elem[j][i]=H.elem[i][j];
		}
	return H;
}

// schnellerer algorithmus fuer M.transpose()*M;
matrix matrix::transp_mult(){
	matrix H(dim_spalten, dim_spalten);
	for(int i=0;i<dim_spalten;++i)
		for(int j=i;j<dim_spalten;++j){
			H.elem[i][j]=0;
			for(int k=0;k<dim_zeilen;++k) H.elem[i][j]+=elem[k][i]*elem[k][j];
			if(i!=j) H.elem[j][i]=H.elem[i][j];
		}
	return H;
}

// -m : unitaerer matrix-operator; negiere komponentenweise
matrix matrix::operator-(){
	matrix H(dim_zeilen, dim_spalten);
	for(int i=0;i<dim_zeilen;++i)
		for(int j=0;j<dim_spalten;++j)
			H.elem[i][j]=elem[i][j]*(-1);
	return H;
}

// m*m = m
matrix matrix::operator*(const matrix& M){
	matrix H(dim_zeilen, M.dim_spalten);
	if(dim_spalten==M.dim_zeilen){
		for(int i=0;i<dim_zeilen;++i)
			for(int j=0;j<M.dim_spalten;++j){
				H.elem[i][j]=0;
				for(int k=0;k<dim_spalten;++k) H.elem[i][j]+=elem[i][k]*M.elem[k][j];
			}
	}else{
		cout <<"error: bei A*B muss gelten: spalten(A)==zeilen(B)"<<endl;
		H.zero();
	}
	return H;
}

// m*s = m (komponentenweise)
matrix matrix::operator*(double d){
	matrix H(*this);
	for(int i=0;i<dim_zeilen;++i)
		for(int j=0;j<dim_spalten;++j)
			H.elem[i][j]*=d;
	return H;
}

// m*v = m (muss u.u. manuell gecastet werden
matrix matrix::operator*(const vektor& v){
	int zeilen=1;
	int spalten=1;
	if(v.transp==0) zeilen=v.dim;
	else{
		if(v.transp==1) spalten=v.dim;
		else cout <<"error: aahhhhh!"<<endl;
	}
	matrix H(zeilen, spalten);
	int k=-1;
	for(int i=0;i<zeilen;++i)
		for(int j=0;j<spalten;++j)
			H.elem[i][j]=v.elem[++k];
	H=*this*H;
	return H;
}

// berechne teilweise v*m oder m*v
vektor matrix::multiply_sub(const vektor& v, int v_begin, int A_row_begin, int A_col_begin, int A_row_end, int A_col_end, char leftright){
	vektor r(v.dim);
	int dist_row=A_row_end-A_row_begin;
	int dist_col=A_col_end-A_col_begin;
	int error=(dist_row<0 || dist_col<0 || A_row_end>=dim_zeilen || A_col_end>=dim_spalten);
	r.transp=(leftright=='l');
	if(error || (r.transp && dist_row+v_begin>=v.dim) || (!r.transp && dist_col+v_begin>=v.dim))
		cout <<"error in matrix::multiply_sub(const vektor&, int, int, int, int, char): zugriff ausserhalb der dimensionen oder so"<<endl;
	else{
		if(leftright=='r'){ // m*v
			for(int i=0;i<=dist_row;++i){
				for(int j=0;j<=dist_col;++j){
					r.elem[i+v_begin]+=elem[i+A_row_begin][j+A_col_begin]*v.elem[j+v_begin];
				}
			}
		}else if(leftright=='l'){ // v*m
			for(int i=0;i<=dist_col;++i){
				for(int j=0;j<=dist_row;++j){
					r.elem[i+v_begin]+=v.elem[j+v_begin]*elem[j+A_row_begin][i+A_col_begin];
				}
			}
		}else cout <<"error in matrix::multiply_sub(const vektor&, int, int, int, int, char): char!='l' && char!='r'"<<endl;
	}
	return r;
}

// s*m = m
matrix operator*(double d, const matrix& M){
	matrix H(M);
	for(int i=0;i<H.dim_zeilen;++i)
		for(int j=0;j<H.dim_spalten;++j)
			H.elem[i][j]*=d;
	return H;
}

// v*m = m (muss manuell gecastet werden)
matrix operator*(vektor& v, const matrix& M){
	int zeilen=1;
	int spalten=1;
	if(v.transp==0) zeilen=v.get_dim();
	else{
		if(v.transp==1) spalten=v.get_dim();
		else cout <<"error: aahhhhh!"<<endl;
	}
	matrix H(zeilen, spalten);
	int k=-1;
	for(int i=0;i<zeilen;++i)
		for(int j=0;j<spalten;++j)
			H.elem[i][j]=v.get(++k);
	H*=M;
	return H;
}

// v*v = m (muss u.u. manuell gecastet werden mit m2s) (tensor-/skalarprodukt)
matrix operator*(vektor& v1, vektor& v2){
	int zeilen_1=1;
	int spalten_1=1;
	int zeilen_2=1;
	int spalten_2=1;
	if(v1.transp==0) zeilen_1=v1.get_dim();
	else if(v1.transp==1) spalten_1=v1.get_dim();
	if(v2.transp==0) zeilen_2=v2.get_dim();
	else if(v2.transp==1) spalten_2=v2.get_dim();
	matrix H1(zeilen_1, spalten_1);
	matrix H2(zeilen_2, spalten_2);
	matrix H(zeilen_1, spalten_2);
	int k=-1;
	for(int i=0;i<zeilen_1;++i)
		for(int j=0;j<spalten_1;++j)
			H1.elem[i][j]=v1.get(++k);
	k=-1;
	for(int i=0;i<zeilen_2;++i)
		for(int j=0;j<spalten_2;++j)
			H2.elem[i][j]=v2.get(++k);
	H=H1*H2;
	H1.deallocate();
	H2.deallocate();
	return H;
}

// m*=m = m
matrix& matrix::operator*=(const matrix& M){
	*this=*this*M;
	return *this;
}

// m*=s = m
matrix& matrix::operator*=(double d){
	for(int i=0;i<dim_zeilen;++i)
		for(int j=0;j<dim_spalten;++j)
			elem[i][j]*=d;
	return *this;
}

// m+m = m
matrix matrix::operator+(const matrix& M){
	matrix H(dim_zeilen, dim_spalten);
	if(dim_zeilen==M.dim_zeilen && dim_spalten==M.dim_spalten){
		for(int i=0;i<dim_zeilen;++i)
			for(int j=0;j<dim_spalten;++j)
				H.elem[i][j]=elem[i][j]+M.elem[i][j];
	}else{
		cout <<"error: dimensionen bei addition muessen gleich sein."<<endl;
		H.zero();
	}
	return H;
}

// m+s*Id = m; addiere nur auf der diagonalen!
matrix matrix::operator+(double d){
	matrix H(*this);
	for(int i=0;i<min(dim_zeilen, dim_spalten);++i) H.elem[i][i]+=d;
	return H;
}

// m+=m = m
matrix& matrix::operator+=(const matrix& M){
	if(dim_zeilen==M.dim_zeilen && dim_spalten==M.dim_spalten){
		for(int i=0;i<dim_zeilen;++i)
			for(int j=0;j<dim_spalten;++j)
				elem[i][j]+=M.elem[i][j];
	}else cout <<"error: dimensionen bei addition muessen gleich sein."<<endl;
	return *this;
}

// m+=s := m+=s*Id = m; addiere nur auf der diagonalen!
matrix& matrix::operator+=(double d){
	for(int i=0;i<min(dim_zeilen, dim_spalten);++i) elem[i][i]+=d;
	return *this;
}

// m-m = m
matrix matrix::operator-(const matrix& M){
	matrix H(dim_zeilen, dim_spalten);
	if(dim_zeilen==M.dim_zeilen && dim_spalten==M.dim_spalten){
		for(int i=0;i<dim_zeilen;++i)
			for(int j=0;j<dim_spalten;++j)
				H.elem[i][j]=elem[i][j]-M.elem[i][j];
	}else{
		cout <<"error: dimensionen bei addition muessen gleich sein."<<endl;
		H.zero();
	}
	return H;
}

// m-s*Id = m; ziehe nur auf der diagonalen ab!
matrix matrix::operator-(double d){
	matrix H(*this);
	for(int i=0;i<min(dim_zeilen, dim_spalten);++i) H.elem[i][i]-=d;
	return H;
}

// m-=m = m
matrix& matrix::operator-=(const matrix& M){
	if(dim_zeilen==M.dim_zeilen && dim_spalten==M.dim_spalten){
		for(int i=0;i<dim_zeilen;++i)
			for(int j=0;j<dim_spalten;++j)
				elem[i][j]-=M.elem[i][j];
	}else cout <<"error: dimensionen bei addition muessen gleich sein."<<endl;
	return *this;
}

// m-=s := m-=s*Id = m;  ziehe nur auf der diagonalen ab!
matrix& matrix::operator-=(double d){
	for(int i=0;i<min(dim_zeilen, dim_spalten);++i) elem[i][i]-=d;
	return *this;
}

// gib komponenten auf dem bildschirm aus
void matrix::ausgabe(int mode){
	// zeilenweises ausgeben der matrix
	if(mode==0){
		for(int i=0;i<dim_zeilen;++i){
			cout <<"[";
			for(int j=0;j<dim_spalten;++j){
				if(elem[i][j]==-0) elem[i][j]=0; // so ein kaes!
				if(j>0) cout <<" ";
				if(elem[i][j]>=0) cout <<" ";
				if(elem[i][j]!=0) printf("%5.4f", elem[i][j]);
				else cout <<"   0  ";
			}
			cout <<"]"<<endl;
		}
	}	else if(mode==1){
		ausgabe_format(1); // aus vektor.h
		for(int i=0;i<dim_zeilen;++i){
			cout <<"[";
			for(int j=0;j<dim_spalten;++j) cout <<setw(11)<<elem[i][j]<<" ";
			cout <<"\b]"<<endl;
		}
	}
}

ostream& matrix::operator<<(ostream& strm){
	ausgabe();
	return strm;
}

#endif