A powerful function parser

///////////////////////////////////////////////////////////////////////////////
//                                                                           //
// Dateiname: VEKTOREN.CPP                                                   //
//                                                                           //
// Autor:     Andreas J�ger, Friedrich-Schiller-Universit�t Jena             //
//                                                                           //
// System:    WIN_RWPM.EXE und WINLRWPM.EXE                                  //
//                                                                           //
// Beschreibung: Vektoren und Matrizen                                       //
//                                                                           //
// Hinweise:                                                                 //
//                                                                           //
///////////////////////////////////////////////////////////////////////////////
#include "stdafx.h"
#include <math.h>
#include <fstream.h>
#include "Collect.h"
#include "mathutil.h"
#include "Vektor.h"
#include "FktParser.h"
#include "Mathstr.h"

// Konstruktoren
CLinePoint::CLinePoint(int d) : Dim(d)
{
	hadr  = (Dim == 0) ? NULL : (new float[Dim]);
	Clear();
}

CLinePoint::CLinePoint(const CVektor& v)
{
	Dim   = v.Dim;
	hadr  = (Dim == 0) ? NULL : (new float[Dim]);
	for (int i=0; i<Dim; i++)
	{
		hadr[i] = (float)v.hadr[i];
	}
}

CLinePoint::CLinePoint(const CLinePoint& v)
{
	Dim   = v.Dim;
	hadr  = (Dim == 0) ? NULL : (new float[Dim]);
	for (int i=0; i<Dim; i++)
	{
		hadr[i] = v.hadr[i];
	}
}

CLinePoint::CLinePoint(long double t, const CVektor& x)
{
	Dim = x.Dim + 1;
	hadr  = new float[Dim];
	hadr[0] = (float)t;
	for (int i=1;i<=x.Dim;i++)
	{
		hadr[i] = (float)x.hadr[i-1];
	}
}

CLinePoint::CLinePoint(long double t, const CLinePoint& x)
{
	Dim = x.Dim + 1;
	hadr  = new float[Dim];
	hadr[0] = (float)t;
	for (int i=1;i<=x.Dim;i++)
	{
		hadr[i] = (float)x.hadr[i-1];
	}
}

void CLinePoint::Clear()
{
	if (hadr && (Dim > 0))
	{
		memset(hadr,0,Dim*sizeof(float));
	}
}

void CLinePoint::SetDim(int dim)
{
	if (dim == Dim) return;
	if (hadr) delete [] hadr;
	Dim = dim;
	if (Dim == 0) {hadr = NULL;return;}
	hadr = new float[Dim];
	Clear();
}

CLinePoint& CLinePoint::operator =(const CLinePoint  v) // ge�ndert
{
	if (&v != this)
	{
		if (v.Dim   != Dim  ) SetDim  (v.Dim);
		if (Dim > 0)
		{
			memmove(hadr,v.hadr,Dim*sizeof(float));
		}
	}
	return (*this);
}

float  CLinePoint::GetC(int t) const
{
	return hadr[t-1];
}

void CLinePoint::SetC(int t,float Val)
{
	hadr[t-1] = Val ;
}

CLinePoint::~CLinePoint()
{
	if (hadr) delete [] hadr;
}

void CLinePoint::Serialize(CArchive& ar)
{
	if (ar.IsStoring())
	{
		ar << Dim;
		ar.Write(hadr, Dim*sizeof(float));
	}
	else
	{
		ar >> Dim;
		delete[] hadr;
		hadr = new float[Dim];
		ar.Read(hadr, Dim*sizeof(float));
	}
}

//***********************CVektor****************************************

// Konstruktoren
CVektor::CVektor(int d, const char* name, bool t) : Name(name)
{
	Test  = 0;
	TDim  = 0;
	Dim   = d;
	sub   = 0;
	Status = (t) ? ZeilenVektor : SpaltenVektor;
	hadr  = (Dim == 0) ? NULL : (new long double[Dim]);
	if (Dim > 0 && hadr == NULL) Status = Invalid;
	Clear();
}

CVektor::CVektor(long double a)
{
	Test  = 0;
	TDim  = 0;
	Dim   = 1;
	sub   = 0;
	Status = SpaltenVektor;
	hadr  = new long double[Dim];
	(*hadr) = a;
}

CVektor::CVektor(long double a, long double b)
{
	Test  = 0;
	TDim  = 0;
	Dim   = 2;
	sub   = 0;
	Status = SpaltenVektor;
	hadr  = new long double[Dim];
	hadr[0] = a;
	hadr[1] = b;
}

CVektor::CVektor(long double a, long double b, long double c)
{
	Test  = 0;
	TDim  = 0;
	Dim   = 3;
	sub   = 0;
	Status = SpaltenVektor;
	hadr  = new long double[Dim];
	hadr[0] = a;
	hadr[1] = b;
	hadr[2] = c;
}

CVektor::CVektor(long double a, long double b, long double c, long double d)
{
	Test  = 0;
	TDim  = 0;
	Dim   = 4;
	sub   = 0;
	Status = SpaltenVektor;
	hadr  = new long double[Dim];
	hadr[0] = a;
	hadr[1] = b;
	hadr[2] = c;
	hadr[3] = d;
}

CVektor::CVektor(const CLinePoint  & v)
{
	Dim   = v.Dim;
	Test  = 0;
	TDim  = 0;
	sub = 0;
	hadr  = (Dim == 0) ? NULL : (new long double[Dim]);
	for (int i=0; i<Dim; i++)
	{
		hadr[i] = v.hadr[i];
	}
}

CVektor::CVektor(const CVektor& v) : Name(v.Name)
{
	Test  = 0;
	TDim  = 0;
	Dim   = v.Dim;
	sub   = v.sub;
	Status = v.Status;
	if (sub == 0)
	{
		hadr  = (Dim == 0) ? NULL : (new long double[Dim]);
		//if (Dim > 0 && hadr == NULL) Status = Invalid;
		//else
		//{
			/*_fmemcpy*/memmove(hadr,v.hadr,Dim*sizeof(long double));
		//}
	}
	else hadr = v.hadr;
}

// Achtung !!! Bei der Konvertierung einer Sub-Matrix in einen CVektor geht
// Achtung !!! deren Referenzeigenschaft verloren !!!

CVektor::CVektor(const CMatrix& M) : Name(M.Name)
{
	Test  = 0;
	TDim  = 0;
	if (M.xDim == 1)
	{
		Dim   = M.yDim;
		Status = ZeilenVektor;
		sub   = 0;
  		hadr  = new long double[Dim];
		//{
			/*_fmemcpy*/memmove(hadr,M.hadr[0],Dim*sizeof(long double));
		//}
	}
	if (M.yDim == 1)
	{
		Dim   = M.xDim;
		Status = SpaltenVektor;
		sub   = 0;
		hadr  = new long double[Dim];
		//{
			for (int i=1;i <= Dim;i++) hadr[i-1] = M.hadr[i-1][0];
		//}
	}
}

CVektor::CVektor(int d, short s, bool t, long double* z)
{
	Test  = 0;
	TDim  = 0;
	Dim   = d;
	sub   = s;
	Status = (t) ? ZeilenVektor : SpaltenVektor;
	if (sub == 0)
	{
		hadr = new long double[Dim];
		Clear();
	}
	else hadr = z;
}

CVektor::CVektor(long double t, const CVektor& x) 
{
	Test  = 0;
	TDim  = 0;
	Dim = x.Dim + 1;
	sub = 0;
	Status = ZeilenVektor;
	hadr  = new long double[Dim];
	hadr[0] = t;
	for (int i=1;i<=x.Dim;i++)
	{
		hadr[i] = x.hadr[i-1];
	}
}

CVektor::CVektor(long double t, const CLinePoint& x)
{
	Test  = 0;
	TDim  = 0;
	Dim = x.Dim + 1;
	sub = 0;
	Status = ZeilenVektor;
	hadr  = new long double[Dim];
	hadr[0] = t;
	for (int i=1;i<=x.Dim;i++)
	{
		hadr[i] = x.hadr[i-1];
	}
}


// Destruktor
CVektor::~CVektor()
{
	if ((sub == 0) && hadr) delete[] hadr;
}

CVektor& CVektor::EinheitsVektor(int k)
{
	Clear();
	if (k>0 && k<=Dim) hadr[k-1] = 1.0L;
	return *this;
}

long double CVektor::GetC(int t) const
{
	ASSERT(t>0 && t<=Dim);
	if (t<1 || t>Dim)
	{
		//MessageBox(NULL,"Error : CVektor::GetC(int t) ","Error in Dimensions !",MB_ICONHAND | MB_OK);
		return 0.0L;
	}
	return hadr[t-1];
}

long double& CVektor::GetC(int t)
{
	if (!Test)
	{
		if ((t < 1) || (t > Dim))
		{
			//MessageBox(NULL,"Error : CVektor::GetC(int t) ","Error in Dimensions !",MB_ICONHAND | MB_OK);
			//return 0.0L;//exit(1);
		}
		else 
			return hadr[t-1];
	}
	//else
	{
		if (t > TDim) TDim = t;
		return hadr[0];
	}
}

void CVektor::SetC(int t, long double Val)
{
	ASSERT(t>0 && t<=Dim);
	if (!Test)
	{
		hadr[t-1] = Val ;
	}
	else
	{
		if (t > TDim) TDim = t;
	}
}

long double&  CVektor::operator ()(int t)
{
	if (!Test)
	{
		if ((t < 1) || (t > Dim))
		{
			char buffer[255+1];
			sprintf(buffer," t = %d , Dim = %d ",t,Dim);
			MessageBox(NULL,"CVektor::operator()(int t) : Error in Dimensions ! ",buffer,MB_ICONHAND | MB_OK);
			CVektor y = *this;
			SetDim(t);
			for (int i=1; i<=y.Dim; i++) 
				SetC(i, y[i]);
			return hadr[t-1];
			//return 0.0L;//exit(1);
		}
		else 
			return hadr[t-1];
	}
	//else
	{
		if (t > TDim) TDim = t;
		return hadr[0];
	}
}

long double&  CVektor::operator [](int t)
{
	if (!Test)
	{
		if ((t < 1) || (t > Dim))
		{
			char buffer[255+1];
			sprintf(buffer," t = %d , Dim = %d ",t,Dim);
			::MessageBox(NULL, _T("CVektor::operator()(int t) : Error in Dimensions ! "), buffer, MB_ICONHAND | MB_OK);
			CVektor y = *this;
			SetDim(t);
			for (int i=1; i<=y.Dim; i++) 
				SetC(i, y[i]);
			return hadr[t-1];
			//return 0.0L;//exit(1);
		}
		else return hadr[t-1];
	}
	//else
	{
		if (t > TDim) TDim = t;
		return hadr[0];
	}
}

long double CVektor::operator ()(int t) const
{
	if (!Test)
	{
		if ((t < 1) || (t > Dim))
		{
			char buffer[255+1];
			sprintf(buffer," t = %d , Dim = %d ",t,Dim);
			MessageBox(NULL,"CVektor::operator()(int t) : Error in Dimensions ! ",buffer,MB_ICONHAND | MB_OK);
			return 0.0L;//exit(1);
		}
		//else
		return hadr[t-1];
	}
	else
	{
//		if (t > TDim) TDim = t;
		MessageBox(NULL,"Kann Dimension nicht testen. Konstantes Objekt.","CVektor",MB_ICONHAND | MB_OK);
		return hadr[0];
	}
}

long double CVektor::operator [](int t) const
{
	if (!Test)
	{
		if ((t < 1) || (t > Dim))
		{
			char buffer[255+1];
			sprintf(buffer," t = %d , Dim = %d ",t,Dim);
			MessageBox(NULL,"CVektor::operator()(int t) : Error in Dimensions ! ","CVektor",MB_ICONHAND | MB_OK);
			return 0.0L;//exit(1);
		}
		//else
		return hadr[t-1];
	}
	else
	{
//		if (t > TDim) TDim = t;
		MessageBox(NULL,"Kann Dimension nicht testen. Konstantes Objekt.","ec",MB_ICONHAND | MB_OK);
		return hadr[0];
	}
}

CMatrix CVektor::Diag() const
{
	CMatrix M(Dim, Dim);
	for (int i=1; i<=Dim; i++)
	{
		M(i, i) = GetC(i);
	}
	return M;
}

CVektor sin(const CVektor& x)
{
	CVektor v(x.Dim);
	for (int i=0; i<x.Dim; i++)
	{
		v.hadr[i] = sinl(x.hadr[i]);
	}
	return v;
}

CVektor cos(const CVektor& x)
{
	CVektor v(x.Dim);
	for (int i=0; i<x.Dim; i++)
	{
		v.hadr[i] = cosl(x.hadr[i]);
	}
	return v;
}

CVektor tan(const CVektor& x)
{
	CVektor v(x.Dim);
	for (int i=0; i<x.Dim; i++)
	{
		v.hadr[i] = tanl(x.hadr[i]);
	}
	return v;
}

CVektor cot(const CVektor& x)
{
	CVektor v(x.Dim);
	for (int i=0; i<x.Dim; i++)
	{
		v.hadr[i] = cosl(x.hadr[i])/sinl(x.hadr[i]);
	}
	return v;
}

CVektor abs(const CVektor& x)
{
	CVektor v(x.Dim);
	for (int i=0; i<x.Dim; i++)
	{
		v.hadr[i] = fabsl(x.hadr[i]);
	}
	return v;
}

CVektor UserFunction(const CVektor& x, const char* FunktionsTerm, const char* Variable)
{
	CString FktTerm(FunktionsTerm);
	CFunction<long double>* fkt = CFunction<long double>::Parse(FktTerm, Variable);
	CVektor v(x.Dim);

	if (fkt==NULL)
		return v;

	for (int i=0; i<x.Dim; i++)
	{
		v.hadr[i] = fkt->Execute(x.hadr[i]);
	}
	delete fkt;
	return v;
}

void CVektor::sin()
{
	for (int i=0; i<Dim; i++)
	{
		hadr[i] = ::sinl(hadr[i]);
	}
}

void CVektor::cos()
{
	for (int i=0; i<Dim; i++)
	{
		hadr[i] = ::cosl(hadr[i]);
	}
}

void CVektor::tan()
{
	for (int i=0; i<Dim; i++)
	{
		hadr[i] = ::tanl(hadr[i]);
	}
}

void CVektor::cot()
{
	for (int i=0; i<Dim; i++)
	{
		hadr[i] = ::cosl(hadr[i]) / ::sinl(hadr[i]);
	}
}

void CVektor::abs()
{
	for (int i=0; i<Dim; i++)
	{
		hadr[i] = fabsl(hadr[i]);
	}
}

void CVektor::UserFunction(const char* FunktionsTerm, const char* Variable)
{
	CString FktTerm(FunktionsTerm);
	CFunction<long double>* fkt = CFunction<long double>::Parse(FktTerm, Variable);

	if (fkt==NULL)
		return;

	for (int i=0; i<Dim; i++)
	{
		hadr[i] = fkt->Execute(hadr[i]);
	}
	delete fkt;
}

void CVektor::SetTrans(bool t)
{
	if (sub == 0)
	{
	  Status = (t) ? ZeilenVektor : SpaltenVektor;
	}
}

void CVektor::Transpose()
{
	if (sub == 0)
	{
		if (Status == ZeilenVektor)
		{
			Status = SpaltenVektor;
			return;
		}
		if (Status == SpaltenVektor)
		{
			Status = ZeilenVektor;
			return;
		}
	}
}

CVektor& CVektor::operator =(const CVektor& v) 
{
	Name = v.Name;
	if (&v != this)
	{
		if (v.Dim != Dim)
			SetDim(v.Dim);
		Status = v.Status;
		if (Dim > 0)
		{
			/*_fmemcpy*/memmove(hadr, v.hadr, Dim*sizeof(long double));
		}
	}
	return *this;
}

CVektor& CVektor::operator =(long double x) 
{
	SetDim(1);
	hadr[0] = x;
	return *this;
}

CVektor& CVektor::operator =(const CString& /*str*/) 
{
	//NotImplemented();
	return *this;
}

CVektor& CVektor::operator +=(const CVektor& v)
{
	for (int i=1; i<=Dim; i++)
		hadr[i-1] += v.hadr[i-1];
	return *this;
}

CVektor& CVektor::operator -=(const CVektor& v)
{
	for (int i=1; i<=Dim; i++)
		hadr[i-1] -= v.hadr[i-1];
	return (*this);
}

CVektor& CVektor::operator *=(long double a)
{
	for (int i=1; i<=Dim; i++)
		hadr[i-1] *= a;
	return (*this);
}

CVektor& CVektor::operator /=(long double a)
{
	for (int i=1; i<=Dim; i++)
		hadr[i-1] /= a;
	return (*this);
}

CVektor& CVektor::operator =(const CMatrix& M)
{
	Name = M.Name;
	if (M.xDim == 1)
	{
		if (Dim != M.yDim) 
			SetDim(M.yDim);
		if ((Status != ZeilenVektor ) && (Dim != 1 )) 
			SetTrans();
		/*_fmemcpy*/memmove(hadr, M.hadr[0], Dim*sizeof(long double));
	}
	if ((M.yDim == 1) && (M.xDim != 1))
	{
		if (Dim != M.xDim) 
			SetDim(M.xDim);
		if (Status != SpaltenVektor) 
			SetTrans(SpaltenVektor);
		for (int i=1; i<=Dim; i++) 
			hadr[i-1] = M.hadr[i-1][0];
	}
	return (*this);
}

CVektor CVektor::operator ()(int i,int j)
{
	return CVektor(j-i+1, 1, Status == ZeilenVektor, &hadr[i-1]);
}

int operator==(const CVektor& v, const CVektor& w)
{
	if (v.Dim != w.Dim)
		return 0;
	for (int i=1; i<=v.Dim; i++)
	{
		if (v[i] != w[i])
			return 0;
	}
	return 1;
}

int operator!=(const CVektor& v, const CVektor& w)
{
	return !(v == w);
}

CVektor CVektor::GetSubVektor(int i, int j) const
{
	ASSERT(i>0 && i<=j && j<=Dim);
	CVektor v(j-i+1);
	v.Status = Status;
	memcpy(v.hadr, &hadr[i-1], v.Dim * sizeof(long double));
	return v;
}

CVektor operator !(const CVektor& v)
{
	CVektor w(v.Dim);
	if (v.Status == CVektor::ZeilenVektor)
		w.Status = CVektor::SpaltenVektor; 
	else
		w.Status = CVektor::ZeilenVektor;
	if (v.Dim>0)
	{
		memmove(w.hadr, v.hadr, v.Dim*sizeof(long double));
	}
	return w;
}

CVektor operator +(const CVektor& v, const CVektor& w)
{
	CVektor s(v.Dim, CString(), v.Status == CVektor::ZeilenVektor);
	for (int i=1; i<=v.Dim; i++) 
		s[i] = v[i] + w[i];
	return s;
}

CVektor operator -(const CVektor& v, const CVektor& w)
{
	CVektor s(v.Dim, CString(), v.Status == CVektor::ZeilenVektor);
	for (int i=1; i<=v.Dim; i++) 
		s[i] = v[i] - w[i];
	return s;
}

CVektor operator *(const long double a, const CVektor& v)
{
	CVektor s(v.Dim, CString(), v.Status == CVektor::ZeilenVektor);
	for (int i=1; i<=v.Dim; i++) 
		s[i] = a*v[i];
	return s;
}

CVektor operator *(const CVektor& v, const long double a)
{
	CVektor s(v.Dim, CString(), v.Status == CVektor::ZeilenVektor);
	for (int i=1; i<=v.Dim; i++) 
		s[i] = a*v[i];
	return s;
}

CVektor operator /(const CVektor& v, const CVektor& w)
{
	CVektor s(v.Dim, CString(), v.Status == CVektor::ZeilenVektor);
	for (int i=1; i<=v.Dim; i++) 
		s[i] = v[i] / w[i];
	return s;
}

CVektor operator /(const CVektor& v, const long double a)
{
	CVektor s(v.Dim, CString(), v.Status == CVektor::ZeilenVektor);
	for (int i=1; i<=v.Dim; i++) 
		s[i] = v[i]/a;
	return s;
}

CMatrix operator *(const CVektor& v, const CVektor& w)
{
	CMatrix p(v.Dim, w.Dim);
	if ((v.Status == CVektor::SpaltenVektor) && (w.Status == CVektor::ZeilenVektor)) 
	{
		for (int i=1; i<=v.Dim; i++) 
			p(i,i,1,w.Dim) = v[i] * w;
	}
	if ((v.Status == CVektor::ZeilenVektor) && (w.Status == CVektor::SpaltenVektor)) 
	{
		long double s = 0.0L;
		for (int i=1; i<=v.Dim; i++) 
			s += v[i] * w[i];
		p = (CMatrix)s;
	}
	return p;
}

void CVektor::Serialize(CArchive& ar)
{
	if (ar.IsStoring())
	{
		ar << Name;
		ar << Dim;
		//ar << Name.GetLength();
		//ar.WriteString(Name);
		ar << Status;
		ar << sub;
		ar.Write(hadr, Dim*sizeof(long double));
	}
	else
	{
		ar >> Name;
		ar >> Dim;
		if (Dim < 0)
			Dim = 0;
		//int i;
		//ar >> i;
		//ar.Read();
		int status;
		ar >> status;
		Status = (CStatus)status;
		ar >> sub;
		TDim = Dim;
		Test = 0;
		delete[] hadr;
		hadr = new long double[Dim];
		if (ar.Read(hadr, Dim*sizeof(long double)) < Dim *sizeof(long double))
		{
			memset(hadr, 0, Dim*sizeof(long double));
		}
	}
}

// File-Format:
// t = [
//    1, 2, 89.873, 24
// ]
istream& Read(istream& is, CVektor& v)
{
	static char temp[1000];
	is.getline(temp, sizeof(temp), '['); // bis '[' �berlesen
//	::MessageBox(NULL, temp, "Vorspann", MB_OK);
	is.getline(temp, sizeof(temp), ']');
//	::MessageBox(NULL, temp, "Eigentlicher Vektor", MB_OK);
	CString s(temp);
	CString t;
	size_t p1 = 0;
	size_t p2 = 0;
	CCollection<CMathString> testcoll(_T(""), 10, 10);
	/*
	while (p2 != -1)
	{
		p2 = s.Find(",", p1);
		if (p2 == -1)
		{
			t = CString(s, p1);
		}
		else
		{
			t = s.Mid(p1, p2-p1);
		}
		testcoll.Insert(new CMathstring(t));
//		::MessageBox(NULL, t.c_str(), "Komponente", MB_OK);
		p1 = p2+1;
	}
	*/
	v.SetDim(testcoll.Count());
	for (int i=1; i<=v.Dim; i++)
	{
		v[i] = testcoll.At(i-1)->Value();
	}
//	::MessageBox(NULL, v, "Vektor gelesen", MB_OK);
	return is;
}

istream& Read(istream& is, CMatrix& M)
{
	static char temp[1000];
	is.getline(temp, sizeof(temp), '['); // bis '[' �berlesen
//	::MessageBox(NULL, temp, "Vorspann", MB_OK);
	is.getline(temp, sizeof(temp), ']');
//	::MessageBox(NULL, temp, "Eigentliche Matrix", MB_OK);
	CString s(temp);
	CString t;
	int i, j;
	size_t p1 = 0;
	size_t p2 = 0;
	CCollection<CMathString> testcoll(_T(""), 10, 10);
	/*
	while (p2 != -1)
	{
		p2 = s.Find(";", p1); // Zeilen-Separator
		if (p2 == -1)
		{
			t = CString(s, p1);
		}
		else
		{
			t = s.Mid(p1, p2-p1);
		}
		testcoll.Insert(new CMatrixhstring(t));
//		::MessageBox(NULL, t.c_str(), "Zeile", MB_OK);
		p1 = p2+1;
	}
	*/
	// Die Anzahl der Zeilen ist testcoll.Count()
	int Zeilen = testcoll.Count();
	int Spalten = 1;
	for (i=1; i<=testcoll.Count(); i++)
	{
		Spalten = 1;
		size_t q1 = 0;
		size_t q2 = 0;
		j=1;
		/*
		while (q2 != -1)
		{
			q2 = testcoll.At(i-1)->Find(",", q1); // Zeilen-Separator
			if (q2 != NPOS)
			{
				Spalten = max(Spalten, j+1);
			}
			q1 = q2+1;
			j++;
		}
		*/
	}
	M.SetDim(Zeilen, Spalten);
	for (i=1; i<=Zeilen; i++)
	{
		size_t r1 = 0;
		size_t r2 = 0;
		j = 1;
		/*
		while (r2 != -1)
		{
			r2 = testcoll.At(i-1)->find(",", r1);
			if (r2 == -1)
			{
				t = CString(*(testcoll.At(i-1)), r1);
			}
			else
			{
				t = CString(*(testcoll.At(i-1)), r1, r2-r1);
			}
			m(i, j) = CMathString(t).value();
//			::MessageBox(NULL, t.c_str(), "Komponente", MB_OK);
			r1 = r2+1;
			j++;
		}
		*/
	}
//	::MessageBox(NULL, m, "Matrix gelesen", MB_OK);
	return is;
}

ostream& operator <<(ostream& os, CVektor& v)
{
	os.setf(os.scientific,os.floatfield);os.precision(18);
	if (os.bad()) throw;//(xmsg("Write failed"));
	os << v.Dim; if (os.bad()) throw;//(xmsg("Write failed"));
	os << "\n"; if (os.bad()) throw;//(xmsg("Write failed"));
	os << (int)v.Status; if (os.bad()) throw;//(xmsg("Write failed"));
	os << "\n"; if (os.bad()) throw;//(xmsg("Write failed"));
	for(int i=1;i <= v.Dim;i++)
	{
		os << v(i); if (os.bad()) throw;//(xmsg("Write failed"));
		os << " " ; if (os.bad()) throw;//(xmsg("Write failed"));
	}
	return os;
}

istream& operator >>(istream& is, CVektor& v)
{
//  is.setf(ios::scientific,ios::floatfield);is.precision(18);
	int _Dim;

	is >> _Dim;
	v.SetDim(_Dim);
	int t;
	is >> t;
	v.Status = (CVektor::CStatus)t;
	for (int i=1; i<=_Dim; i++)
	{
		is >> v[i];
	}
	return is;
}

ostream& CVektor::SaveBinary(ostream& os)
{
	os.write((const char*)(&Dim), sizeof(Dim));
	os.write((const char*)(&Status), sizeof(Status));
	os.write((const char*)(&sub), sizeof(sub));
	int NameLen = Name.GetLength();
	os.write((const char*)(&NameLen), sizeof(NameLen));
	os.write((const char*)(&Name), NameLen);
	os.write((const char*)(&Test), sizeof(Test));
	os.write((const char*)(&TDim), sizeof(TDim));
	os.write((const char*)hadr, Dim*sizeof(long double));
	return os;
}

istream& CVektor::LoadBinary(istream& is)
{
	is.read((char*)(&Dim), sizeof(Dim));
	is.read((char*)(&Status), sizeof(Status));
	is.read((char*)(&sub), sizeof(sub));
	int NameLen;
	is.read((char*)(&NameLen), sizeof(NameLen));
	char* buffer = new char[NameLen+1];
	memset(buffer, 0, NameLen+1);
	is.read(buffer, NameLen);
	Name = buffer;
	is.read((char*)(&Test), sizeof(Test));
	is.read((char*)(&TDim), sizeof(TDim));
	delete [] hadr;
	hadr = new long double[Dim];
	is.read((char*)hadr, Dim*sizeof(long double));
	return is;
}

ostream& operator <<(ostream& os, CLinePoint& v)
{
	os.setf(os.scientific,os.floatfield);os.precision(18);
	if (os.bad()) throw;//(xmsg("Write failed"));
	os << v.Dim; if (os.bad()) throw;//(xmsg("Write failed"));
	os << "\n"; if (os.bad()) throw;//(xmsg("Write failed"));
	for(int i=1;i <= v.Dim;i++)
	{
		os << v[i]; if (os.bad()) throw;//(xmsg("Write failed"));
		os << " " ; if (os.bad()) throw;//(xmsg("Write failed"));
	}
	return os;
}

istream& operator >>(istream& is, CLinePoint& v)
{
//  is.setf(ios::scientific,ios::floatfield);is.precision(18);
	int _Dim;

	is >> _Dim;
	v.SetDim(_Dim);
	for (int i=1; i<=_Dim; i++)
	{
		is >> v[i];
	}
	return is;
}

void CVektor::Clear()
{
	if (hadr && (Dim > 0))
	{
		for (int i=0; i<Dim; i++) hadr[i] = 0.0L;
		//memset(hadr,0,Dim*sizeof(long double));
	}
}

void CVektor::SetDim(int dim)
{
	if (dim == Dim) return;
	if (sub == 0)
	{
		if (hadr) delete [] hadr;
		Dim = dim;
		if (Dim == 0)
		{
			hadr = NULL;
			return;
		}
		if ((hadr = new long double[Dim]) == NULL) Status = Invalid;
		Clear();
	}
}

void CVektor::ForceDimension(int dim)
{
	CVektor v = *this;
	SetDim(dim);
	for (int i=1; i<=Dim && i<=v.Dim; i++)
	{
		hadr[i-1] = v.hadr[i-1];
	}
}

long double CVektor::Norm(int n) const
{
	long double  n1;
	switch (n) 
	{
		case 1:
		{
			n1=0.0L;
			for (int i=1; i<=Dim; i++)
			{
				n1+=fabsl(hadr[i-1]);
			}
			break;
		}
		case 2:
		{
			n1=0.0L;
			for (int i=1; i<=Dim; i++)
			{
				n1+=sqrl(hadr[i-1]);
			}
			n1=sqrtl(n1);
			break;
		}
		case 0:
		case infty:
		{
			n1=((Dim>0) ? fabsl(hadr[0]) : 0.0);;
			for (int i=2; i<=Dim; i++)
			{
				if (fabsl(hadr[i-1]) > n1)
				{
					n1=fabsl(hadr[i-1]);
				}
			}
			break;
		}
		default:
		{
			n1=0.0L;
			for (int i=1; i<=Dim; i++)
				n1 += powl(hadr[i-1], n);
			n1=powl(n1, 1.0L/n);
		}
	}
	return n1;
}

long double Norm(const CVektor& v, int n)
{
	long double  n1;
	switch (n) 
	{
		case 1:
		{
			n1=0.0L;
			for (int i=1; i<=v.Dim; i++)
			{
				n1+=fabsl(v[i]);
			}
			break;
		}
		case 2:
		{
			n1=0.0L;
			for (int i=1; i<=v.Dim; i++)
			{
				n1+=sqrl(v[i]);
			}
			n1=sqrtl(n1);
			break;
		}
		case 0:
		case infty:
		{
			n1=((v.Dim>0) ? fabsl(v[1]) : 0.0);
			for (int i=2; i<=v.Dim; i++)
			{
				if (fabsl(v[i]) > n1)
				{
					n1=fabsl(v[i]);
				}
			}
			break;
		}
		default:
		{
			n1=0.0L;
			for (int i=1; i<=v.Dim; i++)
				n1 += powl(v[i], n);
			n1=powl(n1, 1.0L/n);
		}
	}
	return n1;
}

long double Norm(const CVektor& v1, const CVektor& v2, int n)
{
	long double  n1;
	int Dim = min(v1.Dim, v2.Dim);
	switch (n)
	{
	case 1:
	{
		n1=0.0L;
		for (int i=1;i <= Dim;i++)
		{
			n1+=fabsl(v1[i]-v2[i]);
		}
		break;
	}
	case 2:
	{
		n1=0.0L;
		for (int i=1;i <= Dim;i++)
		{
			n1+=sqrl(v1[i]-v2[i]);
		}
		n1=sqrtl(n1);
		break;
	}
	case 0:
	case infty:
	{
		n1=((Dim>0) ? fabsl(v1[1]-v2[1]) : 0.0);
		for (int i=2;i <= Dim;i++)
		{
			if (fabsl(v1[i]-v2[i]) > n1)
			{
				n1=fabsl(v1[i]-v2[i]);
			}
		}
		break;
	}
	default:
	{
		n1=0.0L;
		for (int i=1; i<=Dim; i++)
			n1 += powl(v1[i]-v2[i], n);
		n1=powl(n1, 1.0L/n);
	}
	}
	return n1;
}

long double CVektor::SqrNorm2() const   // Quadrat der euklid. Norm
{
	long double  n1=0.0L;
	for (int i=1;i <= Dim;i++)
	{
		n1+=sqrl(hadr[i-1]);
	}
	return n1;
}

long double SqrNorm2(const CVektor& v)    // Quadrat der euklid. Norm
{
	long double  n1=0.0L;
	for (int i=1;i <= v.Dim;i++)
	{
		n1+=sqrl(v(i));
	}
	return n1;
}

long double SqrNorm2(const CVektor& v1, const CVektor& v2)    // Quadrat der euklid. Norm
{
	long double  n1=0.0L;
	int Dim = min(v1.Dim, v2.Dim);
	for (int i=1;i <= Dim;i++)
	{
		n1+=sqrl(v1[i]-v2[i]);
	}
	return n1;
}

void CVektor::EinsVektor()
{
	for (int i=0; i<Dim; i++)
		hadr[i] = 1.0;
}

CVektor NullVektor(int dim)
{
	return CVektor(dim);
}

CVektor EinheitsVektor(int dim, int k)
{
	CVektor v(dim);
	v.hadr[k-1] = 1.0;
	return v;
}

CVektor EinsVektor(int dim)
{
	CVektor v(dim);
	for (int i=0; i<v.Dim; i++)
		v.hadr[i] = 1.0;
	return v;
}

CMatrix Diag(const CVektor& v)
{
	CMatrix M(v.Dim, v.Dim);
	for (int i=0; i<v.Dim; i++)
		M.hadr[i][i] = v.hadr[i];
	return M;
}

long double Skalarprodukt(const CVektor& v, const CVektor& w)
{
	ASSERT(v.Dim==w.Dim);
	long double sp = 0.0;
	for (int i=1; i<=v.Dim && i<=w.Dim; i++)
	{
		sp += v[i]*w[i];
	}
	return sp;
}

CVektor Kreuzprodukt(const CVektor& v, const CVektor& w)
{
	ASSERT(v.Dim==w.Dim);
	ASSERT(v.Dim==2 || v.Dim==3);
	CVektor kp(3);
	if (v.Dim==3 && w.Dim==3)
	{
		kp[1] = v[2]*w[3] - v[3]*w[2];
		kp[2] = v[3]*w[1] - v[1]*w[3];
		kp[3] = v[1]*w[2] - v[2]*w[1];
	}
	if (v.Dim==2 && w.Dim==2)
	{
		kp[1] = 0.0;
		kp[2] = 0.0;
		kp[3] = v[1]*w[2] - v[2]*w[1];
	}
	return kp;
}

long double SpatProdukt(const CVektor& u, const CVektor& v, const CVektor& w)
{
	ASSERT(u.Dim==3 && v.Dim==3 && w.Dim==3);
	long double sp = 0.0;
	if (u.Dim==3 && v.Dim==3 && w.Dim==3)
	{
		sp += (u[2]*v[3] - u[3]*v[2]) * w[1];
		sp += (u[3]*v[1] - u[1]*v[3]) * w[2];
		sp += (u[1]*v[2] - u[2]*v[1]) * w[3];
	}
	return sp;
}

void CVektor::Draw(CDC* /*dc*/, CRect& /*rect*/)
{
}

//**********************CMatrix*************************************

CMatrix::CMatrix(int d1,int d2, const char* name) : Name(name)
{
	Status = Normal;
	perm = NULL;
	xDim = d1;
	yDim = d2;
	sub  = 0;
	RSN  = 0;
	hadr = new long double*[xDim];
	for (int i=1;i <= xDim;i++) 
		hadr[i-1] = new long double[yDim];
	Clear();
}

/*CMatrix::CMatrix(int d1)
{
	Status = Normal;
	perm = NULL;
	xDim = d1;
	yDim = d1;
	sub  = 0;
	RSN  = 0;
	hadr = new long double*[xDim];
	for (int i=1;i <= xDim;i++) hadr[i-1] = new long double[yDim];
	Clear();
}*/

CMatrix::CMatrix()
{
	Status = Normal;
	perm = NULL;
	xDim  = 1;
	yDim  = 1;
	sub   = 0;
	RSN   = 0;
	hadr  = new long double*[1];
	*hadr = new long double[1];
	Clear();
}

CMatrix::CMatrix(long double a11)
{
	Status = Normal;
	perm   = NULL;
	xDim   = 1;
	yDim   = 1;
	sub    = 0;
	RSN    = 0;
	hadr   = new long double*[1];
	*hadr  = new long double[1];
	**hadr = a11;
}

CMatrix::CMatrix(long double a11, long double a12, long double a21, long double a22)
{
	Status  = Normal;
	perm    = NULL;
	xDim    = 2;
	yDim    = 2;
	sub     = 0;
	RSN     = 0;
	hadr    = new long double*[2];
	hadr[0] = new long double[2];
	hadr[1] = new long double[2];
	hadr[0][0] = a11; hadr[0][1] = a12;
	hadr[1][0] = a21; hadr[1][1] = a22;
}

CMatrix::CMatrix(long double a11, long double a12, long double a13, long double a21, long double a22, long double a23, long double a31, long double a32, long double a33)
{
	Status  = Normal;
	perm    = NULL;
	xDim    = 3;
	yDim    = 3;
	sub     = 0;
	RSN     = 0;
	hadr    = new long double*[3];
	hadr[0] = new long double[3];
	hadr[1] = new long double[3];
	hadr[2] = new long double[3];
	hadr[0][0] = a11; hadr[0][1] = a12; hadr[0][2] = a13;
	hadr[1][0] = a21; hadr[1][1] = a22; hadr[1][2] = a23;
	hadr[2][0] = a31; hadr[2][1] = a32; hadr[2][2] = a33;
}

CMatrix::CMatrix(long double a11, long double a12, long double a13, long double a14, long double a21, long double a22, long double a23, long double a24, long double a31, long double a32, long double a33, long double a34, long double a41, long double a42, long double a43, long double a44)
{
	Status  = Normal;
	perm    = NULL;
	xDim    = 4;
	yDim    = 4;
	sub     = 0;
	RSN     = 0;
	hadr    = new long double*[4];
	hadr[0] = new long double[4];
	hadr[1] = new long double[4];
	hadr[2] = new long double[4];
	hadr[3] = new long double[4];
	hadr[0][0] = a11; hadr[0][1] = a12; hadr[0][2] = a13; hadr[0][3] = a14;
	hadr[1][0] = a21; hadr[1][1] = a22; hadr[1][2] = a23; hadr[1][3] = a24;
	hadr[2][0] = a31; hadr[2][1] = a32; hadr[2][2] = a33; hadr[2][3] = a34;
	hadr[3][0] = a41; hadr[3][1] = a42; hadr[3][2] = a43; hadr[3][3] = a44;
}

CMatrix::CMatrix(int d1,int d2,int s)
{
	Status = Normal;
	perm = NULL;
	xDim = d1;
	yDim = d2;
	sub  = s;
	RSN  = 0;
	hadr = new long double*[xDim];
	if (sub == 0)
	{
		for (int i=1;i <= xDim;i++) hadr[i-1] = new long double[yDim];
		Clear();
	}
	else
	{
		for (int i=1;i <= xDim;i++) hadr[i-1] = NULL;
	}
}

CMatrix::CMatrix (const CMatrix& M) : Name(M.Name)
{
	Status = M.Status;;
	perm = NULL;
	xDim = M.xDim;
	yDim = M.yDim;
	sub  = M.sub;
	RSN  = 0;
	hadr = new long double*[xDim];
	if (sub == 0)
	{
		for (int i=1;i <= xDim;i++) hadr[i-1] = new long double[yDim];
		for (i=1;i <= xDim ;i++)
		{
	  		memmove(hadr[i-1],M.hadr[i-1],yDim*sizeof(long double));
		}
	}
	else
	{
		for (int i=1;i <= xDim;i++) hadr[i-1] = M.hadr[i-1];
	}
}

// Achtung !!! Bei der Konvertierung eines Sub-Vektors in eine Matrix geht
// Achtung !!! dessen Referenzeigenschaft verloren !!!

CMatrix::CMatrix(const CVektor& v) : Name(v.Name)
{
	Status = Normal;
	perm = NULL;
	if (v.Status == CVektor::ZeilenVektor)
	{
		xDim    = 1;
		yDim    = v.Dim;
		sub     = 0;
		RSN     = 0;
		hadr    = new long double*[xDim];
		hadr[0] = new long double[yDim];
		memmove(hadr[0],v.hadr,yDim*sizeof(long double));
	}
	if (v.Status == CVektor::SpaltenVektor)
	{
		xDim = v.Dim;
		yDim = 1;
		sub  = 0;
		RSN  = 0;
		hadr = new long double*[xDim];
		for (int i=1;i <= xDim;i++) 
			hadr[i-1] = new long double[yDim];
		for (i=1;i <= xDim;i++) 
			hadr[i-1][0] = v.hadr[i-1];
	}
}

CMatrix::CMatrix(const CVektor& v1, const CVektor& v2, bool zeile)
{
	Status = Normal;
	perm = NULL;
	if (zeile)
	{
		xDim    = 2;
		yDim    = v1.Dim;
		sub     = 0;
		RSN     = 0;
		hadr    = new long double*[xDim];
		hadr[0] = new long double[yDim];
		hadr[1] = new long double[yDim];
		memmove(hadr[0],v1.hadr,yDim*sizeof(long double));
		memmove(hadr[1],v2.hadr,yDim*sizeof(long double));
	}
	else
	{
		xDim = v1.Dim;
		yDim = 2;
		sub  = 0;
		RSN  = 0;
		hadr = new long double*[xDim];
		for (int i=1;i <= xDim;i++) 
			hadr[i-1] = new long double[yDim];
		for (i=1;i <= xDim;i++) 
		{
			hadr[i-1][0] = v1.hadr[i-1];
			hadr[i-1][1] = v2.hadr[i-1];
		}
	}
}

CMatrix::CMatrix(const CVektor& v1, const CVektor& v2, const CVektor& v3, bool zeile)
{
	Status = Normal;
	perm = NULL;
	if (zeile)
	{
		xDim    = 3;
		yDim    = v1.Dim;
		sub     = 0;
		RSN     = 0;
		hadr    = new long double*[xDim];
		hadr[0] = new long double[yDim];
		hadr[1] = new long double[yDim];
		hadr[2] = new long double[yDim];
		memmove(hadr[0],v1.hadr,yDim*sizeof(long double));
		memmove(hadr[1],v2.hadr,yDim*sizeof(long double));
		memmove(hadr[2],v3.hadr,yDim*sizeof(long double));
	}
	else
	{
		xDim = v1.Dim;
		yDim = 3;
		sub  = 0;
		RSN  = 0;
		hadr = new long double*[xDim];
		for (int i=1;i <= xDim;i++) 
			hadr[i-1] = new long double[yDim];
		for (i=1;i <= xDim;i++) 
		{
			hadr[i-1][0] = v1.hadr[i-1];
			hadr[i-1][1] = v2.hadr[i-1];
			hadr[i-1][2] = v3.hadr[i-1];
		}
	}
}

CMatrix::CMatrix(const CVektor& v1, const CVektor& v2, const CVektor& v3, const CVektor& v4, bool zeile)
{
	Status = Normal;
	perm = NULL;
	if (zeile)
	{
		xDim    = 4;
		yDim    = v1.Dim;
		sub     = 0;
		RSN     = 0;
		hadr    = new long double*[xDim];
		hadr[0] = new long double[yDim];
		hadr[1] = new long double[yDim];
		hadr[2] = new long double[yDim];
		hadr[3] = new long double[yDim];
		memmove(hadr[0],v1.hadr,yDim*sizeof(long double));
		memmove(hadr[1],v2.hadr,yDim*sizeof(long double));
		memmove(hadr[2],v3.hadr,yDim*sizeof(long double));
		memmove(hadr[3],v4.hadr,yDim*sizeof(long double));
	}
	else
	{
		xDim = v1.Dim;
		yDim = 4;
		sub  = 0;
		RSN  = 0;
		hadr = new long double*[xDim];
		for (int i=1;i <= xDim;i++) 
			hadr[i-1] = new long double[yDim];
		for (i=1;i <= xDim;i++) 
		{
			hadr[i-1][0] = v1.hadr[i-1];
			hadr[i-1][1] = v2.hadr[i-1];
			hadr[i-1][2] = v3.hadr[i-1];
			hadr[i-1][3] = v4.hadr[i-1];
		}
	}
}

/*
CMatrix::CMatrix(CVektorList  &  L)
{
	xDim = L.CN().Dim;
	yDim = L.Nr_of_Nodes;
	sub  = 0;
	RSN  = 0;
	hadr = new long double*[xDim];
	for (int i=1;i <= xDim;i++) 
	hadr[i-1] = new long double[yDim];
	Clear();

	CVektorListIterator j(L);i=1;
	while (j) { (*this)(1,xDim,i,i) = j.Current();j++;i++;}
}
*/

CMatrix::~CMatrix()
{
	if (sub == 0)
	{
		for (int i=1;i <= xDim;i++) delete [] hadr[i-1];
		delete [] hadr;
	}
	else delete [] hadr;
	if (perm) delete [] perm;
}

long double& CMatrix::GetC(int i,int j)
{
	if ((i<1) || (i>xDim) || (j<1) || (j>yDim))
	{
		return hadr[0][0]; // Fehler
	}
	return hadr[i-1][j-1];
}

long double  CMatrix::GetC(int i,int j) const
{
	if ((i<1) || (i>xDim) || (j<1) || (j>yDim))
	{
		return 0.0L;
	}
	return hadr[i-1][j-1];
}

void CMatrix::SetC(int i,int j,long double Val)
{
	if ((i<1) || (i>xDim) || (j<1) || (j>yDim))
	{
  		return;
	}
	hadr[i-1][j-1] = Val ;
}

long double& CMatrix::operator ()(int i,int j)
{
	if (i<1 || i>xDim || j<1 || j>yDim)
	{
		return hadr[0][0];
	}
	return hadr[i-1][j-1];
}

long double CMatrix::operator ()(int i,int j) const
{
	if (i<1 || i>xDim || j<1 || j>yDim)
	{
      return 0.0L;
	}
	return hadr[i-1][j-1];
}

void CMatrix::SetDim(int d1,int d2)
{
	if (sub == 0)
	{
		for (int i=1;i <= xDim;i++) 
			delete [] hadr[i-1];
		delete [] hadr;
		xDim = d1;
		yDim = d2;
		hadr = new long double*[xDim];
		for (i=1;i <= xDim;i++) 
			hadr[i-1] = new long double[yDim];
		Clear();
	}
}

void CMatrix::ForceDimension(int d1,int d2)
{
	CMatrix A = *this;
	SetDim(d1, d2);
	for (int i=1; i<=Zeilen() && i<=A.Zeilen(); i++)
	{
		for (int j=1; j<=Spalten() && j<=A.Spalten(); j++)
		{
			SetC(i, j, A(i, j));
		}
	}
}

CMatrix CMatrix::operator() (int i,int j,int k,int l)
{
	ASSERT(i<=j);
	ASSERT(k<=l);
	ASSERT(i>=1);
	ASSERT(j<=xDim);
	ASSERT(k>=1);
	ASSERT(l<=yDim);
	CMatrix M(j-i+1,l-k+1,1);
	for (int n=i;n <= j;n++) 
		M.hadr[n-i] = &hadr[n-1][k-1];
	return M;
}

CMatrix CMatrix::GetSubMatrix(int i,int j,int k,int l) const
{
	ASSERT(i<=j);
	ASSERT(k<=l);
	ASSERT(i>=1);
	ASSERT(j<=xDim);
	ASSERT(k>=1);
	ASSERT(l<=yDim);
	CMatrix M(j-i+1,l-k+1,1);
	for (int n=i; n<=j; n++) 
	{
		for (int m=k; m<=l; m++)
		{
			M.hadr[n-i][m-1] = hadr[n-1][k+m-1];
		}
	}
	return M;
}

CMatrix& CMatrix::operator =(const CMatrix& M)
{
	Status = M.Status;
	Name = M.Name;
	if (&M != this)
	{
		if ((M.xDim != xDim) || (M.yDim != yDim)) SetDim(M.xDim,M.yDim);
		for (int i=1;i <= M.xDim;i++)
		{
			memmove(hadr[i-1],M.hadr[i-1],yDim*sizeof(long double)) ;
		}
	}
	return *this;
}

CMatrix& CMatrix::operator =(long double x)
{
	SetDim(1,1);
	hadr[0][0] = x;
	return *this;
}

CMatrix& CMatrix::operator +=(const CMatrix& M)
{
	ASSERT(xDim==M.xDim);
	ASSERT(yDim==M.yDim);
	for (int i=1;i <= xDim;i++)
	{
		for (int j=1;j <= yDim;j++)
		{
			hadr[i-1][j-1] += M.hadr[i-1][j-1];
		}
	}
	return (*this);
}

CMatrix& CMatrix::operator -=(const CMatrix& m)
{
	ASSERT(xDim==m.xDim);
	ASSERT(yDim==m.yDim);
	for (int i=1;i <= xDim;i++)
	{
		for (int j=1;j <= yDim;j++)
		{
			hadr[i-1][j-1] -= m.hadr[i-1][j-1];
		}
	}
	return (*this);
}

CMatrix& CMatrix::operator *=(long double a)
{
	for (int i=1;i <= xDim;i++)
	{
		for (int j=1;j <= yDim;j++)
		{
			hadr[i-1][j-1] *= a;
		}
	}
	return (*this);
}

CMatrix& CMatrix::operator /=(long double a)
{
	for (int i=1;i <= xDim;i++)
	{
		for (int j=1;j <= yDim;j++)
		{
			hadr[i-1][j-1] /= a;
		}
	}
	return (*this);
}

CMatrix& CMatrix::operator =(const CVektor& v)
{
	Status = Normal;
	Name = v.Name;
	if (v.Status == CVektor::ZeilenVektor)
	{
		if ((yDim  != v.Dim) || (xDim != 1)) SetDim(1,v.Dim);
		memmove(hadr[0],v.hadr,yDim*sizeof(long double));
	}
	if (v.Status == CVektor::SpaltenVektor)
	{
		if ((xDim  != v.Dim) || (yDim != 1)) SetDim(v.Dim,1);
		for (int i=1;i <= xDim;i++) 
			hadr[i-1][0] = v.hadr[i-1];
	}
	return (*this);
}

int operator==(const CMatrix& A, const CMatrix& B)
{
	if (A.Zeilen() != B.Zeilen())
		return 0;
	if (A.Spalten() != B.Spalten())
		return 0;
	for (int i=1; i<=A.Zeilen(); i++)
	{
		for (int j=1; j<=A.Spalten(); j++)
		{
			if (A(i, j) != B(i,j))
				return 0;
		}
	}
	return 1;
}

int operator!=(const CMatrix& A, const CMatrix& B)
{
	return !(A == B);
}

void CMatrix::Clear()
{
	Status = Normal;
	for (int i=1;i <= xDim;i++) 
		memset(hadr[i-1],0,yDim*sizeof(long double));
}

void CMatrix::EinheitsMatrix()
{
	Clear();
	for (int i=1; i <= min(xDim, yDim); i++)
	{
		hadr[i-1][i-1] = 1.0L;
	}
}

void CMatrix::EinsMatrix()
{
	Clear();
	for (int i=1; i <= xDim; i++)
	{
		for (int j=1; j <= yDim; j++)
		{
			hadr[i-1][j-1] = 1.0L;
		}
	}
}

void CMatrix::NegEinheitsMatrix()
{
	Clear();
	for (int i=1; i <= min(xDim, yDim); i++)
	{
		hadr[i-1][i-1] = -1.0L;
	}
}

void CMatrix::OberesDreieck()
{
	for (int i=2; i <= xDim; i++)
	{
		for (int j=1; j<i && j<=yDim; j++)
		{
			hadr[i-1][j-1] = 0.0L;
		}
	}
}

void CMatrix::UnteresDreieck()
{
	for (int i=1; i < xDim; i++)
	{
		for (int j=i; j<=yDim; j++)
		{
			hadr[i-1][j-1] = 0.0L;
		}
	}
}

CMatrix NullMatrix(const CMatrix& M)
{
	return CMatrix(M.Zeilen(), M.Spalten());
}

CMatrix NullMatrix(int dim)
{
	return CMatrix(dim, dim);
}

CMatrix NullMatrix(int zeilen, int spalten)
{
	return CMatrix(zeilen, spalten);
}


CMatrix EinheitsMatrix(const CMatrix& M)
{
	CMatrix N(M.Zeilen(), M.Spalten());
	N.EinheitsMatrix();
	return N;
}

CMatrix EinheitsMatrix(int dim)
{
	CMatrix N(dim, dim);
	N.EinheitsMatrix();
	return N;
}

CMatrix EinheitsMatrix(int zeilen, int spalten)
{
	CMatrix N(zeilen, spalten);
	N.EinheitsMatrix();
	return N;
}

CMatrix NegEinheitsMatrix(const CMatrix& M)
{
	CMatrix N(M.Zeilen(), M.Spalten());
	N.NegEinheitsMatrix();
	return N;
}

CMatrix NegEinheitsMatrix(int dim)
{
	CMatrix N(dim, dim);
	N.NegEinheitsMatrix();
	return N;
}

CMatrix NegEinheitsMatrix(int zeilen, int spalten)
{
	CMatrix N(zeilen, spalten);
	N.NegEinheitsMatrix();
	return N;
}


CMatrix OberesDreieck(const CMatrix& M)
{
	CMatrix N(M);
	N.OberesDreieck();
	return N;
}

CMatrix UnteresDreieck(const CMatrix& M)
{
	CMatrix N(M);
	N.UnteresDreieck();
	return N;
}

void CVektor::Random()
{
	srand( (unsigned)time( NULL ) );
	for (int i=1; i<=Dim; i++)
	{
		hadr[i-1] = rand();
	}
}

CVektor Diag(const CMatrix& M)
{
	CVektor v(min(M.Zeilen(), M.Spalten()));
	return v;
}

// Gaxpy-Operation
CVektor Gaxpy(const CMatrix& A, const CVektor& x, const CVektor& y)
{
	CVektor v = y;
	CVektor w = A*x;
	v += w;
	return v;
}


	// Ausgabe in einen Ger�tekontext
void CMatrix::Draw(CDC* /*dc*/, CRect& /*rect*/)
{
}

void CMatrix::ResizeRows(int zeilen)
{
	if (zeilen > Zeilen()) return;
	for (int i=zeilen; i<Zeilen(); i++)
	{
		delete[] hadr[i];
		hadr[i] = 0;
	}
	xDim = zeilen;
}

CMatrix operator !(const CMatrix& M)
{
	CMatrix W(M.yDim, M.xDim);
	for (int i=1; i<=M.xDim;i++)
	{
		for (int j=1; j<=M.yDim; j++)
		{
			W(j,i) = M(i,j);
		}
	}
	return W;
}

CMatrix operator +(const CMatrix& M, const CMatrix& N)
{
	ASSERT(M.xDim==N.xDim);
	ASSERT(M.yDim==N.yDim);
	CMatrix S(M.xDim, M.yDim);
	for (int i=1; i<=M.xDim; i++)
	{
		for (int j=1; j<=M.yDim; j++)
		{
			S(i,j) = M(i,j) + N(i,j);
		}
	}
	return S;
}

CMatrix operator -(const CMatrix& M, const CMatrix& N)
{
	ASSERT(M.xDim==N.xDim);
	ASSERT(M.yDim==N.yDim);
	CMatrix S(M.xDim, M.yDim);
	for (int i=1; i<=M.xDim; i++)
	{
		for (int j=1; j<=M.yDim; j++)
		{
			S(i,j) = M(i,j) - N(i,j);
		}
	}
	return S;
}

CMatrix operator *(long double a, const CMatrix& M)
{
	CMatrix S(M.xDim, M.yDim);
	for (int i=1; i<=M.xDim; i++)
	{
		for (int j=1; j<=M.yDim; j++)
		{
			S(i,j) = a * M(i,j);
		}
	}
	return S;
}

CMatrix operator *(const CMatrix& M, long double a)
{
	CMatrix S(M.xDim, M.yDim);
	for (int i=1; i<=M.xDim; i++)
	{
		for (int j=1; j<=M.yDim; j++)
		{
			S(i,j) = a * M(i,j);
		}
	}
	return S;
}

CMatrix operator /(const CMatrix& M, long double a)
{
	CMatrix S(M.xDim, M.yDim);
	for (int i=1; i<=M.xDim; i++)
	{
		for (int j=1; j<=M.yDim; j++)
		{
			S(i,j) = M(i,j) / a;
		}
	}
	return S;
}

CMatrix operator *(const CMatrix& M, const CMatrix& N)
{
	ASSERT(M.yDim==N.xDim);
	long double s = 0.0L;
	CMatrix P(M.xDim, N.yDim);

	for (int i=1; i<=M.xDim; i++)
	{
		 for (int j=1; j<=N.yDim; j++)
		 {
				s = 0.0L;
				for (int k=1; k<=M.yDim; k++)
				{
					s += M(i,k) * N(k,j);
				}
				P(i,j) = s;
		 }
	}
	return P;
}

CVektor operator *(const CMatrix& M, const CVektor& v)
{
	ASSERT(M.yDim==v.Dim);
	long double s = 0.0L;
	CVektor p(M.xDim);
	for (int i=1; i<=M.xDim; i++)
	{
		 s = 0.0L;
		 for (int j=1; j<=M.yDim; j++)
		 {
				s += M(i,j) * v(j);
		 }
		 p(i) = s;
	}
	return p;
}

CVektor operator *(const CVektor& v, const CMatrix& M)
{
	ASSERT(v.Dim==M.yDim);
	long double s = 0.0L;
	CVektor p(M.yDim);
	for (int i=1; i<=M.yDim; i++)
	{
		 s = 0.0L;
		 for (int j=1; j<=M.xDim; j++)
		 {
				s += M(j,i) * v(j);
		 }
		 p(i) = s;
	}
	return p;
}

/*
ostream& operator <<(ostream& os,CMatrix& m)
{
	os.setf(ios::scientific,ios::floatfield);os.precision(18);
	os << "\n";
	for(int i=1;i <= m.xDim;i++)
	{
		os << " [ ";
		for (int j=1;j <= m.yDim;j++)
		{
			os  <<  m(i,j) << " " ;
		}
		os << "]" << "\n";
	}
	os.setf(ios::fixed,ios::floatfield);os.precision(0);
	return os;
}
*/

void CMatrix::Serialize(CArchive& ar)
{
	if (ar.IsStoring())
	{
		ar << Name;
		ar << xDim;
		ar << yDim;
		ar << Status;
		for (int i=1; i<=xDim; i++)
			ar.Write(&hadr[i-1], yDim*sizeof(long double));
	}
	else
	{
		for (int i=1; i<=xDim; i++)
			delete [] hadr[i-1];
		delete [] hadr;
		ar >> Name;
		ar >> xDim;
		ar >> yDim;
		int status;
		ar >> status;
		Status = (CStatus)status;
		for (i=1; i<=xDim; i++)
		{
			if (ar.Read(&hadr[i-1], yDim*sizeof(long double)) < yDim*sizeof(long double))
			{
				memset(&hadr[i-1], 0, yDim*sizeof(long double));
			}
		}
	}
}

istream&  operator >>(istream& is,CMatrix& M)
{ 
	is.setf(ios::scientific, ios::floatfield);
	is.precision(18);

	int xd,yd;

	is >> xd;
	is >> yd;

	M.SetDim(xd,yd);

	is >> M.RSN;

	for (int i=1; i<=xd; i++)
	{
		for (int j=1; j<=yd; j++)
		{
			is >> M(i,j);
		}
	}

	return is;
}

ostream&  operator <<(ostream& os, CMatrix& M)
{ 
	os.setf(ios::scientific, ios::floatfield);
	os.precision(18);
	if (os.bad()) throw;//(xmsg("Write failed"));
	if (M.sub == 1)
	{
		::MessageBox(NULL, _T("Error : operator <<(ostream  &  os,CMatrix  &  m )"), _T("Cannot stream SubMatrix"), MB_OK);
		return os;//exit(1);
	}

	os << M.xDim; if (os.bad()) throw;//(xmsg("Write failed"));
	os << "\n"; if (os.bad()) throw;//(xmsg("Write failed"));
	os << M.yDim; if (os.bad()) throw;//(xmsg("Write failed"));
	os << "\n"; if (os.bad()) throw;//(xmsg("Write failed"));
	os << M.RSN; if (os.bad()) throw;//(xmsg("Write failed"));
	os << "\n"; if (os.bad()) throw;//(xmsg("Write failed"));
	
	for (int i=1;i <= M.xDim;i++)
	{
		for(int j=1;j <= M.yDim;j++)
		{
			 os << M(i,j); if (os.bad()) throw;//(xmsg("Write failed"));
			 os << "\n"; if (os.bad()) throw;//(xmsg("Write failed"));
		}
	}

	return os;
}

ostream& CMatrix::SaveBinary(ostream& os)
{
	os.write((const char*)(&xDim), sizeof(xDim));
	os.write((const char*)(&yDim), sizeof(yDim));
	os.write((const char*)(&Status), sizeof(Status));
	os.write((const char*)(&sub), sizeof(sub));
	int NameLen = Name.GetLength();
	os.write((const char*)(&NameLen), sizeof(NameLen));
	os.write((const char*)(&Name), NameLen);
	os.write((const char*)(&RSN), sizeof(RSN));
	for (int i=1;i <= xDim;i++)
		os.write((const char*)hadr[i-1], yDim*sizeof(long double));
	return os;
}

istream& CMatrix::LoadBinary(istream& is)
{
	is.read((char*)(&xDim), sizeof(xDim));
	is.read((char*)(&yDim), sizeof(yDim));
	is.read((char*)(&Status), sizeof(Status));
	is.read((char*)(&sub), sizeof(sub));
	int NameLen;
	is.read((char*)(&NameLen), sizeof(NameLen));
	char* buffer = new char[NameLen+1];
	memset(buffer, 0, NameLen+1);
	is.read(buffer, NameLen);
	Name = buffer;
	is.read((char*)(&RSN), sizeof(RSN));
	delete [] perm;
	perm = NULL;
	delete [] hadr;
	hadr = new long double*[xDim];
	for (int i=1;i <= xDim;i++) 
	{
		hadr[i-1] = new long double[yDim];
		is.read((char*)hadr[i-1], yDim*sizeof(long double));
	}
	return is;
}

void CMatrix::RowSwap(int a,int b)
{
	ASSERT(a>=1 && a<=xDim);
	ASSERT(b>=1 && b<=yDim);
	long double* p;
	if ( a != b )
	{
		RSN++;
		p = hadr[a-1];
		hadr[a-1] = hadr[b-1];
		hadr[b-1] = p;
	}
}

/*
//************************CVektorListNode*************************************

CVektorListNode::CVektorListNode(CVektor& Data)
{ ContainerData = Data;}

CVektorListNode::CVektorListNode()
{ }

CVektorListNode::~CVektorListNode()
{ }

CVektor& CVektorListNode::GetC()
{ return ContainerData;}

void CVektorListNode::SetC(CVektor& Data)
{ ContainerData = Data;}

//***************************CVektorList******************************************


void CVektorList::Insert(CVektorListNode *x,CVektor& Data)
{  CVektorListNode *t,*Nlong double;
   t = new CVektorListNode(Data);
   Nlong double       = x->Nlong double;
   x->Nlong double    = t;
   Nlong double->Prev = t;
   t->Prev    = x;
   t->Nlong double    = Nlong double;
	Nr_of_Nodes++;
}

void CVektorList::InsertCNNlong double(CVektor& Data)
{ Insert(CN,Data);}

void CVektorList::InsertCNPrev(CVektor& Data)
{ Insert(CN->Prev,Data);}

void CVektorList::InsertCNNlong doubleCN(CVektor& Data)
{ Insert(CN,Data);CN = CN->Nlong double;}

void CVektorList::InsertCNPrevCN(CVektor& Data)
{ Insert(CN->Prev,Data);CN = CN->Prev;}

void CVektorList::Kill(CVektorListNode *x)
{  CVektorListNode *Nlong double,*OverNlong double;
	Nlong double=x->Nlong double;
   if (Nlong double != Head) {
       OverNlong double=Nlong double->Nlong double;
       x->Nlong double=OverNlong double;
       OverNlong double->Prev=x;
       delete Nlong double;
       Nr_of_Nodes--;
   }
}

void CVektorList::KillCNNlong double()
{ Kill(CN);}

void CVektorList::KillCNPrev()
{ Kill(CN->Prev->Prev);}

void CVektorList::KillCN()
{ CN=CN->Nlong double;Kill(CN->Prev);}

CVektor CVektorList::Del(CVektorListNode *x)
{  CVektorListNode *Nlong double,*OverNlong double;
   CVektor Data;
   Nlong double=x->Nlong double;
   if (Nlong double != Head){
      OverNlong double=Nlong double->Nlong double;
      x->Nlong double=OverNlong double;
      OverNlong double->Prev=x;
      Data = Nlong double->GetC();
      delete Nlong double;
		Nr_of_Nodes--;
   }
   return Data;
}

CVektor CVektorList::DelCNNlong double()
{ return Del(CN);}

CVektor CVektorList::DelCNPrev()
{ return Del(CN->Prev->Prev);}

CVektor CVektorList::DelCN()
{ CN=CN->Nlong double;return Del(CN->Prev);}

int CVektorList::Empty()
{ return Nr_of_Nodes==0 ? 1 : 0;}

int CVektorList::CN_equal_Head()
{ return Head==CN ? 1 : 0 ;}

int CVektorList::CN_equal_PrevHead()
{ return Head->Prev==CN ? 1 : 0;}

void CVektorList::Clear()
{ while (Empty() != 1) Kill(Head);
  Nr_of_Nodes=0;
  CN=Head;
}

int CVektorList::Get_Nr_of_Nodes()
{ return Nr_of_Nodes;}

CVektor& CVektorList::First_GetC()
{ if (Empty()==1) Error("Liste ist leer !");
  return (Head->Nlong double->GetC());
}

void CVektorList::First_SetC(CVektor& Data)
{ if (Empty()==1) Error("Liste ist leer !");
  Head->Nlong double->SetC(Data);
}

CVektor& CVektorList::Last_GetC()
{ if (Empty()==1) Error("Liste ist leer !");
  return (Head->Prev->GetC());
}

void CVektorList::Last_SetC(CVektor& Data)
{ if (Empty()==1) Error("Liste ist leer !");
  Head->Prev->SetC(Data);
}

CVektor& CVektorList::CN_GetC()
{ if (CN == Head)  Error("Lesender Datenzugriff auf Listenkopf");
  return (CN->GetC());
}

void CVektorList::CN_SetC(CVektor& Data)
{ if (CN == Head) Error("Schreibender Datenzugriff auf Listenkopf !");
  CN->SetC(Data);
}

CVektor& CVektorList::CNN_GetC()
{ if (CN->Nlong double == Head) Error("Lesender Datenzugriff auf Listenkopf !");
  return (CN->Nlong double->GetC());
}

void CVektorList::CNN_SetC(CVektor& Data)
{ if (CN->Nlong double == Head) Error("Schreibender Datenzugriff auf Listenkopf !");
  CN->SetC(Data);
}

void CVektorList::ResetCN()
{ CN = Head->Nlong double;}

void CVektorList::Nlong doubleNode()
{ CN = CN->Nlong double;}

void CVektorList::PrevNode()
{ CN = CN->Prev;}

void CVektorList::FirstNode()
{ CN = Head->Nlong double;}

void CVektorList::LastNode()
{ CN = Head->Prev;}

void CVektorList::SetCN(CVektorListNode *p)
{ CN = p;}

CVektorList::CVektorList()
{ Head = new CVektorListNode();
  Head->Prev  = Head;
  Head->Nlong double  = Head;
  Nr_of_Nodes = 0;
  ResetCN();
}

CVektorList::CVektorList(int nr)
{ CVektor a;
  Head = new CVektorListNode();
  Head->Prev  = Head;
  Head->Nlong double  = Head;
  Nr_of_Nodes = 0;
  ResetCN();
  for(int i=1;i <= nr;i++) InsertCNNlong doubleCN(a);
  ResetCN();
}

CVektorList::CVektorList(CVektorList  & L)
{ Head = new CVektorListNode();
  Head->Prev  = Head;
  Head->Nlong double  = Head;
  Nr_of_Nodes = 0;
  CVektorListIterator i(L);
  while (i) { InsertCNNlong doubleCN(i.Current());i++;}
  ResetCN();
}

CVektorList::~CVektorList()
{ Clear();
  delete Head;
}

CVektorList  & CVektorList::operator =(CVektorList  & L)
{ if (Nr_of_Nodes == L.Nr_of_Nodes) {
     CVektorListIterator i(*this);
     CVektorListIterator j(L    );
	  while (i) { i.Current() = j.Current();i++,j++;}}
  else {
	  Error("Error : CVektorList::operator =(CVektorList& L) : Number of Nodes does not CMatrixch ! ");
  }
  return L;
}

ostream& operator <<(ostream& os,CVektorList& L)
{ CVektorListIterator i(L);
  while (i) { os << i.Current() << "  ";i++;}
  return os;
}

//************************TListIterator*************************************

CVektorListIterator::CVektorListIterator(CVektorList  & L)
{ List = &L;
  Restart();
}

CVektorListIterator::CVektorListIterator(CVektorListIterator  & It)
{ List      = It.List;
  ItCurrent = It.ItCurrent;
}

CVektorListIterator  & CVektorListIterator::operator =(CVektorListIterator  & It)
{ List      = It.List;
  ItCurrent = It.ItCurrent;
  return It;
}

CVektor& CVektorListIterator::Current()
{ if (!(ItCurrent==List->Head)) return ItCurrent->GetC();
  else return Dummy;
}

void CVektorListIterator::Restart()
{ ItCurrent=List->Head->Nlong double;}

CVektorListIterator::operator int()
{ return (ItCurrent==List->Head) ? 0 : 1;}

CVektorListIterator  & CVektorListIterator::operator ++( int )
{ ItCurrent=ItCurrent->Nlong double;return (*this);}

CVektorListIterator CVektorListIterator::operator ++()
{ CVektorListIterator temp(*this);ItCurrent=ItCurrent->Nlong double;return temp;}

//***************Normen f�r Vektorlisten**********************************

long double Norm(CVektorList  & L,int n)
{  long double  n1;
	CVektorListIterator i(L);
	switch (n) {
		case 1:{ n1 = 0.0L;
			 while (i) { n1 += Norm(i.Current(),1);i++;}
			 break;}
		case 2:{ n1 = 0.0L;
			 while (i) { n1+=  SqrNorm2(i.Current());i++;}
			 n1 = sqrtl(n1);
			 break;}
		case 3:{ n1 = 0.0L;
			 while (i) { n1 = max(n1,Norm(i.Current(),3));i++;}
			 break;
		  }
	}
   return n1;
}

long double SqrNorm2(CVektorList  & L)
{ long double n1 = 0.0L;
  CVektorListIterator i(L);
  while (i) { n1+=SqrNorm2(i.Current());i++;}
  return n1;
}
*/

// Gleichungssysteme
CVektor LU_Loesung(CMatrix& A, const CVektor& b, bool ignorestatus)
{
	CVektor x(b);
	if (LU(A, ignorestatus) && Vorw_Subst1(A, b, x, ignorestatus) && Rueckw_Subst(A, b, x, ignorestatus)) return x;
	return x;
}

bool CMatrix::LU(bool ignorestatus)
{
	if (!ignorestatus)
	{
		if (Status == LU_Faktorisiert) return true; // schon faktorisiert
		if (Status != Normal) return false;// nicht unterst�tzt, insbesondere Invalid
	}
//	RechtecksCMatrixrizen werden unterst�tzt if (Zeilen() != Spalten()) return false;
	int m=xDim;
	int n=yDim;
	// Permutations-"Matrix"
	if (perm) delete [] perm;
	perm = new int[m];
	int i,j,k,piv;
	for (i=1; i<=m; i++) perm[i-1] = i;
	RSN = 0;
	// Beginn der Faktorisierung
	for (k=1;k<=n-1;k++)
	{
		// Pivotsuche
		piv=k;
		for (i=k+1;i<=m;i++)
		{
			if (fabsl(hadr[i-1][k-1])>fabsl(hadr[piv-1][k-1])) piv=i;
		}
		if (piv>k)
		{
			RowSwap(piv,k);
			int tmp = perm[piv-1];
			perm[piv-1]=perm[k-1];
			perm[k-1]=tmp;
		}
		// lik(k,n)
		if (fabsl(hadr[k-1][k-1]) < 1E-19) return false;
		for (i=k+1;i<=m;i++)
		{
			hadr[i-1][k-1] /= hadr[k-1][k-1];
		}
		// mult(k,n)
		for (j=k+1;j<=n;j++)
		{
			for (i=k+1;i<=m;i++)
			{
				hadr[i-1][j-1] -= hadr[i-1][k-1]*hadr[k-1][j-1];
			}
		}
	}
	// Faktorisierung erfolgreich abgeschlossen
	Status = LU_Faktorisiert;
	return true;
}

bool Vorw_Subst1(CMatrix& A, const CVektor& b, CVektor& x, bool /*ignorestatus*/)
{
	int n=A.Spalten();
	int k;
	x.SetDim(n);
	for (k=1;k<=n;k++)
	{
		if (A.perm)
			x(k)=b(A.perm[k-1]);
		else
			x(k)=b[k];
	}
	for (k=1;k<=n;k++)
	{
		for (int j=1;j<k;j++)
		{
			x(k) -= A(k,j)*x(j);
		}
	}
	return true;
}

bool Vorw_Subst_(CMatrix& A, const CVektor& b, CVektor& x, bool /*ignorestatus*/)
{
	int n=A.Spalten();
	int k;
	x.SetDim(n);
	for (k=1;k<=n;k++)
	{
		if (A.perm)
			x(k)=b(A.perm[k-1]);
		else
			x(k)=b[k];
	}
	for (k=1;k<=n;k++)
	{
		for (int j=1;j<k;j++)
		{
			x(k) -= A(k,j)*x(j);
		}
		x(k) /= A(k, k);
	}
	return true;
}

bool Rueckw_Subst(CMatrix& A, const CVektor& /*b*/, CVektor& x, bool /*ignorestatus*/)
{
	int n=A.Spalten();
	if (x.GetDim() != n) x.SetDim(n);
	for (int k=n;k>=1;k--)
	{
		for (int j=k+1;j<=n;j++)
		{
			x(k) -= A(k,j)*x(j);
		}
		if (fabsl(A(k,k)) < 1E-19) return false;
		x(k) /= A(k,k);
	}
	return true;
}

bool CMatrix::Householder_QR(CMatrix* Q, bool ignorestatus)
{
	if (!ignorestatus)
	{
		if (Status == QR_Faktorisiert) return true; // schon faktorisiert
		if (Status == LU_Faktorisiert) return false;// nicht unterst�tzt
	}
	int m=Zeilen();
	int n=Spalten();
	int i,j,k;
	long double beta;
	CVektor v(m);
	CVektor w(m);//�nderung, fr�her w(n)
	for (j=1;j<=n;j++)
	{
		long double my=0.0L;
		for (i=j;i<=m;i++)
		{
			my += sqrl((*this)(i,j));
		}
		my = sqrtl(my);
		for (i=j;i<=m;i++)
		{
			v(i) = (*this)(i,j);
		}
		beta = (*this)(j,j) + signl((*this)(j,j))*my;
		for (i=j+1;i<=m;i++)
		{
			v(i) /= beta;
		}
		v(j)=1.0L;
		beta=0.0L;
		for (i=j;i<=m;i++)
		{
			beta += sqrl(v(i));
		}
		beta = -2.0L/beta;
		w.SetDim(n);
		w.Clear();
		for (i=j;i<=n;i++)
		{
			for (k=j;k<=m;k++)
			{
				w(i) += beta * (*this)(k,i) * v(k);
			}
		}
		for (i=j;i<=m;i++)
		{
			for (k=j;k<=n;k++)
			{
				(*this)(i,k) += v(i)*w(k);
			}
		}
		if (j<m)
		{
			for (k=j+1;k<=m;k++)
			{
				(*this)(k,j)=v(k);
			}
		}
	}
	if (Q)
	{
		// noch nicht implementiert
		Q->SetDim(m,m);
		Q->Clear();
		//Q->EinheitsMatrix();
		for (i=1; i<=m; i++)
		{
			(*Q)(i,i) = 1.0L;
		}
		for (j=n;j>=1;j--)
		{
			v(j)=1.0L;
			for (i=j+1;i<=m;i++)
			{
				v(i) = (*this)(i,j);
			}
			beta=0.0L;
			for (i=j;i<=m;i++)
			{
				beta += sqrl(v(i));
			}
			beta = -2.0L/beta;
			w.SetDim(m);
			w.Clear();
			for (i=j;i<=m;i++)
			{
				for (k=j;k<=m;k++)
				{
					w(i) += beta * (*Q)(k,i) * v(k);
				}
			}
			for (i=j;i<=m;i++)
			{
				for (k=j;k<=m;k++)
				{
					(*Q)(i,k) += v(i)*w(k);
				}
			}
		}
	}
	Status = QR_Faktorisiert;
	return true;
}

CMatrix inv(CMatrix& A)
{
	int n = A.xDim;
	CMatrix B(n, n);
	CVektor e(n);
	CVektor y(n);
	if (A.Zeilen() != A.Spalten()) return B;
	for (int i=1; i<=n; i++)
	{
		e.EinheitsVektor(i);
		LU_Loesung(A, e, y);
		for (int j=1; j<=n; j++)   //effektiver?
		{
			B(j, i) = y(j);
		}
	}
	return B;
}

void CMatrix::Random()
{
	srand( (unsigned)time( NULL ) );
	for (int i=1; i<=Zeilen(); i++)
	{
		for (int j=1; j<=Spalten(); j++)
		{
			hadr[i-1][j-1] = rand();
		}
	}
}

CMatrix SquareRoot(const CMatrix& A)
{
	CMatrix B = A;
	return B;
}

CMatrix Exp(const CMatrix& A)
{
	CMatrix B = A;
	return B;
}


void DDX_Vektor(CDataExchange* pDX, CListCtrl& control, CVektor& v, const char* Title, bool trans)
{
	int i;
	if (control.m_hWnd==NULL) return;
	CHeaderCtrl* pHeader = (CHeaderCtrl*)control.GetDlgItem(0);
	int nColumnCount = pHeader->GetItemCount();
	if (!pDX->m_bSaveAndValidate)
	{
		// Initialisieren
		control.DeleteAllItems();
		for (i=nColumnCount; i>=0; i--)
		{
			control.DeleteColumn(i);
		}
		control.InsertColumn(0, (Title==NULL) ? _T("Vektor") : Title, LVCFMT_LEFT, 30);
		if (trans)
		{
			for (i=1; i<=v.Dim; i++)
			{
				control.InsertColumn(i, CMathString(i), LVCFMT_RIGHT, 50);
			}
			control.InsertItem(1, CMathString(1));
			for (i=1; i<=v.Dim; i++)
			{
				control.SetItemText(0, i, CMathString(v[i], 15, true));
			}
		}
		else
		{
			control.InsertColumn(1, _T("A"), LVCFMT_RIGHT, 50);
			control.SetItemCount(v.Dim);
			for (i=1; i<=v.Dim; i++)
			{
				control.InsertItem(i, CMathString(i));
				control.SetItemText(i-1, 1, CMathString(v[i], 15, true));
			}
		}
	}
	else
	{
		if (trans)
		{
			v.SetDim(nColumnCount-1);
			for (i=1; i<=v.Dim; i++)
			{
				CMathString str;
				str = control.GetItemText(0, i);
				v[i] = str.Value();
			}
		}
		else
		{
			v.SetDim(control.GetItemCount());
			for (i=1; i<=v.Dim; i++)
			{
				CMathString str;
				str = control.GetItemText(i-1, 1);
				v[i] = str.Value();
			}
		}
	}
}

void DDX_Matrix(CDataExchange* pDX, CListCtrl& control, CMatrix& M, const char* Title)
{
	int i, j;
	ASSERT(control.m_hWnd != NULL);
	if (control.m_hWnd==NULL) return;

	CHeaderCtrl* pHeader = (CHeaderCtrl*)control.GetDlgItem(0);
	ASSERT(pHeader != NULL);
	int nColumnCount = pHeader->GetItemCount();
	if (!pDX->m_bSaveAndValidate)
	{
		// Initialisieren
		for (i=nColumnCount; i>=0; i--)
		{
			control.DeleteColumn(i);
		}
		control.DeleteAllItems();
		control.InsertColumn(0, (Title==NULL) ? _T("Matrix") : Title, LVCFMT_LEFT, 30);
		control.SetItemCount(M.Zeilen());
		for (i=1; i<=M.Spalten(); i++)
		{
			control.InsertColumn(i, CMathString(i), LVCFMT_RIGHT, 50);
		}
		for (j=1; j<=M.Zeilen(); j++)
		{
			control.InsertItem(j, CMathString(j));
			for (i=1; i<=M.Spalten(); i++)
			{
				control.SetItemText(j-1, i, CMathString(M(j, i), 15, true));
			}
		}
	}
	else
	{
		// Get the number of columns
		M.SetDim(control.GetItemCount(), nColumnCount-1);
		for (j=1; j<=M.Zeilen(); j++)
		{
			for (i=1; i<=M.Spalten(); i++)
			{
				CMathString str;
				str = control.GetItemText(j-1, i);
				M(j, i) = str.Value();
			}
		}
	}
}

void CMatrix::Define_Matrix(int d1, ...)
{
	int z, s;
	long double x;
	va_list marker;

	SetDim(d1, d1);

	va_start(marker, d1);     /* Initialize variable arguments. */
	for (z=0; z<d1; z++)
	{
		for (s=0; s<d1; s++)
		{
			x = va_arg(marker, long double);
			hadr[z][s] = x;
		}
	}
	va_end(marker);              /* Reset variable arguments.      */
}

void CMatrix::Define_UMatrix(int d1, ...)
{
	int z, s;
	long double x;
	va_list marker;

	SetDim(d1, d1);

	va_start(marker, d1);     /* Initialize variable arguments. */
	for (z=0; z<d1; z++)
	{
		for (s=0; s<z; s++)
		{
			x = va_arg(marker, long double);
			hadr[z][s] = x;
		}
	}
	va_end(marker);              /* Reset variable arguments.      */
}

long double CMatrix::FrobeniusNorm() const
{
	long double nrm = 0.0L;
	for (int i=1; i<=Zeilen(); i++)
	{
		for (int j=1; j<=Spalten(); j++)
		{
			nrm += sqrl(fabsl(GetC(i, j)));
		}
	}
	return sqrtl(nrm);
}

long double CMatrix::ZeilensummenNorm() const
{
	long double nrm = 0.0L;
	for (int i=1; i<=Zeilen(); i++)
	{
		long double summe = 0.0L;
		for (int j=1; j<=Spalten(); j++)
		{
			summe += fabsl(GetC(i, j));
		}
		if (summe>=nrm)
			nrm = summe;
	}
	return nrm;
}

long double CMatrix::SpaltensummenNorm() const
{
	long double nrm = 0.0L;
	for (int j=1; j<=Spalten(); j++)
	{
		long double summe = 0.0L;
		for (int i=1; i<=Zeilen(); i++)
		{
			summe += fabsl(GetC(i, j));
		}
		if (summe>=nrm)
			nrm = summe;
	}
	return nrm;
}

long double CMatrix::GesamtNorm() const
{
	long double nrm = 0.0L;
	for (int j=1; j<=Spalten(); j++)
	{
		for (int i=1; i<=Zeilen(); i++)
		{
			nrm += fabsl(GetC(i, j));
		}
	}
	return nrm;
}

long double CMatrix::SpektralNorm() const
{
	long double nrm = 0.0L;
	return nrm;
}


long double FrobeniusNorm(const CMatrix& A)
{
	long double nrm = 0.0L;
	for (int i=1; i<=A.Zeilen(); i++)
	{
		for (int j=1; j<=A.Spalten(); j++)
		{
			nrm += sqrl(fabsl(A(i, j)));
		}
	}
	return sqrtl(nrm);
}

long double ZeilensummenNorm(const CMatrix& A)
{
	long double nrm = 0.0L;
	for (int i=1; i<=A.Zeilen(); i++)
	{
		long double summe = 0.0L;
		for (int j=1; j<=A.Spalten(); j++)
		{
			summe += fabsl(A(i, j));
		}
		if (summe>=nrm)
			nrm = summe;
	}
	return nrm;
}

long double SpaltensummenNorm(const CMatrix& A)
{
	long double nrm = 0.0L;
	for (int j=1; j<=A.Spalten(); j++)
	{
		long double summe = 0.0L;
		for (int i=1; i<=A.Zeilen(); i++)
		{
			summe += fabsl(A(i, j));
		}
		if (summe>=nrm)
			nrm = summe;
	}
	return nrm;
}

long double GesamtNorm(const CMatrix& A)
{
	long double nrm = 0.0L;
	for (int j=1; j<=A.Spalten(); j++)
	{
		for (int i=1; i<=A.Zeilen(); i++)
		{
			nrm += fabsl(A(i, j));
		}
	}
	return nrm;
}

long double SpektralNorm(const CMatrix& A)
{
	long double nrm = 0.0L;
	return nrm;
}

CVektor CMatrix::GetRow(int row) const
{
	ASSERT(row>=1 && row<=Zeilen());
	CVektor v(Spalten(), _T("Zeile"), true);
	for (int i=1; i<=Spalten(); i++)
	{
		v[i] = GetC(row, i);
	}
	return v;
}

CVektor CMatrix::GetCol(int col) const
{
	ASSERT(col>=1 && col<=Spalten());
	CVektor v(Zeilen());
	for (int i=1; i<=Zeilen(); i++)
	{
		v[i] = GetC(i, col);
	}
	return v;
}

CVektor CMatrix::ZeilenSummen() const
{
	CVektor v(Zeilen());
	for (int i=1; i<=Zeilen(); i++)
	{
		v[i] = 0.0L;
		for (int j=1; j<=Spalten(); j++)
		{
			v[i] += GetC(i, j);
		}
	}
	return v;
}

CVektor CMatrix::SpaltenSummen() const
{
	CVektor v(Spalten());
	for (int i=1; i<=Spalten(); i++)
	{
		v[i] = 0.0L;
		for (int j=1; j<=Zeilen(); j++)
		{
			v[i] += GetC(j, i);
		}
	}
	return v;
}

CVektor CMatrix::Diag() const
{
	int dim = min(Zeilen(), Spalten());
	CVektor v(dim);
	for (int i=1; i<=dim; i++)
	{
		v[i] = GetC(i, i);
	}
	return v;
}

CMatrix sin(const CMatrix& M)
{
	CMatrix N(M);
	N.sin();
	return N;
}

CMatrix cos(const CMatrix& M)
{
	CMatrix N(M);
	N.cos();
	return N;
}

CMatrix tan(const CMatrix& M)
{
	CMatrix N(M);
	N.tan();
	return N;
}

CMatrix cot(const CMatrix& M)
{
	CMatrix N(M);
	N.cot();
	return N;
}

CMatrix abs(const CMatrix& M)
{
	CMatrix N(M);
	N.abs();
	return N;
}

CMatrix UserFunction(const CMatrix& M, const char* FunktionsTerm, const char* Variable)
{
	CMatrix N(M);
	N.UserFunction(FunktionsTerm, Variable);
	return N;
}

void CMatrix::sin()
{
	for (int i=0; i<xDim; i++)
	{
		for (int j=0; j<yDim; j++)
		{
			hadr[i][j] = ::sinl(hadr[i][j]);
		}
	}
}

void CMatrix::cos()
{
	for (int i=0; i<xDim; i++)
	{
		for (int j=0; j<yDim; j++)
		{
			hadr[i][j] = ::cosl(hadr[i][j]);
		}
	}
}

void CMatrix::tan()
{
	for (int i=0; i<xDim; i++)
	{
		for (int j=0; j<yDim; j++)
		{
			hadr[i][j] = ::tanl(hadr[i][j]);
		}
	}
}

void CMatrix::cot()
{
	for (int i=0; i<xDim; i++)
	{
		for (int j=0; j<yDim; j++)
		{
			hadr[i][j] = ::cosl(hadr[i][j]) / ::sinl(hadr[i][j]);
		}
	}
}

void CMatrix::abs()
{
	for (int i=0; i<xDim; i++)
	{
		for (int j=0; j<yDim; j++)
		{
			hadr[i][j] = ::fabsl(hadr[i][j]);
		}
	}
}

void CMatrix::UserFunction(const char* FunktionsTerm, const char* Variable)
{
	CString FktTerm(FunktionsTerm);
	CFunction<long double>* fkt = CFunction<long double>::Parse(FktTerm, Variable);

	if (fkt==NULL)
		return;

	for (int i=0; i<xDim; i++)
	{
		for (int j=0; j<yDim; j++)
		{
			hadr[i][j] = fkt->Execute(hadr[i][j]);
		}
	}
	delete fkt;
}

bool LU(CMatrix& A, CMatrix& L, CMatrix& U, bool ignorestatus)
{
	int i,j;
	int m = A.xDim;
	int n = A.yDim;
	if (!LU(A, ignorestatus)) return false;
	// Aufbau von U
	U.Clear();
	for (i=1;i <= m;i++)
	{
		for (int j=i;j <= n;j++)
		{
			U(i,j) = A(i,j);
		}
	}
	// Aufbau von L
	L.Clear();
	for (i=1;i <= m;i++)
	{
		for (j=1;j < i;j++)
		{
			L(i,j) = A(i,j);
		}
	}
	for (i=1;i <= m;i++) L(i,i) = 1.0L;
	return true;
}

long double Det(CMatrix& A)
{  
	long double d = 1.0L;
	if (!LU(A)) return 0.0L;
	for (int i=1; i<=A.xDim; i++)
	{
		d *= A(i,i);
	}
	if (A.RSN % 2 != 0) d = -d;
	return d;
}

long double CMatrix::Det()
{  
	long double d = 1.0L;
	if (!LU()) return 0.0L;
	for (int i=1; i<=xDim; i++)
	{
		d *= GetC(i,i);
	}
	if (RSN % 2 != 0) d = -d;
	return d;
}

// QR-Faktorisierung der Matrix A
bool _QR(CMatrix& A, bool /*ignorestatus*/)
{
	int i;
	int m = A.xDim;
	int n = A.yDim;
	CVektor v(m);
	for (i=1;i <= n;i++)
	{
		 if (i < m)
		 {
				house(CVektor(A(i,m,i,i)),v(i,m));
				rowhouse(A(i,m,i,n),v(i,m));
				A(i+1,m,i,i) = v(i+1,m);
		 }
	}
	return true;
}

bool _QR(CMatrix& A, CMatrix* Q, CMatrix* R, bool /*ignorestatus*/)
{
	int i;
	int m = A.xDim;
	int n = A.yDim;
	CVektor v(m);
	// Berechnung der QR-Faktorisierung von A
	_QR(A);
	if (!Q) return true;
	// R�ckw�rtsakkumulation von Q
	Q->SetDim(m,m);
	Q->Clear();
	for (i=1;i <= m;i++) (*Q)(i,i) = 1.0L;
	for (i=n;i >= 1;i--)
	{
		 if (i < m)
		 {
			v(i) = 1.0L;
			v(i+1,m) = A(i+1,m,i,i);
			rowhouse((*Q)(i,m,i,m),v(i,m));
		 }
	}
	if (!R) return true;
	// Aufbau von R
	R->SetDim(m, n);
	R->Clear();
	for (i=1;i <= m;i++)
	{
		 for (int j=i;j <= n;j++)
		 {
				(*R)(i,j) = A(i,j);
		 }
	}
	return true;
}

bool QR_Solve(CMatrix& A, CVektor& b, CVektor& c, long double TOL)
{
	int m = A.xDim;
	int n = A.yDim;

	CVektor v(m);
	CVektor y(m);

	// gegebenenfalls Berechnung der QR-Faktorisierung
	if (!A.QR()) return false;
	// Initialisierung des CVektors y mit der rechten Seite b
	y = b;
	// Berechnung von !Q * y
	for (int i=1;i <= n;i++)
	{
		if (i < m)
		{
			v(i+1,m) = A(i+1,m,i,i);
			v(i) = 1.0L;
			rowhouse(y(i,m),v(i,m));
		}
	}
	// Initialisierung des CVektors c mit y
	c = y;
	// L�sung von R * c = y;
	for (int k = m;k >= 1;k--)
	{
		for (int j = k+1;j <= m;j++)
		{
			c(k) -= A(k,j) * c(j);
		}
		if (fabsl(A(k,k)) < TOL) { A.SetStatus(CMatrix::Invalid);return false;}
		else c(k) /= A(k,k);
	}
	return true;
}

// Transformation von A auf Hessenbergform mittels Householdertransformationen
bool HQR(CMatrix& A)
{
	if (A.Status == CMatrix::Hessenbergform) return true;
	if (A.Status != CMatrix::Normal) return false;
	int i;
	int m = A.xDim;
	int n = A.yDim;
	CVektor v(m);
	for (i=1;i <= n;i++)
	{
		if (i < m-1)
		{
			house(CVektor(A(i+1,m,i,i)),v(i+1,m));
			rowhouse(A(i+1,m,i,n),v(i+1,m));
			colhouse(A(1,m,i+1,n),v(i+1,n));
			A(i+2,m,i,i) = v(i+2,m);
		}
	}
	A.Status = CMatrix::Hessenbergform;
	return true;
}

// Transformation von A auf Hessenbergform mittels Householdertransformationen
bool HQR(CMatrix& A, CMatrix& Q, CMatrix& R)
{
	int i;
	int m = A.xDim;
	int n = A.yDim;
	CVektor v(m);
	// Berechnung der Hessenbergform von A mittels Householdertransformationen
	if (!HQR(A)) return false;
	// R�ckw�rtsakkumulation von Q
	Q.Clear();
	for (i=1;i <= m;i++) Q(i,i) = 1.0L;
	for (i=n;i >= 1;i--)
	{
		 if (i < m-1)
		 {
			v(i+1) = 1.0L;
			v(i+2,m) = A(i+2,m,i,i);
			rowhouse(Q(i+1,m,i+1,m),v(i+1,m));
		 }
	}
	// Aufbau von R
	R.Clear();
	for (int j=1;j <= n;j++) R(1,j) = A(1,j);
	for (i=2;i <= m;i++)
	{
		 for (int j=i-1;j <= n;j++)
		 {
				R(i,j) = A(i,j);
		 }
	}
	return true;
}

// GG-Faktorisierung der symmmetrischen positiv definiten Matrix A
bool GG(CMatrix& A, bool ignorestatus)
{
	if (!ignorestatus)
	{
		if (A.Status == CMatrix::GG_Faktorisiert) return true;
		if (A.Status != CMatrix::Normal) return false;
	}
	int i,j,k;
	int n = A.yDim;
	for (k=1;k <= n;k++)
	{
		 A(k,k) = sqrtl(A(k,k));
		 // A(k+1,n,k,k) = A(k+1,n,k,k) / A(k,k);
		 for (j=k+1;j <= n;j++)
		 {
				A(j,k) /= A(k,k);
		 }
		 for (j=k+1;j <= n;j++)
		 {
			// A(j,n,j,j) = A(j,n,j,j) - A(j,n,k,k)*A(j,k);
			for (i=j;i <= n;i++)
			{
				A(i,j) -= A(i,k) * A(j,k);
			}
		 }
	}
	A.Status = CMatrix::GG_Faktorisiert;
	return true;
}

bool GG(CMatrix& A, CMatrix& G, bool ignorestatus)
{
	int i,j;
	int n = A.yDim;
	// GG-Faktorisierung von A
	if (!GG(A, ignorestatus)) return false;
	// Aufbau von G
	G.Clear();
	for (i=1;i <= n;i++)
	{
		for (j=1;j <= i;j++)
		{
			G(i,j) = A(i,j);
		}
	}
	return true;
}

bool GG_Solve(CMatrix& A, CVektor& b, CVektor& c, long double TOL)
{
	int j,k;
	int n = A.xDim;
	CMatrix G(n,n);
	CVektor y(n);

	// gegebenenfalls Berechnung der GG-Faktorisierung
	if (!GG(A,G)) {A.SetStatus(CMatrix::Invalid); return false;}
	// Initialisierung von y mit b
	y = b;
	// Berechnung der L�sung von G * y = b
	for (k=1;k <= n;k++){
		 for (j=1;j <= (k-1);j++)
		 {
				y(k) -= A(k,j) * y(j);
		 }
		 if (fabsl(A(k,k)) < TOL) {A.SetStatus(CMatrix::Invalid);return false;}
		 else y(k) /= A(k,k);
	}
	// Initialisierung von c mit y
	c = y;
	// L�sung von !G * c = y
	for (k = n;k >= 1;k--)
	{
		 for (j = k+1;j <= n;j++)
		 {
				c(k) -= A(j,k) * c(j);
		 }
		 if (fabsl(A(k,k)) < TOL) {A.SetStatus(CMatrix::Invalid);return false;}
		 else c(k) /= A(k,k);
	}
	return true;
}

//Bi-Diagonalisierung einer Matrix A = U * B * V
bool BIQR(CMatrix& A)
{
	if (A.Status == CMatrix::QR_Bidiagonalisiert) return true;
	if (A.Status != CMatrix::Normal) return false;
	int m = A.xDim;
	int n = A.yDim;
	CVektor v(m);
	CVektor w(n,CString(),true);

	for (int j=1;j <= n;j++)
	{
		if (j < m)
		{
			house(CVektor(A(j,m,j,j)),v(j,m));
			rowhouse(A(j,m,j,n),v(j,m));
			A(j+1,m,j,j) = v(j+1,m);
		}
		if (j <= n-2)
		{
			house(CVektor(A(j,j,j+1,n)),w(j+1,n));
			colhouse(A(j,m,j+1,n),w(j+1,n));
			A(j,j,j+2,n) = w(j+2,n);
		}
	}
	A.Status = CMatrix::QR_Bidiagonalisiert;
	return true;
}

bool BIQR(CMatrix& A, CMatrix& U, CMatrix& B, CMatrix& V)
{
	int i;
	int m = A.xDim;
	int n = A.yDim;
	CVektor v(m);
	CVektor w(n,CString(),1);
	// Berechnung der Bi-Diagonalisierung von A
	if (!BIQR(A)) return false;
	// R�ckw�rtsakkumulation von U
	U.Clear();
	for (i=1;i <= m;i++) U(i,i) = 1.0L;
	for (i=n;i >= 1;i--)
	{
		 if (i < m)
		 {
			v(i) = 1.0L;
			v(i+1,m) = A(i+1,m,i,i);
			rowhouse(U(i,m,i,m),v(i,m));
		 }
	}
	// R�ckw�rtsakkumulation von V
	V.Clear();
	for (i=1;i <= n;i++) V(i,i) = 1.0L;
	for (i=n;i >= 1;i--)
	{
		 if (i <= n-2)
		 {
			w(i+1)   = 1.0L;
			w(i+2,n) = A(i,i,i+2,n);
			colhouse(V(i+1,n,i+1,n),w(i+1,n));
		 }
	}
	// Aufbau von B
	B.Clear();
	for (i=1;i <= n-1;i++)
	{
		 B(i,i)   = A(i,i);
		 B(i,i+1) = A(i,i+1);
	}
	B(n,n) = A(n,n);
	return true;
}

// Berechnung der Singul�rwerte der quadratischen Matrix A
bool SVD(CMatrix& A, CVektor& v, long double TOL)
{
	if (A.Status != CMatrix::Normal) return false;
	int n    = A.xDim;
	int Test = 1;
	int It   = 0;
	long double k    = 0.0L;

	long double a_old,a_new;
	long double a1,a2,a3,a4;
	long double c,s,r;

	CMatrix U(n,n);
	CMatrix V(n,n);
	CMatrix B(n,n);

	// Transformation der Matrix A auf Bidiagonalform
	if (!BIQR(A,U,B,V)) return false;;

	a_old = B(n,n);
	a_new = B(n,n);

	do 
	{
		if ((n == 1) || (n==2))
		{
			if (n==1) { v(1) = fabsl(B(1,1));Test = 0;}
			if (n==2)
			{
				a1   = sqrl(B(1,1)) + sqrl(B(1,2)) + sqrl(B(2,2)) ;
				a2   = sqrl( B(1,1)*B(2,2) );
				v(1) = sqrtl( 0.5L*a1 + sqrtl(sqrl(0.5L*a1)-a2) );
				v(2) = sqrtl( 0.5L*a1 - sqrtl(sqrl(0.5L*a1)-a2) );
				Test = 0;
			}
		}
		else
		{
			c = sqrl(B(1,1)) - k;
			s = B(1,2)*B(1,1);
			r = sqrtl(sqrl(c)+sqrl(s));
			c /= r;
			s /= r;

			giv(B(1,n,1,2),c,s);
			rebidiag(B(1,n,1,n));

			a1 = fabsl(B(n-1,n  ));a2 = fabsl(B(n-2,n-1));
			a3 = fabsl(B(n  ,n  ));a4 = fabsl(B(n-1,n-1));

			// falls  B(n,n-1) oder B(n-1,n-2) gen�gend klein sind ...
			if ( min(a1,a2) <= TOL*(a3 + a4) )
			{
				if (a1 < a2)
				{
					v(n) = fabsl(B(n,n));
					n      = n - 1;
				}
				else
				{
					a1 = sqrl(B(n-1,n-1)) + sqrl(B(n-1,n)) + sqrl(B(n,n)) ;
					a2 = sqrl( B(n-1,n-1)*B(n,n) );
					v(n-1) = sqrtl( 0.5L*a1 + sqrtl(sqrl(0.5L*a1)-a2) );
					v(n  ) = sqrtl( 0.5L*a1 - sqrtl(sqrl(0.5L*a1)-a2) );
					n = n - 2;
				}
				// Shiftparameter zur�cksetzen
				k     = 0.0L;
				a_old = B(n,n);
				a_new = B(n,n);
			}
			else
			{
				a_new = B(n,n);
				// Shifparameter setzen
				if (fabsl(1.0L - a_old/a_new) <= 0.3L) k = sqrl(a_new);
				a_old = a_new;
			}
			It++;
		}
	} while ((Test != 0) && (It < 200));
	if (It > 200) return false;/*Error*/;
	return true;
}

bool CMatrix::Cholesky_SVD()
{
	return true;
}

// Berechnung des Householder Vektors v zum Vektor x
void house(CVektor& x, CVektor& v)
{
	long double m = Norm(x,2);
	long double b;
	if (m != 0.0L)
	{
		b = x(1) + signl(x(1)) * m;
		// v(2,x.Dim) = x(2,x.Dim)/b;
		for (int i=2;i<=x.GetDim();i++) 
			v(i) = x(i) / b;
	}
	v(1) = 1.0L;
}

// Berechnung des Produktes der Householder Matrix P mit der Matrix A;
// P ist im Householder Vektor v codiert
void rowhouse(CMatrix& A, const CVektor& v)
{
	CVektor w(A.yDim);
	long double s = 0.0L;
	//  b = -2.0L / ((!v) * v)
	long double b = -2.0L / SqrNorm2(v);
	//  w = b * ((!A) * v)
	for (int i=1;i <= A.yDim;i++)
	{
		s = 0.0L;
		for (int j=1;j <= A.xDim;j++) s += A(j,i)*v(j);
		w(i) = b * s;
	}
	// A = A + v * (!w);
	for (i=1;i <= A.xDim;i++)
	{
		s = v(i);
		for (int j=1;j <= A.yDim;j++)
		{
			A(i,j) += s * w(j);
		}
	}
}

// Berechnung des Produktes der Householder Matrix P mit dem CVektor a
// P ist im Householder Vektor v codiert
void rowhouse(CVektor& a, const CVektor& v)
{
	long double w = 0.0L;
	//  b = -2.0L * (!v * v)
	long double b = -2.0L / SqrNorm2(v);
	//  w = b * (!a * v)
	for (int j=1;j <= a.GetDim();j++) 
		w += a(j) * v(j);
	w *= b;
	// a = a + v * w;
	for (int i=1;i <= a.GetDim();i++) 
		a(i) += w * v(i);
}

// Berechnung des Produktes der Matrix A mit der Householder Matrix P
// P ist im Householder Vektor v codiert
void colhouse(CMatrix& A, const CVektor& v)
{
	CVektor w(A.xDim);
	long double s = 0.0L;
	//  b = -2.0L / ((!v) * v)
	long double b = -2.0L / SqrNorm2(v);
	//w = b * (A * v);
	for (int i=1;i <= A.xDim;i++)
	{
		s = 0.0L;
		for (int j=1;j <= A.yDim;j++) s += A(i,j)*v(j);
		w(i) = b * s;
	}
	// A = A + w * (!v);
	for (i=1;i <= A.xDim;i++)
	{
		s = w(i);
		for (int j=1;j <= A.yDim;j++) 
			A(i,j) += s * v(j);
	}
}

// Berechnung der Givenskoeffinzienten c und s zu a und b
void givens(long double a,long double b,long double& c,long double& s)
{
	long double t;
	if (b==0.0L)
	{
		c = 1.0L;
		s = 0.0L;
	}
	else
	{
		if (fabsl(b) > fabsl(a))
		{
			t = a/b;
			s = 1.0L/sqrtl(1.0L+sqrl(t));
			c = s*t;
		}
		else
		{
			t = b/a;
			c = 1.0L/sqrtl(1.0L+sqrl(t));
			s = c*t;
		}
	}
}

// Berechnung des Produktes der Givensrotation P mit der Matrix A
void rowgivens(CMatrix& A, long double c, long double s)
{
	long double t1,t2;
	for (int j=1;j <= A.yDim;j++)
	{
		t1 = A(1,j);
		t2 = A(2,j);
		A(1,j) = c*t1 + s*t2;
		A(2,j) = c*t2 - s*t1;
	}
}

// Berechnung des Produktes der Matrix A mit der Givensrotation P
void colgivens(CMatrix& A, long double c, long double s)
{
	long double t1,t2;
	for (int j=1;j <= A.xDim;j++)
	{
		t1 = A(j,1);
		t2 = A(j,2);
		A(j,1) = c*t1 + s*t2;
		A(j,2) = c*t2 - s*t1;
	}
}

// Berechnung des Produktes der Matrix A mit der Givensrotation P
// verwendet bei der Berechnung der Singul�rwerte von A
void giv(CMatrix& A, long double c, long double s)
{
	long double t1,t2;
	for (int i=1;i <= A.xDim;i++)
	{
		t1 = A(i,1);
		t2 = A(i,2);
		A(i,1) = c*t1 + s*t2;
		A(i,2) = s*t1 - c*t2;
	}
}

// QR-Faktorisierung von A = Q*R ( A liegt in Hessenbergform vor )
// mittels Givensrotationen
// und Berechnung von R*Q ( in A zur�ckgegeben )
void eiggiv(CMatrix& A)
{
	long double c,s;
	int n = A.yDim;
	for (int i=1;i <= n-1;i++)
	{
		 givens(A(i,i),A(i+1,i),c,s);
		 rowgivens(A(i,i+1,1,n),c,s);
		 colgivens(A(1,n,i,i+1),c,s);
	}
}

// Rebidiagonalisierung der quadratischen Matrix A
// verwendet bei der Berechnung der Sigul�rwerte von A
void rebidiag(CMatrix& A)
{
	int m = A.xDim;
	CVektor v(2);
	CVektor w(2,CString(),1);
	for (int i=1;i <= m;i++)
	{
		 if (i < m)
		 {
			house(CVektor(A(i,i+1,i,i)),v);
			rowhouse(A(i,i+1,i,m),v);
			//   A(i+1,i) = v(2);
		 }
		 if (i <= m-2)
		 {
			house(CVektor(A(i,i,i+1,i+2)),w);
			colhouse(A(i,m,i+1,i+2),w);
			//   A(i,i+2) = w(2);
		 }
	}
}