/* fMatObj.H

  object-oriented interface for OptiVec matrix functions of type "float"

  Copyright (c) 1996-2008 by OptiCode
  All Rights Reserved.
*/

#ifndef FMATOBJ_H

#include <VecObj.h>
#include <fVecObj.h>
#include <uVecObj.h>
#include <iVecObj.h>

#define cmta const matrix<float>&
#define cvta const vector<float>&
#define cta  const float&
#define cua  const unsigned&
#define cuint const unsigned

#define mta  matrix<float>&
#define vta  vector<float>&
#define ta   float&
#define ua   unsigned&

#define mt   matrix<float>
#define vt   vector<float>
#define t    float
#define tcplx  fComplex
#define tMatrix fMatrix
#define tVector fVector

#define tMatObj fMatObj
#define tVecObj fVecObj
#ifdef __BORLANDC__
#pragma warn -inl
#endif

#if defined _MSC_VER || (defined __BORLANDC__ && __BORLANDC__ >= 0x500)
namespace OptiVec{
#endif

	VECOBJ_NEW_TEMPLATE class matrix<t>
	{
#ifdef _MSC_VER
		public:
		tMatrix  mx_matr;
		tVector  vx_m0;
		unsigned u_ht, u_len;
		ui       ui_size;
#else
		tMatrix  mx_matr;
		tVector  vx_m0;
		unsigned u_ht, u_len;
		ui       ui_size;
		friend class matrix<char>;
		friend class matrix<short>;
		friend class matrix<int>;
		friend class matrix<long>;
		friend class matrix<quad>;
		friend class matrix<unsigned char>;
		friend class matrix<unsigned short>;
		friend class matrix<unsigned int>;
		friend class matrix<unsigned long>;
		friend class matrix<ui>;
		friend class matrix<float>;
		friend class matrix<double>;
		friend class matrix<extended>;
		friend class matrix<fComplex>;
		friend class matrix<dComplex>;
		friend class matrix<eComplex>;
		friend class matrix<fPolar>;
		friend class matrix<dPolar>;
		friend class matrix<ePolar>;
		friend class vector<float>;
		friend class vector<double>;
		friend class vector<extended>;
		friend class vector<fComplex>;
		friend class vector<dComplex>;
		friend class vector<eComplex>;
#endif
		public:
////////////////////////////////////////////////////////////////////////////////////////////////////
//constructors and destructors
//
		matrix(){ mx_matr=0; vx_m0=0; u_ht=0; u_len=0; ui_size=0;}

		matrix(cua ht, cua len)
		{	if((ht > 0) && (len > 0))
			{	mx_matr=MF_matrix(ht,len); vx_m0=mx_matr[0];
				u_ht=ht; u_len=len; ui_size=((ui)ht)*len;
			}
			else{	mx_matr=0; vx_m0=0; u_ht=0; u_len=0; ui_size=0;}
		}

		matrix(cua ht, cua len, cta fill)
		{	if(!ht || !len) //size == 0
			{	mx_matr=0; vx_m0=0; u_ht=0; u_len=0; ui_size=0;}
			else
			{
				u_ht=ht; u_len=len; ui_size=((ui)ht)*len;
				if(!fill) //fill == 0.0
				{	mx_matr=MF_matrix0(ht,len); vx_m0=mx_matr[0];}
				else
				{	mx_matr=MF_matrix(ht,len);  vx_m0=mx_matr[0];
					VF_equC(vx_m0,ui_size,fill);}
			}
		}

		matrix(cmta init)
		{	if(init.ui_size)
			{	mx_matr= MF_matrix(init.u_ht, init.u_len); vx_m0=mx_matr[0];
				u_ht=init.u_ht; u_len=init.u_len; ui_size=init.ui_size;
				VF_equV(vx_m0, init.vx_m0, ui_size);
			} else
			{ mx_matr=0; vx_m0=0; u_ht=0; u_len=0; ui_size=0;}
		}

		~matrix(){ if(mx_matr) M_free(mx_matr);}


////////////////////////////////////////////////////////////////////////////////////////////////////
//management
//
		void enforceSize()const{ if(!ui_size) throw OptiVec_EmptySizeErr();}
		void enforceSize(unsigned ht, unsigned len)const{
			if(u_len != len || u_ht != ht || !ui_size) throw OptiVec_SizeErr();}

		void enforceLen(unsigned len)const{ if(u_len != len) throw OptiVec_SizeErr();}
		void enforceHt(unsigned ht)const{ if(u_ht != ht) throw OptiVec_SizeErr();}

		void enforcePow2()const{ if(!isipow2(u_len) || !isipow2(u_ht)) throw OptiVec_Pow2SizeErr();}

		void setSize(cua ht, cua len)
		{	ui size;
			size = ((ui)ht)*len;
			if(!size) dealloc(); //can set size to zero
			if(u_len != len || u_ht != ht) //only allocate if different sizes
			{	if(mx_matr) M_free(mx_matr);
				mx_matr= MF_matrix(ht, len); vx_m0=mx_matr[0];
				u_len=len; u_ht=ht; ui_size= size;
			}
		}

		void setSize(cmta X)
		{	if(!X.ui_size) dealloc(); //can set size to zero
			if(u_len != X.u_len || u_ht != X.u_ht) //only allocate if different sizes
			{	if(mx_matr) M_free(mx_matr);
				mx_matr= MF_matrix(X.u_ht, X.u_len); vx_m0=mx_matr[0];
				u_len=X.u_len; u_ht=X.u_ht; ui_size= X.ui_size;
			}
		}

		ui getHt()   const { return u_ht; }
		ui getLen()  const { return u_len; }
		ui getSize() const { return ui_size; }
		tMatrix getMatrix() const { return mx_matr; }
		tVector getM0()     const { return vx_m0; }

#if !defined (_CRTDBG_MAP_ALLOC) // crtdbg.h defines free as a macro
		void free()
		{	if(mx_matr) M_free(mx_matr);
			mx_matr=0; vx_m0=0; u_ht=0; u_len=0; ui_size=0;
		}
#endif

		void dealloc()   // identical to free
		{	if(mx_matr) M_free(mx_matr);
			mx_matr=0; vx_m0=0; u_ht=0; u_len=0; ui_size=0;
		}

		ta at(unsigned m, unsigned n) const
		{
			if(u_ht<=m || u_len<=n) throw OptiVec_SmallSizeErr();
			return mx_matr[m][n];
		}


///////////////////////////////////////////////////////////////////////////////////////////////////
// Overloaded operators
//
		mta operator=(cmta matr)
		{	setSize(matr.u_ht, matr.u_len);
			enforceSize();
			VF_equV(vx_m0, matr.vx_m0, ui_size);
			return *this;}

		mta operator=(cta scalar)
		{	enforceSize();
			VF_equC(vx_m0, ui_size, scalar);
			return *this;}

		//this is the only really dangerous operator
		tVector operator[](const int& index){return mx_matr[index];}

		mta operator+=(cmta X)
		{
			enforceSize(X.u_ht, X.u_len);
			VF_accV(vx_m0, X.vx_m0, ui_size);
			return *this;}

		mta operator+=(cta scalar)
		{
			enforceSize();
			VF_addC(vx_m0, vx_m0, ui_size, scalar);
			return *this;}

		mta operator*=(cta scalar)
		{
			enforceSize();
			VF_mulC(vx_m0,vx_m0,ui_size,scalar);
			return *this;}

		mt operator+(cmta rhs)
		{	enforceSize(rhs.u_ht, rhs.u_len);
			tMatObj tmp(this->u_ht, this->u_len);
			VF_addV(tmp.vx_m0,this->vx_m0,rhs.vx_m0,tmp.ui_size);
			return tmp;}

		mt operator+(cta C)
		{	tMatObj tmp(this->u_ht, this->u_len);
			tmp.enforceSize();
			VF_addC(tmp.vx_m0,this->vx_m0,tmp.ui_size,C);
			return tmp;}

		mt operator-(cmta rhs)
		{	enforceSize(rhs.u_ht, rhs.u_len);
			tMatObj tmp(this->u_ht, this->u_len);
			VF_subV(tmp.vx_m0,this->vx_m0,rhs.vx_m0,tmp.ui_size);
			return tmp;}

		mt operator-(cta C)
		{	tMatObj tmp(this->u_ht, this->u_len);
			tmp.enforceSize();
			VF_subC(tmp.vx_m0,this->vx_m0,tmp.ui_size,C);
			return tmp;}

		mt operator*(cmta rhs)
		{	enforceLen(rhs.u_ht);
			tMatObj tmp(this->u_ht, rhs.u_len);
			tmp.enforceSize();
			MFmulM(tmp.vx_m0,this->vx_m0,rhs.vx_m0,this->u_ht, this->u_len,rhs.u_len);
			return tmp;}

		vt operator*(cvta rhs)
		{	enforceLen(rhs.ui_size);
			tVecObj tmp(this->u_ht);
			tmp.enforceSize();
			MFmulV(tmp.vx_vect,this->vx_m0,rhs.vx_vect,this->u_ht, this->u_len);
			return tmp;}

		// operator*(cvta VX, cmta MY) is member of fVecObj, see below

		mt operator*(cta C)
		{	tMatObj tmp(this->u_ht, this->u_len);
			tmp.enforceSize();
			VF_mulC(tmp.vx_m0,this->vx_m0,tmp.ui_size,C);
			return tmp;}

		mt operator/(cta C)
		{	tMatObj tmp(this->u_ht, this->u_len);
			tmp.enforceSize();
			VF_divC(tmp.vx_m0,this->vx_m0,tmp.ui_size,C);
			return tmp;}

		friend mt operator+(cta C, cmta X)
		{	tMatObj tmp(X.u_ht,X.u_len);
			VF_addC(tmp.vx_m0,X.vx_m0,tmp.ui_size,C);
			return tmp;}

		friend mt operator-(cta C, cmta X)
		{	tMatObj tmp(X.u_ht,X.u_len);
			VF_subrC(tmp.vx_m0,X.vx_m0,tmp.ui_size,C);
			return tmp;}

		friend mt operator*(cta C, cmta X)
		{	tMatObj tmp(X.u_ht,X.u_len);
			VF_mulC(tmp.vx_m0,X.vx_m0,tmp.ui_size,C);
			return tmp;}

/////////////////////////////////////////////////////////////////////////////////////////
/////////////////////////////////////////////////////////////////////////////////////////
/////////////////////////////////////////////////////////////////////////////////////////
//OptiVec functions in member format
//from MXstd.h
//

/////////////////////////////////////////////////////////////////////////////////////////
//Generation
//
		//these are not recommended, better to use constructors
		void _matrix(cuint ht, cuint len) //underscore necessary to avoid confusing with constructor
		{
			dealloc();
			ui_size=((ui)ht)*len;
			if(ui_size){ u_ht=ht;u_len=len;mx_matr= MF_matrix(ht,len);vx_m0=mx_matr[0];}
		}

		void _matrix0(cuint ht, cuint len) //underscore necessary to avoid confusing with constructor
		{
			dealloc();
			ui_size=((ui)ht)*len;
			if(ui_size){ u_ht=ht;u_len=len;mx_matr= MF_matrix0(ht,len);vx_m0=mx_matr[0];}
		}


/////////////////////////////////////////////////////////////////////////////////////////
//Addressing single matrix elements
//
		t* Pelement(cuint m, cuint n) const
		{
			if(u_ht<=m || u_len<=n) throw OptiVec_SmallSizeErr();
			return MFPelement(vx_m0,u_ht,u_len,m,n);
		}

		t element(cuint m, cuint n) const
		{
			if(u_ht<=m || u_len<=n) throw OptiVec_SmallSizeErr();
			return MFelement(vx_m0,u_ht,u_len,m,n);
		}


/////////////////////////////////////////////////////////////////////////////////////////
//Initialization
//

		void equ0()
		{
			enforceSize();
			VF_equ0(vx_m0,ui_size);
		}

		void equ1()
		{
			enforceSize();
			MFequ1(vx_m0,u_len);
		}

		void equm1()
		{
			enforceSize();
			MFequm1(vx_m0,u_len);
		}

		long random(const long& seed, cta Min, cta Max)
		{
			enforceSize();
			return VF_random(vx_m0, ui_size, seed, Min, Max);
		}

		void outerprod(cvta X, cvta Y)
		{
			enforceSize( X.ui_size, Y.ui_size);
			MFouterprod(vx_m0, X.vx_vect, Y.vx_vect, u_ht, u_len);
		}

		void Row_equ0(cuint iRow)
		{
			if( !u_len || (u_ht <= iRow) )throw OptiVec_SmallSizeErr();
			MFRow_equ0(vx_m0, u_ht, u_len, iRow);
		}

		void Col_equ0(cuint iCol)
		{
			if( !u_ht || (u_len <= iCol) )throw OptiVec_SmallSizeErr();
			MFCol_equ0(vx_m0, u_ht, u_len, iCol);
		}

		void Dia_equ0()
		{
			if( !u_len || (u_len != u_ht) )throw OptiVec_SmallSizeErr();
			MFDia_equ0(vx_m0, u_len);
		}

		void Row_equC(cuint iRow, cta C)
		{
			if( !u_len || (u_ht <= iRow) )throw OptiVec_SmallSizeErr();
			MFRow_equC(vx_m0, u_ht, u_len, iRow, C);
		}

		void Col_equC(cuint iCol, cta C)
		{
			if( !u_ht || (u_len <= iCol) )throw OptiVec_SmallSizeErr();
			MFCol_equC(vx_m0, u_ht, u_len, iCol, C);
		}

		void Dia_equC(cta C)
		{
			if( !u_len || (u_len != u_ht) )throw OptiVec_SmallSizeErr();
			MFDia_equC(vx_m0, u_len, C);
		}

		void Row_equV(cuint iRow, cvta X)
		{
			if( !u_len || (u_ht <= iRow) )throw OptiVec_SmallSizeErr();
			MFRow_equV(vx_m0, u_ht, u_len, iRow, X.vx_vect);
		}

		void Col_equV(cuint iCol, cvta X)
		{
			if( !u_ht || (u_len <= iCol) )throw OptiVec_SmallSizeErr();
			MFCol_equV(vx_m0, u_ht, u_len, iCol, X.vx_vect);
		}

		void Dia_equV(cvta X)
		{
			if( !u_len || (u_len != u_ht) )throw OptiVec_SmallSizeErr();
			MFDia_equV(vx_m0, u_len, X.vx_vect);
		}

		void Trd_equM(cmta MTrd)
		{
			if( MTrd.u_ht != 3 || u_ht != u_len || u_len != MTrd.u_len )throw OptiVec_SizeErr();
			MFTrd_equM(vx_m0, MTrd.vx_m0, u_len);
		}

		void Trd_extract(cmta MA)
		{
			if( u_ht != 3 || MA.u_ht != MA.u_len || MA.u_len != u_len )throw OptiVec_SizeErr();
			MFTrd_extract(vx_m0, MA.vx_m0, u_len);
		}

		void equM(cmta X)
		{
			enforceSize(X.u_ht, X.u_len);
			MFequM(vx_m0, X.vx_m0, u_ht, u_len);
		}

		void neg(cmta X)
		{
			enforceSize(X.u_ht, X.u_len);
			MFneg(vx_m0, X.vx_m0, u_ht, u_len);
		}

		void UequL()
		{
			if( !u_len || u_len != u_ht )throw OptiVec_SizeErr();
			MFUequL(vx_m0, u_len);
		}

		void LequU()
		{
			if( !u_len || u_len != u_ht )throw OptiVec_SizeErr();
			MFLequU(vx_m0, u_len);
		}

		void Parzen()
		{
			enforceSize();
			MFParzen(vx_m0, u_ht, u_len);
		}

		void Hanning()
		{
			enforceSize();
			MFHanning(vx_m0, u_ht, u_len);
		}

		void Welch()
		{
			enforceSize();
			MFWelch(vx_m0, u_ht, u_len);
		}

/////////////////////////////////////////////////////////////////////////////////////////
//Data-type interconversions
//
		//defined by including dMatObj.h
		void DtoF(const dMatObj& X);

		//defined by including eMatObj.h
		void EtoF(const eMatObj& X);



/////////////////////////////////////////////////////////////////////////////////////////
//Submatrices
//
		void submatrix(cmta MSrce, cuint firstRowInCol, cuint sampInCol, cuint firstColInRow, cuint sampInRow )
		{
			if( (MSrce.u_ht  <= firstRowInCol+(u_ht-1) *sampInCol) ||
			   (MSrce.u_len <= firstColInRow+(u_len-1)*sampInRow) ) throw OptiVec_SizeErr();
			MFsubmatrix(vx_m0, u_ht, u_len, MSrce.vx_m0, MSrce.u_ht, MSrce.u_len,
			      firstRowInCol, sampInCol, firstColInRow, sampInRow );
		}

		void submatrix_equM(cuint firstRowInCol, cuint sampInCol, cuint firstColInRow, cuint sampInRow, cmta MSrce)
		{
			if( (u_ht  <= firstRowInCol+(MSrce.u_ht-1) *sampInCol) ||
			    (u_len <= firstColInRow+(MSrce.u_len-1)*sampInRow) ) throw OptiVec_SizeErr();
			MFsubmatrix_equM(vx_m0, u_ht, u_len,
			    firstRowInCol, sampInCol, firstColInRow, sampInRow,
			    MSrce.vx_m0, MSrce.u_ht, MSrce.u_len );
		}

		void equMblock(cmta MSrce, cuint firstRowInCol, cuint firstColInRow )
		{
			if( (MSrce.u_ht  < firstRowInCol+u_ht) ||
			    (MSrce.u_len < firstColInRow+u_len) ) throw OptiVec_SizeErr();
			MFequMblock(vx_m0, u_ht, u_len, MSrce.vx_m0, MSrce.u_ht, MSrce.u_len,
			            firstRowInCol, firstColInRow );
		}

		void equMblockT(cmta MSrce, cuint firstRowInCol, cuint firstColInRow )
		{
			if( (MSrce.u_ht  < firstRowInCol+u_len) ||
			    (MSrce.u_len < firstColInRow+u_ht) ) throw OptiVec_SizeErr();
			MFequMblockT(vx_m0, u_ht, u_len, MSrce.vx_m0, MSrce.u_ht, MSrce.u_len,
			             firstRowInCol, firstColInRow );
		}

		void block_equM(cuint firstRowInCol, cuint firstColInRow, cmta MSrce )
		{
			if( (u_ht  < firstRowInCol+MSrce.u_ht) ||
			     (u_len < firstColInRow+MSrce.u_len) ) throw OptiVec_SizeErr();
			MFblock_equM(vx_m0, u_ht, u_len, firstRowInCol, firstColInRow,
			             MSrce.vx_m0, MSrce.u_ht, MSrce.u_len );
		}

		void block_equMT(cuint firstRowInCol, cuint firstColInRow, cmta MSrce )
		{
			if( (u_ht  < firstRowInCol+MSrce.u_len) ||
			    (u_len < firstColInRow+MSrce.u_ht) ) throw OptiVec_SizeErr();
			MFblock_equMT(vx_m0, u_ht, u_len, firstRowInCol, firstColInRow,
			              MSrce.vx_m0, MSrce.u_ht, MSrce.u_len );
		}

/////////////////////////////////////////////////////////////////////////////////////////
// Deleting or inserting a row or a column
//

		void Row_delete(cmta MA, cuint iRow)
		{
			if( MA.u_ht <= iRow ) throw OptiVec_SmallSizeErr();
			if( mx_matr == MA.mx_matr ) // overwrite MA?
			{	ui_size -= u_len;  // don't reallocate, as row pointers remain intact
				MFRow_delete(vx_m0, vx_m0, u_ht--, u_len, iRow);
			}
			else
			{	enforceSize(MA.u_ht-1, MA.u_len);
				MFRow_delete(vx_m0, MA.vx_m0, MA.u_ht, MA.u_len, iRow);
			}
		}

		void Col_delete(cmta MA, cuint iCol)
		{
			if( MA.u_len <= iCol ) throw OptiVec_SmallSizeErr();
			if( mx_matr == MA.mx_matr ) // overwrite MA?
			{
				tMatrix tmp = MF_matrix( u_ht, u_len-1 );
				MFCol_delete(tmp[0], MA.vx_m0, u_ht, u_len--, iCol);
				M_free( mx_matr );
				vx_m0 = tmp[0];
				ui_size -= u_ht;
				mx_matr = tmp;
			}
			else
			{
				enforceSize(MA.u_ht, MA.u_len-1);
				MFCol_delete(vx_m0, MA.vx_m0, u_ht, u_len, iCol);
			}
		}

		void Row_insert(cmta MA, cuint iRow, cvta X)
		{
			if( u_ht <= iRow ) throw OptiVec_SmallSizeErr();
			if( mx_matr == MA.mx_matr ) // overwrite MA?
			{
				tMatrix tmp = MF_matrix( ++u_ht, u_len );
				MFRow_insert(tmp[0], MA.vx_m0, u_ht, u_len, iRow, X.vx_vect);
				M_free( mx_matr );
				vx_m0 = tmp[0];
				ui_size += u_len;
				mx_matr = tmp;
			}
			else
			{
				enforceSize(MA.u_ht+1, MA.u_len);
				MFRow_insert(vx_m0, MA.vx_m0, u_ht, u_len, iRow, X.vx_vect);
			}
		}

		void Col_insert(cmta MA, cuint iCol, cvta X)
		{
			if( u_len <= iCol ) throw OptiVec_SmallSizeErr();
			if( mx_matr == MA.mx_matr ) // overwrite MA?
			{
				tMatrix tmp = MF_matrix( u_ht, ++u_len );
				MFCol_insert(tmp[0], MA.vx_m0, u_ht, u_len, iCol, X.vx_vect);
				M_free( mx_matr );
				vx_m0 = tmp[0];
				ui_size += u_ht;
				mx_matr = tmp;
			}
			else
			{
				enforceSize(MA.u_ht, MA.u_len+1);
				MFCol_insert(vx_m0, MA.vx_m0, u_ht, u_len, iCol, X.vx_vect);
			}
		}


/////////////////////////////////////////////////////////////////////////////////////////
//Basic arithmetic operations performed on one single row, or one single column
// of any matrix, or on the diagonal of a square matrix
//

		void Row_addC(cuint iRow, cta C)
		{
			if( !u_len || (u_ht <= iRow) )throw OptiVec_SmallSizeErr();
			MFRow_addC(vx_m0, u_ht, u_len, iRow, C);
		}

		void Col_addC(cuint iCol, cta C)
		{
			if( !u_ht || (u_len <= iCol) )throw OptiVec_SmallSizeErr();
			MFCol_addC(vx_m0, u_ht, u_len, iCol, C);
		}

		void Dia_addC( cta C)
		{
			if( !u_len || (u_len != u_ht) )throw OptiVec_SmallSizeErr();
			MFDia_addC(vx_m0, u_len, C);
		}

		void Row_addV(cuint iRow, cvta X)
		{
			if( !u_len || (u_ht <= iRow) || (X.ui_size < u_len))throw OptiVec_SmallSizeErr();
			MFRow_addV(vx_m0, u_ht, u_len, iRow, X.vx_vect);
		}

		void Col_addV(cuint iCol, cvta X)
		{
			if( !u_ht || (u_len <= iCol) || (X.ui_size < u_ht) )throw OptiVec_SmallSizeErr();
			MFCol_addV(vx_m0, u_ht, u_len, iCol, X.vx_vect);
		}

		void Dia_addV( cvta X)
		{
			if( !u_len || (u_len != u_ht) || (X.ui_size < u_len) )throw OptiVec_SmallSizeErr();
			MFDia_addV(vx_m0, u_len, X.vx_vect);
		}

		void Row_subC(cuint iRow, cta C)
		{
			if( !u_len || (u_ht <= iRow) )throw OptiVec_SmallSizeErr();
			MFRow_subC(vx_m0, u_ht, u_len, iRow, C);
		}

		void Col_subC(cuint iCol, cta C)
		{
			if( !u_ht || (u_len <= iCol) )throw OptiVec_SmallSizeErr();
			MFCol_subC(vx_m0, u_ht, u_len, iCol, C);
		}

		void Dia_subC( cta C)
		{
			if( !u_len || (u_len != u_ht) )throw OptiVec_SmallSizeErr();
			MFDia_subC(vx_m0, u_len, C);
		}

		void Row_subV(cuint iRow, cvta X)
		{
			if( !u_len || (u_ht <= iRow) || (X.ui_size < u_len))throw OptiVec_SmallSizeErr();
			MFRow_subV(vx_m0, u_ht, u_len, iRow, X.vx_vect);
		}

		void Col_subV(cuint iCol, cvta X)
		{
			if( !u_ht || (u_len <= iCol) || (X.ui_size < u_ht) )throw OptiVec_SmallSizeErr();
			MFCol_subV(vx_m0, u_ht, u_len, iCol, X.vx_vect);
		}

		void Dia_subV( cvta X)
		{
			if( !u_len || (u_len != u_ht) || (X.ui_size < u_len) )throw OptiVec_SmallSizeErr();
			MFDia_subV(vx_m0, u_len, X.vx_vect);
		}

		void Row_subrC(cuint iRow, cta C)
		{
			if( !u_len || (u_ht <= iRow) )throw OptiVec_SmallSizeErr();
			MFRow_subrC(vx_m0, u_ht, u_len, iRow, C);
		}

		void Col_subrC(cuint iCol, cta C)
		{
			if( !u_ht || (u_len <= iCol) )throw OptiVec_SmallSizeErr();
			MFCol_subrC(vx_m0, u_ht, u_len, iCol, C);
		}

		void Dia_subrC( cta C)
		{
			if( !u_len || (u_len != u_ht) )throw OptiVec_SmallSizeErr();
			MFDia_subrC(vx_m0, u_len, C);
		}

		void Row_subrV(cuint iRow, cvta X)
		{
			if( !u_len || (u_ht <= iRow) || (X.ui_size < u_len))throw OptiVec_SmallSizeErr();
			MFRow_subrV(vx_m0, u_ht, u_len, iRow, X.vx_vect);
		}

		void Col_subrV(cuint iCol, cvta X)
		{
			if( !u_ht || (u_len <= iCol) || (X.ui_size < u_ht) )throw OptiVec_SmallSizeErr();
			MFCol_subrV(vx_m0, u_ht, u_len, iCol, X.vx_vect);
		}

		void Dia_subrV( cvta X)
		{
			if( !u_len || (u_len != u_ht) || (X.ui_size < u_len) )throw OptiVec_SmallSizeErr();
			MFDia_subrV(vx_m0, u_len, X.vx_vect);
		}

		void Row_mulC(cuint iRow, cta C)
		{
			if( !u_len || (u_ht <= iRow) )throw OptiVec_SmallSizeErr();
			MFRow_mulC(vx_m0, u_ht, u_len, iRow, C);
		}

		void Col_mulC(cuint iCol, cta C)
		{
			if( !u_ht || (u_len <= iCol) )throw OptiVec_SmallSizeErr();
			MFCol_mulC(vx_m0, u_ht, u_len, iCol, C);
		}

		void Dia_mulC( cta C)
		{
			if( !u_len || (u_len != u_ht) )throw OptiVec_SmallSizeErr();
			MFDia_mulC(vx_m0, u_len, C);
		}

		void Row_mulV(cuint iRow, cvta X)
		{
			if( !u_len || (u_ht <= iRow) || (X.ui_size < u_len))throw OptiVec_SmallSizeErr();
			MFRow_mulV(vx_m0, u_ht, u_len, iRow, X.vx_vect);
		}

		void Col_mulV(cuint iCol, cvta X)
		{
			if( !u_ht || (u_len <= iCol) || (X.ui_size < u_ht) )throw OptiVec_SmallSizeErr();
			MFCol_mulV(vx_m0, u_ht, u_len, iCol, X.vx_vect);
		}

		void Dia_mulV( cvta X)
		{
			if( !u_len || (u_len != u_ht) || (X.ui_size < u_len) )throw OptiVec_SmallSizeErr();
			MFDia_mulV(vx_m0, u_len, X.vx_vect);
		}

		void Row_divC(cuint iRow, cta C)
		{
			if( !u_len || (u_ht <= iRow) )throw OptiVec_SmallSizeErr();
			MFRow_divC(vx_m0, u_ht, u_len, iRow, C);
		}

		void Col_divC(cuint iCol, cta C)
		{
			if( !u_ht || (u_len <= iCol) )throw OptiVec_SmallSizeErr();
			MFCol_divC(vx_m0, u_ht, u_len, iCol, C);
		}

		void Dia_divC( cta C)
		{
			if( !u_len || (u_len != u_ht) )throw OptiVec_SmallSizeErr();
			MFDia_divC(vx_m0, u_len, C);
		}

		void Row_divV(cuint iRow, cvta X)
		{
			if( !u_len || (u_ht <= iRow) || (X.ui_size < u_len))throw OptiVec_SmallSizeErr();
			MFRow_divV(vx_m0, u_ht, u_len, iRow, X.vx_vect);
		}

		void Col_divV(cuint iCol, cvta X)
		{
			if( !u_ht || (u_len <= iCol) || (X.ui_size < u_ht) )throw OptiVec_SmallSizeErr();
			MFCol_divV(vx_m0, u_ht, u_len, iCol, X.vx_vect);
		}

		void Dia_divV( cvta X)
		{
			if( !u_len || (u_len != u_ht) || (X.ui_size < u_len) )throw OptiVec_SmallSizeErr();
			MFDia_divV(vx_m0, u_len, X.vx_vect);
		}

		void Row_divrC(cuint iRow, cta C)
		{
			if( !u_len || (u_ht <= iRow) )throw OptiVec_SmallSizeErr();
			MFRow_divrC(vx_m0, u_ht, u_len, iRow, C);
		}

		void Col_divrC(cuint iCol, cta C)
		{
			if( !u_ht || (u_len <= iCol) )throw OptiVec_SmallSizeErr();
			MFCol_divrC(vx_m0, u_ht, u_len, iCol, C);
		}

		void Dia_divrC( cta C)
		{
			if( !u_len || (u_len != u_ht) )throw OptiVec_SmallSizeErr();
			MFDia_divrC(vx_m0, u_len, C);
		}

		void Row_divrV(cuint iRow, cvta X)
		{
			if( !u_len || (u_ht <= iRow) || (X.ui_size < u_len))throw OptiVec_SmallSizeErr();
			MFRow_divrV(vx_m0, u_ht, u_len, iRow, X.vx_vect);
		}

		void Col_divrV(cuint iCol, cvta X)
		{
			if( !u_ht || (u_len <= iCol) || (X.ui_size < u_ht) )throw OptiVec_SmallSizeErr();
			MFCol_divrV(vx_m0, u_ht, u_len, iCol, X.vx_vect);
		}

		void Dia_divrV( cvta X)
		{
			if( !u_len || (u_len != u_ht) || (X.ui_size < u_len) )throw OptiVec_SmallSizeErr();
			MFDia_divrV(vx_m0, u_len, X.vx_vect);
		}


/////////////////////////////////////////////////////////////////////////////////////////
//One-dimensional vector operations
//performed along all rows or columns, or along the diagonal of a square matrix
//(some of these functions belong to fVecObj.h)

		t Dia_max() const
		{
			if( !u_len || (u_len != u_ht) )throw OptiVec_SizeErr();
			return MFDia_max(vx_m0, u_len);
		}

		t Dia_min() const
		{
			if( !u_len || (u_len != u_ht) )throw OptiVec_SizeErr();
			return MFDia_min(vx_m0, u_len);
		}

		t Dia_absmax() const
		{
			if( !u_len || (u_len != u_ht) )throw OptiVec_SizeErr();
			return MFDia_absmax(vx_m0, u_len);
		}

		t Dia_absmin() const
		{
			if( !u_len || (u_len != u_ht) )throw OptiVec_SizeErr();
			return MFDia_absmin(vx_m0, u_len);
		}

		t Dia_sum() const
		{
			if( !u_len || (u_len != u_ht) )throw OptiVec_SizeErr();
			return MFDia_sum(vx_m0, u_len);
		}

		t Dia_prod() const
		{
			if( !u_len || (u_len != u_ht) )throw OptiVec_SizeErr();
			return MFDia_prod(vx_m0, u_len);
		}

		void Rows_runsum()
		{
			enforceSize();
			MFRows_runsum(vx_m0, u_ht, u_len);
		}

		void Cols_runsum()
		{
			enforceSize();
			MFCols_runsum(vx_m0, u_ht, u_len);
		}

		void Rows_runprod()
		{
			enforceSize();
			MFRows_runprod(vx_m0, u_ht, u_len);
		}

		void Cols_runprod()
		{
			enforceSize();
			MFCols_runprod(vx_m0, u_ht, u_len);
		}

		void Rows_rotate(const int pos)
		{
			enforceSize();
			MFRows_rotate(vx_m0, u_ht, u_len, pos);
		}

		void Cols_rotate(const int pos)
		{
			enforceSize();
			MFCols_rotate(vx_m0, u_ht, u_len, pos);
		}

		void Rows_rotate_buf(const int pos, cmta MBuf)
		{
			enforceSize();
			MFRows_rotate_buf(vx_m0, u_ht, u_len, pos, MBuf.vx_m0);
		}

		void Cols_rotate_buf(const int pos, cmta MBuf)
		{
			enforceSize();
			MFCols_rotate_buf(vx_m0, u_ht, u_len, pos, MBuf.vx_m0);
		}

		void Rows_reflect()
		{
			enforceSize();
			MFRows_reflect(vx_m0, u_ht, u_len);
		}

		void Cols_reflect()
		{
			enforceSize();
			MFCols_reflect(vx_m0, u_ht, u_len);
		}


/////////////////////////////////////////////////////////////////////////////////////////
//Operations involving two rows or two colums of one matrix
//
		void Rows_exchange(cuint row1, cuint row2)
		{
			if( !u_len || (u_ht <= row1) || (u_ht <= row2))throw OptiVec_SizeErr();
			MFRows_exchange(vx_m0, u_ht, u_len, row1, row2);
		}

		void Cols_exchange(cuint col1, cuint col2)
		{
			if( !u_len || (u_ht <= col1) || (u_ht <= col2))throw OptiVec_SizeErr();
			MFCols_exchange(vx_m0, u_ht, u_len, col1, col2);
		}

		void Rows_add(cuint destRow, cuint srceRow)
		{
			if( !u_len || (u_ht <= destRow) || (u_ht <= srceRow))throw OptiVec_SizeErr();
			MFRows_add(vx_m0, u_ht, u_len, destRow, srceRow);
		}

		void Cols_add(cuint destCol, cuint srceCol)
		{
			if( !u_len || (u_ht <= destCol) || (u_ht <= srceCol))throw OptiVec_SizeErr();
			MFCols_add(vx_m0, u_ht, u_len, destCol, srceCol);
		}

		void Rows_sub(cuint destRow, cuint srceRow)
		{
			if( !u_len || (u_ht <= destRow) || (u_ht <= srceRow))throw OptiVec_SizeErr();
			MFRows_sub(vx_m0, u_ht, u_len, destRow, srceRow);
		}

		void Cols_sub(cuint destCol, cuint srceCol)
		{
			if( !u_len || (u_ht <= destCol) || (u_ht <= srceCol))throw OptiVec_SizeErr();
			MFCols_sub(vx_m0, u_ht, u_len, destCol, srceCol);
		}

		void Rows_Cadd(cuint destRow, cuint srceRow, cta C)
		{
			if( !u_len || (u_ht <= destRow) || (u_ht <= srceRow))throw OptiVec_SizeErr();
			MFRows_Cadd(vx_m0, u_ht, u_len, destRow, srceRow, C);
		}

		void Cols_Cadd(cuint destCol, cuint srceCol, cta C)
		{
			if( !u_len || (u_ht <= destCol) || (u_ht <= srceCol))throw OptiVec_SizeErr();
			MFCols_Cadd(vx_m0, u_ht, u_len, destCol, srceCol, C);
		}

		void Rows_lincomb(cuint destRow, cta destC, cuint srceRow, cta srceC)
		{
			if( !u_len || (u_ht <= destRow) || (u_ht <= srceRow))throw OptiVec_SizeErr();
			MFRows_lincomb(vx_m0, u_ht, u_len, destRow, destC, srceRow, srceC);
		}

		void Cols_lincomb(cuint destCol, cta destC, cuint srceCol, cta srceC)
		{
			if( !u_len || (u_ht <= destCol) || (u_ht <= srceCol)) throw OptiVec_SizeErr();
			MFCols_lincomb(vx_m0, u_ht, u_len, destCol, destC, srceCol, srceC);
		}

/////////////////////////////////////////////////////////////////////////////////////////
//Center of Gravity
//

		fComplex centerOfGravityInd() const
		{
			enforceSize();
			return MFcenterOfGravityInd(vx_m0,u_ht, u_len);
		}

		fComplex centerOfGravityV(cvta X, cvta Y) const
		{
            if( u_ht != Y.ui_size || u_len != X.ui_size || !ui_size )
                throw OptiVec_SizeErr();
			return MFcenterOfGravityV(X.vx_vect,Y.vx_vect,vx_m0,u_ht,u_len);
		}

/////////////////////////////////////////////////////////////////////////////////////////
//Matrix transposition
//
		void transpose(cmta MA)
		{
			enforceSize(MA.u_len, MA.u_ht);
			MFtranspose(vx_m0, MA.vx_m0, u_ht, u_len);
		}

/////////////////////////////////////////////////////////////////////////////////////////
//Interpolation of Matrix Elements
//
		void polyinterpol(cvta X, cvta Y, cvta Xtab, cvta Ytab, cmta MZtab,
                          unsigned degX, unsigned degY)
		{
            if( u_ht != Y.ui_size || u_len != X.ui_size || !ui_size ||
                MZtab.u_ht != Ytab.ui_size || MZtab.u_len != Xtab.ui_size
                || !MZtab.ui_size )    throw OptiVec_SizeErr();
			if((degX>10) || (degY>10)) throw OptiVec_Err();
			if((Xtab.ui_size<degX) || (Ytab.ui_size<degY)) throw OptiVec_SmallSizeErr();
			MFpolyinterpol(vx_m0,X.vx_vect,Y.vx_vect,u_ht,u_len,
				Xtab.vx_vect,Ytab.vx_vect,MZtab.vx_m0,MZtab.u_ht,MZtab.u_len,degX,degY);
		}

		void ratinterpol(cvta X, cvta Y, cvta Xtab, cvta Ytab, cmta MZtab,
                          unsigned degX, unsigned degY)
		{
            if( u_ht != Y.ui_size || u_len != X.ui_size || !ui_size ||
                MZtab.u_ht != Ytab.ui_size || MZtab.u_len != Xtab.ui_size
                || !MZtab.ui_size )    throw OptiVec_SizeErr();
			if((degX>10) || (degY>10)) throw OptiVec_Err();
			if((Xtab.ui_size<degX) || (Ytab.ui_size<degY)) throw OptiVec_SmallSizeErr();
			MFratinterpol(vx_m0,X.vx_vect,Y.vx_vect,u_ht,u_len,
				Xtab.vx_vect,Ytab.vx_vect,MZtab.vx_m0,MZtab.u_ht,MZtab.u_len,degX,degY);
		}

		void natCubSplineInterpol(cvta X, cvta Y, cvta Xtab, cvta Ytab, cmta MZtab)
		{
            if( u_ht != Y.ui_size || u_len != X.ui_size || !ui_size ||
                MZtab.u_ht != Ytab.ui_size || MZtab.u_len != Xtab.ui_size
                || !MZtab.ui_size )    throw OptiVec_SizeErr();
			if((Xtab.ui_size<3) || (Ytab.ui_size<3)) throw OptiVec_SmallSizeErr();
			MFnatCubSplineInterpol(vx_m0,X.vx_vect,Y.vx_vect,u_ht,u_len,
				Xtab.vx_vect,Ytab.vx_vect,MZtab.vx_m0,MZtab.u_ht,MZtab.u_len);
		}

/////////////////////////////////////////////////////////////////////////////////////////
//Matrix arithmetics
//(The following functions are more efficient than the arithmetic operators defined above)
//
		void addM(cmta MA, cmta MB)
		{
			enforceSize(MA.u_ht, MA.u_len);
			enforceSize(MB.u_ht, MB.u_len);
			MFaddM(vx_m0, MA.vx_m0, MB.vx_m0, u_ht, u_len);
		}

		void addMT(cmta MA, cmta MB)
		{
			enforceSize(MA.u_ht,  MA.u_len);
			enforceSize(MB.u_len, MB.u_ht);
			MFaddMT(vx_m0, MA.vx_m0, MB.vx_m0, u_ht, u_len);
		}

		void subM(cmta MA, cmta MB)
		{
			enforceSize(MA.u_ht, MA.u_len);
			enforceSize(MB.u_ht, MB.u_len);
			MFsubM(vx_m0, MA.vx_m0, MB.vx_m0, u_ht, u_len);
		}

		void subMT(cmta MA, cmta MB)
		{
			enforceSize(MA.u_ht,  MA.u_len);
			enforceSize(MB.u_len, MB.u_ht);
			MFsubMT(vx_m0, MA.vx_m0, MB.vx_m0, u_ht, u_len);
		}

		void subrM(cmta MA, cmta MB)
		{
			enforceSize(MA.u_ht, MA.u_len);
			enforceSize(MB.u_ht, MB.u_len);
			MFsubM(vx_m0, MB.vx_m0, MA.vx_m0, u_ht, u_len);
		}

		void subrMT(cmta MA, cmta MB)
		{
			enforceSize(MA.u_ht,  MA.u_len);
			enforceSize(MB.u_len, MB.u_ht);
			MFsubrMT(vx_m0, MA.vx_m0, MB.vx_m0, u_ht, u_len);
		}

		void mulC(cmta MA, cta C)
		{
			enforceSize(MA.u_ht, MA.u_len);
			MFmulC(vx_m0, MA.vx_m0, u_ht, u_len, C);
		}

		void divC(cmta MA, cta C)
		{
			enforceSize(MA.u_ht, MA.u_len);
			MFdivC(vx_m0, MA.vx_m0, u_ht, u_len, C);
		}

		void s_addM(cmta MA, cmta MB, cta C)
		{
			enforceSize(MA.u_ht, MA.u_len);
			enforceSize(MB.u_ht, MB.u_len);
			MFsaddM(vx_m0, MA.vx_m0, MB.vx_m0, u_ht, u_len, C);
		}

		void s_addMT(cmta MA, cmta MB, cta C)
		{
			enforceSize(MA.u_ht, MA.u_len);
			enforceSize(MB.u_len, MB.u_ht);
			MFsaddMT(vx_m0, MA.vx_m0, MB.vx_m0, u_ht, u_len, C);
		}

		void s_subM(cmta MA, cmta MB, cta C)
		{
			enforceSize(MA.u_ht, MA.u_len);
			enforceSize(MB.u_ht, MB.u_len);
			MFssubM(vx_m0, MA.vx_m0, MB.vx_m0, u_ht, u_len, C);
		}

		void s_subMT(cmta MA, cmta MB, cta C)
		{
			enforceSize(MA.u_ht, MA.u_len);
			enforceSize(MB.u_len, MB.u_ht);
			MFssubMT(vx_m0, MA.vx_m0, MB.vx_m0, u_ht, u_len, C);
		}

		void s_subrM(cmta MA, cmta MB, cta C)
		{
			enforceSize(MA.u_ht, MA.u_len);
			enforceSize(MB.u_ht, MB.u_len);
			MFssubM(vx_m0, MB.vx_m0, MA.vx_m0, u_ht, u_len, C);
		}

		void s_subrMT(cmta MA, cmta MB, cta C)
		{
			enforceSize(MA.u_ht, MA.u_len);
			enforceSize(MB.u_len, MB.u_ht);
			MFssubrMT(vx_m0, MA.vx_m0, MB.vx_m0, u_ht, u_len, C);
		}

		void lincomb(cmta MA, cmta MB, cta CA, cta CB)
		{
			enforceSize(MA.u_ht, MA.u_len);
			enforceSize(MB.u_ht, MB.u_len);
			MFlincomb(vx_m0, MA.vx_m0, MB.vx_m0, u_ht, u_len, CA, CB);
		}

		//mixed multiplications matrix * vector belong to fVecObj
		void mulM( cmta MA, cmta MB )
		{
			if( u_ht != MA.u_ht || u_len != MB.u_len || MB.u_ht != MA.u_len )throw OptiVec_SizeErr();
			MFmulM( vx_m0, MA.vx_m0, MB.vx_m0, MA.u_ht, MA.u_len, MB.u_len);
		}

		void mulMT( cmta MA, cmta MB )
		{
			if( u_ht != MA.u_ht || u_len != MB.u_ht || MB.u_len != MA.u_len )throw OptiVec_SizeErr();
			MFmulMT( vx_m0, MA.vx_m0, MB.vx_m0, MA.u_ht, MA.u_len, MB.u_ht);
		}

		void TmulM( cmta MA, cmta MB )
		{
			if( u_ht != MA.u_len || u_len != MB.u_len || MB.u_ht != MA.u_ht )throw OptiVec_SizeErr();
			MFTmulM( vx_m0, MA.vx_m0, MB.vx_m0, MA.u_ht, MA.u_len, MB.u_len);
		}

		void TmulMT( cmta MA, cmta MB )
		{
			if( u_ht != MA.u_len || u_len != MB.u_ht || MB.u_len != MA.u_ht )throw OptiVec_SizeErr();
			MFTmulMT( vx_m0, MA.vx_m0, MB.vx_m0, MA.u_ht, MA.u_len, MB.u_ht);
		}

		void dia_mulM( cvta MADia, cmta MB )
		{
			if( u_ht != MADia.ui_size || u_len != MB.u_len || MB.u_ht != MADia.ui_size )throw OptiVec_SizeErr();
			MFdiamulM( vx_m0, MADia.vx_vect, MB.vx_m0, MB.u_ht, MB.u_len );
		}

		void dia_mulMT( cvta MADia, cmta MB )
		{
			if( u_ht != MADia.ui_size || u_len != MB.u_ht || MB.u_len != MADia.ui_size )throw OptiVec_SizeErr();
			MFdiamulMT( vx_m0, MADia.vx_vect, MB.vx_m0, MB.u_ht, MB.u_len );
		}

		void mulMdia( cmta MA, cvta MBDia )
		{
			if( u_ht != MA.u_ht || u_len != MBDia.ui_size || MBDia.ui_size != MA.u_len )throw OptiVec_SizeErr();
			MFmulMdia( vx_m0, MA.vx_m0, MBDia.vx_vect, MA.u_ht, MA.u_len );
		}

		void TmulMdia( cmta MA, cvta MBDia )
		{
			if( u_ht != MA.u_len || u_len != MBDia.ui_size || MBDia.ui_size != MA.u_ht )throw OptiVec_SizeErr();
			MFTmulMdia( vx_m0, MA.vx_m0, MBDia.vx_vect, MA.u_ht, MA.u_len );
		}

/////////////////////////////////////////////////////////////////////////////////////////
//Linear Algebra
//
		int LUdecompose( uVecObj Ind, cmta MA )
		{
			if( u_ht != MA.u_ht || u_len != MA.u_len || u_ht != u_len || u_ht != Ind.ui_size )
					throw OptiVec_SizeErr();
			return MFLUdecompose( vx_m0, Ind.vx_vect, MA.vx_m0, u_len );
		}

		int LUdecompose( uVecObj* Ind, cmta MA )
		{
			if( u_ht != MA.u_ht || u_len != MA.u_len || u_ht != u_len || u_ht != Ind->ui_size )
					throw OptiVec_SizeErr();
			return MFLUdecompose( vx_m0, Ind->vx_vect, MA.vx_m0, u_len );
		}

		int  LUDresult() const
		{
			return MF_LUDresult();
		}

		void LUDsetEdit(cta Thresh) const
		{
			MF_LUDsetEdit(Thresh);
		}

		t LUDgetEdit() const
		{
			return MF_LUDgetEdit();
		}

		void LUinv( cmta MLU, const uVecObj& Ind )
		{
			if( u_ht != MLU.u_ht || u_len != MLU.u_len || u_ht != u_len || u_ht != Ind.ui_size )
					throw OptiVec_SizeErr();
			MFLUinv( vx_m0, MLU.vx_m0, Ind.vx_vect, u_len );
		}

		t LUdet( int permut ) const
		{
			if( u_ht != u_len ) throw OptiVec_SizeErr();
			return MFLUdet( vx_m0, u_len, permut );
		}

		int SVdecompose( tMatObj *MV, tVecObj *W, cmta MA )
		{
			if( (MA.u_ht < MA.u_len ? (u_ht != u_len) : (u_ht != MA.u_ht))
				|| u_len != MA.u_len || u_len != MV->u_len || u_len != MV->u_ht
				|| u_len != W->ui_size )   throw OptiVec_SizeErr();
			return MFSVdecompose( vx_m0, MV->vx_m0, W->vx_vect,
			          MA.vx_m0, MA.u_ht, MA.u_len );
		}

		int SVdecompose( tMatObj MV, tVecObj W, cmta MA )
		{
			if( (MA.u_ht < MA.u_len ? (u_ht != u_len) : (u_ht != MA.u_ht))
			  || u_len != MA.u_len || u_len != MV.u_len || u_len != MV.u_ht
			  || u_len != W.ui_size )   throw OptiVec_SizeErr();
			return MFSVdecompose( vx_m0, MV.vx_m0, W.vx_vect,
			         MA.vx_m0, MA.u_ht, MA.u_len );
		}

		void SVDsetEdit(cta Thresh) const
		{
			MF_SVDsetEdit(Thresh);
		}

		t SVDgetEdit() const
		{
			return MF_SVDgetEdit();
		}

		int inv( cmta MA )
		{
			if( u_ht != MA.u_ht || u_len != MA.u_len || u_ht != u_len )
				throw OptiVec_SizeErr();
			return MFinv( vx_m0, MA.vx_m0, u_len );
		}

		t det() const
		{
			if( u_ht != u_len ) throw OptiVec_SizeErr();
			return MFdet( vx_m0, u_len );
		}

/////////////////////////////////////////////////////////////////////////////////////////
// Eigenvalues, Matrix Square-Root
//
// sym_eigenval is a member of tVecObj

		int sym_sqrt( cmta MA )
		{
			if( u_ht != MA.u_ht || u_len != MA.u_len || u_ht != u_len )
			throw OptiVec_SizeErr();
			return MFs_sqrt( vx_m0, MA.vx_m0, u_len );
		}


/////////////////////////////////////////////////////////////////////////////////////////
//Fourier Transforms, Convolutions, Filtering
//

		tcplx getRspEdit() const
		{
			return VF_getRspEdit();
		}

		void setRspEdit(const tcplx& trunc) const
		{
			VF_setRspEdit(trunc);
		}

		void p_FFT(cmta X, int dir=1)
		{
			enforcePow2();
			MFpFFT(vx_m0,X.vx_m0,u_ht,u_len,dir);
		}

		//p_FFTtoC member of cfMatObj

		void p_filter(cmta X, cmta Flt)
		{
			enforcePow2();
			MFpfilter(vx_m0, X.vx_m0, Flt.vx_m0, u_ht, u_len);
		}

		void p_autocorr(cmta X)
		{
			enforcePow2();
			MFpautocorr(vx_m0, X.vx_m0, u_ht, u_len);
		}

		void p_xcorr(cmta X, cmta Y)
		{
			enforcePow2();
			MFpxcorr(vx_m0, X.vx_m0, Y.vx_m0, u_ht, u_len);
		}

		void p_spectrum(cmta X, cmta Win)
		{
			unsigned htSpec=u_ht-1, lenSpec=u_len-1;
			if(!isipow2(htSpec) || !isipow2(lenSpec)) throw OptiVec_Pow2SizeErr();
			if(X.u_ht < 2*htSpec || X.u_len < 2*lenSpec) throw OptiVec_BigSizeErr();
			if(Win.u_ht != 2*htSpec || Win.u_len != 2*lenSpec) throw OptiVec_SizeErr();
			MFpspectrum(vx_m0, htSpec, lenSpec, X.vx_m0, X.u_ht, X.u_len, Win.vx_m0);
		}

		void p_convolve(tMatObj* Flt, cmta X, cmta Rsp)
		{
			enforcePow2();
			if((u_ht!=Flt->u_ht)||(u_len!=Flt->u_len)||(u_ht!=X.u_ht)||(u_len!=X.u_len)||
				(u_ht!=Rsp.u_ht)||(u_len!=Rsp.u_len)) throw OptiVec_SizeErr();
			if((mx_matr==Rsp.mx_matr)||(Flt->mx_matr==X.mx_matr))
				throw OptiVec_Err();
			MFpconvolve(vx_m0, Flt->vx_m0, X.vx_m0, Rsp.vx_m0, u_ht, u_len);
		}

		 void p_convolve(tMatObj Flt, cmta X, cmta Rsp)
		{
			enforcePow2();
			if((u_ht!=Flt.u_ht)||(u_len!=Flt.u_len)||(u_ht!=X.u_ht)||(u_len!=X.u_len)||
				(u_ht!=Rsp.u_ht)||(u_len!=Rsp.u_len)) throw OptiVec_SizeErr();
			if((mx_matr==Rsp.mx_matr)||(Flt.mx_matr==X.mx_matr))
				throw OptiVec_Err();
			MFpconvolve(vx_m0, Flt.vx_m0, X.vx_m0, Rsp.vx_m0, u_ht, u_len);
		}

		void p_deconvolve(tMatObj* Flt, cmta X, cmta Rsp)
		{
			enforcePow2();
			if((u_ht!=Flt->u_ht)||(u_len!=Flt->u_len)||(u_ht!=X.u_ht)||(u_len!=X.u_len)||
				(u_ht!=Rsp.u_ht)||(u_len!=Rsp.u_len)) throw OptiVec_SizeErr();
			if((mx_matr==Rsp.mx_matr)||(Flt->mx_matr==X.mx_matr))
				throw OptiVec_Err();
			MFpdeconvolve(vx_m0, Flt->vx_m0, X.vx_m0, Rsp.vx_m0, u_ht, u_len);
		}

		void p_deconvolve(tMatObj Flt, cmta X, cmta Rsp)
		{
			enforcePow2();
			if((u_ht!=Flt.u_ht)||(u_len!=Flt.u_len)||(u_ht!=X.u_ht)||(u_len!=X.u_len)||
				(u_ht!=Rsp.u_ht)||(u_len!=Rsp.u_len)) throw OptiVec_SizeErr();
			if((mx_matr==Rsp.mx_matr)||(Flt.mx_matr==X.mx_matr))
				throw OptiVec_Err();
			MFpdeconvolve(vx_m0, Flt.vx_m0, X.vx_m0, Rsp.vx_m0, u_ht, u_len);
		}

		void l_FFT(cmta X, int dir=1)
		{
			enforcePow2();
			MFlFFT(vx_m0,X.vx_m0,u_ht,u_len,dir);
		}

		//l_FFTtoC member of cfMatObj

		void l_filter(cmta X, cmta Flt)
		{
			enforcePow2();
			MFlfilter(vx_m0, X.vx_m0, Flt.vx_m0, u_ht, u_len);
		}

		void l_autocorr(cmta X)
		{
			enforcePow2();
			MFlautocorr(vx_m0, X.vx_m0, u_ht, u_len);
		}

		void l_xcorr(cmta X, cmta Y)
		{
			enforcePow2();
			MFlxcorr(vx_m0, X.vx_m0, Y.vx_m0, u_ht, u_len);
		}

		void l_spectrum(cmta X, cmta Win)
		{
			unsigned htSpec=u_ht-1, lenSpec=u_len-1;
			if(!isipow2(htSpec) || !isipow2(lenSpec)) throw OptiVec_Pow2SizeErr();
			if(X.u_ht < 2*htSpec || X.u_len < 2*lenSpec) throw OptiVec_BigSizeErr();
			if(Win.u_ht != 2*htSpec || Win.u_len != 2*lenSpec) throw OptiVec_SizeErr();
			MFlspectrum(vx_m0, htSpec, lenSpec, X.vx_m0, X.u_ht, X.u_len, Win.vx_m0);
		}

		void l_convolve(tMatObj* Flt, cmta X, cmta Rsp)
		{
			enforcePow2();
			if((u_ht!=Flt->u_ht)||(u_len!=Flt->u_len)||(u_ht!=X.u_ht)||(u_len!=X.u_len)||
				(u_ht!=Rsp.u_ht)||(u_len!=Rsp.u_len)) throw OptiVec_SizeErr();
			if((mx_matr==Rsp.mx_matr)||(Flt->mx_matr==X.mx_matr))
				throw OptiVec_Err();
			MFlconvolve(vx_m0, Flt->vx_m0, X.vx_m0, Rsp.vx_m0, u_ht, u_len);
		}

		void l_convolve(tMatObj Flt, cmta X, cmta Rsp)
		{
			enforcePow2();
			if((u_ht!=Flt.u_ht)||(u_len!=Flt.u_len)||(u_ht!=X.u_ht)||(u_len!=X.u_len)||
				(u_ht!=Rsp.u_ht)||(u_len!=Rsp.u_len)) throw OptiVec_SizeErr();
			if((mx_matr==Rsp.mx_matr)||(Flt.mx_matr==X.mx_matr))
				throw OptiVec_Err();
			MFlconvolve(vx_m0, Flt.vx_m0, X.vx_m0, Rsp.vx_m0, u_ht, u_len);
		}

		void l_deconvolve(tMatObj* Flt, cmta X, cmta Rsp)
		{
			enforcePow2();
			if((u_ht!=Flt->u_ht)||(u_len!=Flt->u_len)||(u_ht!=X.u_ht)||(u_len!=X.u_len)||
				(u_ht!=Rsp.u_ht)||(u_len!=Rsp.u_len)) throw OptiVec_SizeErr();
			if((mx_matr==Rsp.mx_matr)||(Flt->mx_matr==X.mx_matr))
				throw OptiVec_Err();
			MFldeconvolve(vx_m0, Flt->vx_m0, X.vx_m0, Rsp.vx_m0, u_ht, u_len);
		}

		void l_deconvolve(tMatObj Flt, cmta X, cmta Rsp)
		{
			enforcePow2();
			if((u_ht!=Flt.u_ht)||(u_len!=Flt.u_len)||(u_ht!=X.u_ht)||(u_len!=X.u_len)||
				(u_ht!=Rsp.u_ht)||(u_len!=Rsp.u_len)) throw OptiVec_SizeErr();
			if((mx_matr==Rsp.mx_matr)||(Flt.mx_matr==X.mx_matr))
				throw OptiVec_Err();
			MFldeconvolve(vx_m0, Flt.vx_m0, X.vx_m0, Rsp.vx_m0, u_ht, u_len);
		}

		void s_FFT(cmta X, int dir=1)
		{
			enforcePow2();
			MFsFFT(vx_m0,X.vx_m0,u_ht,u_len,dir);
		}

		//s_FFTtoC member of cfMatObj

		void s_filter(cmta X, cmta Flt)
		{
			enforcePow2();
			MFsfilter(vx_m0, X.vx_m0, Flt.vx_m0, u_ht, u_len);
		}

		void s_autocorr(cmta X)
		{
			enforcePow2();
			MFsautocorr(vx_m0, X.vx_m0, u_ht, u_len);
		}

		void s_xcorr(cmta X, cmta Y)
		{
			enforcePow2();
			MFsxcorr(vx_m0, X.vx_m0, Y.vx_m0, u_ht, u_len);
		}

		void s_spectrum(cmta X, cmta Win)
		{
			unsigned htSpec=u_ht-1, lenSpec=u_len-1;
			if(!isipow2(htSpec) || !isipow2(lenSpec)) throw OptiVec_Pow2SizeErr();
			if(X.u_ht < 2*htSpec || X.u_len < 2*lenSpec) throw OptiVec_BigSizeErr();
			if(Win.u_ht != 2*htSpec || Win.u_len != 2*lenSpec) throw OptiVec_SizeErr();
			MFsspectrum(vx_m0, htSpec, lenSpec, X.vx_m0, X.u_ht, X.u_len, Win.vx_m0);
		}

		void s_convolve(tMatObj* Flt, cmta X, cmta Rsp)
		{
			enforcePow2();
			if((u_ht!=Flt->u_ht)||(u_len!=Flt->u_len)||(u_ht!=X.u_ht)||(u_len!=X.u_len)||
				(u_ht!=Rsp.u_ht)||(u_len!=Rsp.u_len)) throw OptiVec_SizeErr();
			if((mx_matr==Rsp.mx_matr)||(Flt->mx_matr==X.mx_matr))
				throw OptiVec_Err();
			MFsconvolve(vx_m0, Flt->vx_m0, X.vx_m0, Rsp.vx_m0, u_ht, u_len);
		}

		void s_convolve(tMatObj Flt, cmta X, cmta Rsp)
		{
			enforcePow2();
			if((u_ht!=Flt.u_ht)||(u_len!=Flt.u_len)||(u_ht!=X.u_ht)||(u_len!=X.u_len)||
				(u_ht!=Rsp.u_ht)||(u_len!=Rsp.u_len)) throw OptiVec_SizeErr();
			if((mx_matr==Rsp.mx_matr)||(Flt.mx_matr==X.mx_matr))
				throw OptiVec_Err();
			MFsconvolve(vx_m0, Flt.vx_m0, X.vx_m0, Rsp.vx_m0, u_ht, u_len);
		}

		void s_deconvolve(tMatObj* Flt, cmta X, cmta Rsp)
		{
			enforcePow2();
			if((u_ht!=Flt->u_ht)||(u_len!=Flt->u_len)||(u_ht!=X.u_ht)||(u_len!=X.u_len)||
				(u_ht!=Rsp.u_ht)||(u_len!=Rsp.u_len)) throw OptiVec_SizeErr();
			if((mx_matr==Rsp.mx_matr)||(Flt->mx_matr==X.mx_matr))
				throw OptiVec_Err();
			MFsdeconvolve(vx_m0, Flt->vx_m0, X.vx_m0, Rsp.vx_m0, u_ht, u_len);
		}

		void s_deconvolve(tMatObj Flt, cmta X, cmta Rsp)
		{
			enforcePow2();
			if((u_ht!=Flt.u_ht)||(u_len!=Flt.u_len)||(u_ht!=X.u_ht)||(u_len!=X.u_len)||
				(u_ht!=Rsp.u_ht)||(u_len!=Rsp.u_len)) throw OptiVec_SizeErr();
			if((mx_matr==Rsp.mx_matr)||(Flt.mx_matr==X.mx_matr))
				throw OptiVec_Err();
			MFsdeconvolve(vx_m0, Flt.vx_m0, X.vx_m0, Rsp.vx_m0, u_ht, u_len);
		}

		void FFT(cmta X, int dir=1)
		{
			enforcePow2();
			MFFFT(vx_m0,X.vx_m0,u_ht,u_len,dir);
		}

		//FFTtoC member of cfMatObj

		void filter(cmta X, cmta Flt)
		{
			enforcePow2();
			MFfilter(vx_m0, X.vx_m0, Flt.vx_m0, u_ht, u_len);
		}

		void autocorr(cmta X)
		{
			enforcePow2();
			MFautocorr(vx_m0, X.vx_m0, u_ht, u_len);
		}

		void xcorr(cmta X, cmta Y)
		{
			enforcePow2();
			MFxcorr(vx_m0, X.vx_m0, Y.vx_m0, u_ht, u_len);
		}

		void spectrum(cmta X, cmta Win)
		{
			unsigned htSpec=u_ht-1, lenSpec=u_len-1;
			if(!isipow2(htSpec) || !isipow2(lenSpec)) throw OptiVec_Pow2SizeErr();
			if(X.u_ht < 2*htSpec || X.u_len < 2*lenSpec) throw OptiVec_BigSizeErr();
			if(Win.u_ht != 2*htSpec || Win.u_len != 2*lenSpec) throw OptiVec_SizeErr();
			MFspectrum(vx_m0, htSpec, lenSpec, X.vx_m0, X.u_ht, X.u_len, Win.vx_m0);
		}

		void convolve(tMatObj* Flt, cmta X, cmta Rsp)
		{
			enforcePow2();
			if((u_ht!=Flt->u_ht)||(u_len!=Flt->u_len)||(u_ht!=X.u_ht)||(u_len!=X.u_len)||
				(u_ht!=Rsp.u_ht)||(u_len!=Rsp.u_len)) throw OptiVec_SizeErr();
			if((mx_matr==Rsp.mx_matr)||(Flt->mx_matr==X.mx_matr))
				throw OptiVec_Err();
			MFconvolve(vx_m0, Flt->vx_m0, X.vx_m0,
				Rsp.vx_m0, u_ht, u_len);
		}

		void convolve(tMatObj Flt, cmta X, cmta Rsp)
		{
			enforcePow2();
			if((u_ht!=Flt.u_ht)||(u_len!=Flt.u_len)||(u_ht!=X.u_ht)||(u_len!=X.u_len)||
				(u_ht!=Rsp.u_ht)||(u_len!=Rsp.u_len)) throw OptiVec_SizeErr();
			if((mx_matr==Rsp.mx_matr)||(Flt.mx_matr==X.mx_matr))
				throw OptiVec_Err();
			MFconvolve(vx_m0, Flt.vx_m0, X.vx_m0, Rsp.vx_m0, u_ht, u_len);
		}

		void deconvolve(tMatObj* Flt, cmta X, cmta Rsp)
		{
			enforcePow2();
			if((u_ht!=Flt->u_ht)||(u_len!=Flt->u_len)||(u_ht!=X.u_ht)||(u_len!=X.u_len)||
				(u_ht!=Rsp.u_ht)||(u_len!=Rsp.u_len)) throw OptiVec_SizeErr();
			if((mx_matr==Rsp.mx_matr)||(Flt->mx_matr==X.mx_matr))
				throw OptiVec_Err();
			MFdeconvolve(vx_m0, Flt->vx_m0, X.vx_m0, Rsp.vx_m0, u_ht, u_len);
		}

		void deconvolve(tMatObj Flt, cmta X, cmta Rsp)
		{
			enforcePow2();
			if((u_ht!=Flt.u_ht)||(u_len!=Flt.u_len)||(u_ht!=X.u_ht)||(u_len!=X.u_len)||
				(u_ht!=Rsp.u_ht)||(u_len!=Rsp.u_len)) throw OptiVec_SizeErr();
			if((mx_matr==Rsp.mx_matr)||(Flt.mx_matr==X.mx_matr))
				throw OptiVec_Err();
			MFdeconvolve(vx_m0, Flt.vx_m0, X.vx_m0, Rsp.vx_m0, u_ht, u_len);
		}

/////////////////////////////////////////////////////////////////////////////////////////
//Data fitting
// All of the fitting functions themselves belong to fVecObj.
// Only a few helper functions follow here:

		void setLinfitNeglect( const t& Thresh ) const
		{ VF_setLinfitNeglect( Thresh ); }

		t getLinfitNeglect() const
		{  return VF_getLinfitNeglect(); }

		t nonlinfit_getChi2() const
		{ return MF_nonlinfit_getChi2();  }

		t nonlinfitwW_getChi2() const
		{ return MF_nonlinfitwW_getChi2();  }

		unsigned nonlinfit_getTestRun() const
		{   return MF_nonlinfit_getTestRun(); }

		unsigned nonlinfitwW_getTestRun() const
		{   return MF_nonlinfitwW_getTestRun(); }

		unsigned nonlinfit_getTestPar() const
		{ return MF_nonlinfit_getTestPar(); }

		unsigned nonlinfitwW_getTestPar() const
		{ return MF_nonlinfitwW_getTestPar(); }

		int  nonlinfit_getTestDir() const
		{ return MF_nonlinfit_getTestDir(); }

		int  nonlinfitwW_getTestDir() const
		{ return MF_nonlinfitwW_getTestDir(); }

		void  nonlinfit_stop() const
		{  MF_nonlinfit_stop(); }

		void  nonlinfitwW_stop() const
		{  MF_nonlinfitwW_stop(); }

/////////////////////////////////////////////////////////////////////////////////////////
//Input and Output
//
		void fprint(FILE* stream, unsigned linewidth) const
		{
			enforceSize();
			MFfprint(stream,vx_m0,u_ht, u_len,linewidth);
		}

#if !defined _Windows || defined __FLAT__ || defined _WIN32
		void cprint() const
		{
			enforceSize();
			MFcprint(vx_m0,u_ht,u_len);
		}
#endif
		void print() const
		{
			fprint(stdout,80);
		}

		void write(FILE* stream) const
		{
			enforceSize();
			MFwrite(stream,vx_m0,u_ht,u_len);
		}

		void read(FILE* stream)
		{
			enforceSize();
			MFread(vx_m0,u_ht,u_len,stream);
		}

		void setWriteFormat(char* formatString) const
		{
			VF_setWriteFormat(formatString);
		}

		void setWriteSeparate(char* sepString) const
		{
			VF_setWriteSeparate(sepString);
		}

		void store(FILE* stream) const
		{
			enforceSize();
			fwrite(&u_ht,sizeof(unsigned),1,stream);
			fwrite(&u_len,sizeof(unsigned),1,stream);
			VF_store(stream,vx_m0,ui_size);
		}

		void recall(FILE* stream)
		{
			fread(&u_ht,sizeof(unsigned),1,stream);
			fread(&u_len,sizeof(unsigned),1,stream);
			setSize(u_ht,u_len);
			VF_recall(vx_m0,ui_size,stream);
		}

		#ifdef __MGRAPH_H
			void xyzAutoDensityMap( cvta X, cvta Y, COLORREF mincol, COLORREF maxcol ) const
			{
				enforceSize((unsigned)X.ui_size,(unsigned)Y.ui_size);
				MFxyzAutoDensityMap( X.vx_vect, Y.vx_vect, vx_m0, u_ht, u_len,
														mincol, maxcol );
			}

			void xyzDataDensityMap( cvta X, cvta Y ) const
			{
				enforceSize((unsigned)X.ui_size,(unsigned)Y.ui_size);
				MFxyzDataDensityMap( X.vx_vect, Y.vx_vect, vx_m0, u_ht, u_len);
			}

			void zAutoDensityMap( COLORREF mincol, COLORREF maxcol ) const
			{
				enforceSize();
				MFzAutoDensityMap( vx_m0, u_ht, u_len, mincol, maxcol );
			}

			void zDataDensityMap() const
			{
				enforceSize();
				MFzDataDensityMap( vx_m0, u_ht, u_len);
			}

		#endif   // __MGRAPH_H
	};  // end of class matrix<T>


	inline vt tVecObj::operator*(cmta MY)
	{
		tVecObj tmp(this->ui_size);
		tmp.enforceSize();
		VFmulM(tmp.vx_vect,vx_vect,MY.vx_m0,MY.u_ht,MY.u_len);
		return tmp;
	}

	inline void tVecObj::Row_extract(cmta MA, const unsigned iRow)
	{
		if( MA.u_ht <= iRow || ui_size < MA.u_len || !ui_size ) throw OptiVec_SizeErr();
		MFRow_extract( vx_vect, MA.vx_m0, MA.u_ht, MA.u_len, iRow );
	}

	inline void tVecObj::Col_extract(cmta MA, const unsigned iCol)
	{
		if( MA.u_len <= iCol || ui_size < MA.u_ht || !ui_size ) throw OptiVec_SizeErr();
		MFCol_extract( vx_vect, MA.vx_m0, MA.u_ht, MA.u_len, iCol );
	}

	inline void tVecObj::Dia_extract(cmta MA)
	{
		if( !ui_size || ui_size != MA.u_len || (MA.u_len != MA.u_ht) )throw OptiVec_SizeErr();
		MFDia_extract( vx_vect, MA.vx_m0, MA.u_len );
	}

	inline void tVecObj::Rows_max(cmta MA)
	{
		if( !MA.u_ht || !MA.u_len || ui_size < MA.u_ht ) throw OptiVec_SizeErr();
		MFRows_max( vx_vect, MA.vx_m0, MA.u_ht, MA.u_len );
	}

	inline void tVecObj::Cols_max(cmta MA)
	{
		if( !MA.u_ht || !MA.u_len || ui_size < MA.u_len ) throw OptiVec_SizeErr();
		MFCols_max( vx_vect, MA.vx_m0, MA.u_ht, MA.u_len );
	}

	inline void tVecObj::Rows_min(cmta MA)
	{
		if( !MA.u_ht || !MA.u_len || ui_size < MA.u_ht ) throw OptiVec_SizeErr();
		MFRows_min( vx_vect, MA.vx_m0, MA.u_ht, MA.u_len );
	}

	inline void tVecObj::Cols_min(cmta MA)
	{
		if( !MA.u_ht || !MA.u_len || ui_size < MA.u_len ) throw OptiVec_SizeErr();
		MFCols_min( vx_vect, MA.vx_m0, MA.u_ht, MA.u_len );
	}

	inline void tVecObj::Rows_absmax(cmta MA)
	{
		if( !MA.u_ht || !MA.u_len || ui_size < MA.u_ht ) throw OptiVec_SizeErr();
		MFRows_absmax( vx_vect, MA.vx_m0, MA.u_ht, MA.u_len );
	}

	inline void tVecObj::Cols_absmax(cmta MA)
	{
		if( !MA.u_ht || !MA.u_len || ui_size < MA.u_len ) throw OptiVec_SizeErr();
		MFCols_absmax( vx_vect, MA.vx_m0, MA.u_ht, MA.u_len );
	}

	inline void tVecObj::Rows_absmin(cmta MA)
	{
		if( !MA.u_ht || !MA.u_len || ui_size < MA.u_ht ) throw OptiVec_SizeErr();
		MFRows_absmin( vx_vect, MA.vx_m0, MA.u_ht, MA.u_len );
	}

	inline void tVecObj::Cols_absmin(cmta MA)
	{
		if( !MA.u_ht || !MA.u_len || ui_size < MA.u_len ) throw OptiVec_SizeErr();
		MFCols_absmin( vx_vect, MA.vx_m0, MA.u_ht, MA.u_len );
	}

	inline void tVecObj::Rows_sum(cmta MA)
	{
		if( !MA.u_ht || !MA.u_len || ui_size < MA.u_ht ) throw OptiVec_SizeErr();
		MFRows_sum( vx_vect, MA.vx_m0, MA.u_ht, MA.u_len );
	}

	inline void tVecObj::Cols_sum(cmta MA)
	{
		if( !MA.u_ht || !MA.u_len || ui_size < MA.u_len ) throw OptiVec_SizeErr();
		MFCols_sum( vx_vect, MA.vx_m0, MA.u_ht, MA.u_len );
	}

	inline void tVecObj::Rows_prod(cmta MA)
	{
		if( !MA.u_ht || !MA.u_len || ui_size < MA.u_ht ) throw OptiVec_SizeErr();
		MFRows_prod( vx_vect, MA.vx_m0, MA.u_ht, MA.u_len );
	}

	inline void tVecObj::Cols_prod(cmta MA)
	{
		if( !MA.u_ht || !MA.u_len || ui_size < MA.u_len ) throw OptiVec_SizeErr();
		MFCols_prod( vx_vect, MA.vx_m0, MA.u_ht, MA.u_len );
	}

	inline void tVecObj::MmulV( cmta MA, cvta X )
	{
		if( ui_size != MA.u_ht || X.ui_size != MA.u_len ) throw OptiVec_SizeErr();
		MFmulV( vx_vect, MA.vx_m0, X.vx_vect, MA.u_ht, MA.u_len );
	}

	inline void tVecObj::MTmulV( cmta MA, cvta X )
	{
		if( ui_size != MA.u_len || X.ui_size != MA.u_ht ) throw OptiVec_SizeErr();
		MFTmulV( vx_vect, MA.vx_m0, X.vx_vect, MA.u_ht, MA.u_len );
	}

	inline void tVecObj::mulM( cvta X, cmta MA )
	{
		if( ui_size != MA.u_len || X.ui_size != MA.u_ht ) throw OptiVec_SizeErr();
		VFmulM( vx_vect, X.vx_vect, MA.vx_m0, MA.u_ht, MA.u_len );
	}

	inline void tVecObj::mulMT( cvta X, cmta MA )
	{
		if( ui_size != MA.u_ht || X.ui_size != MA.u_len ) throw OptiVec_SizeErr();
		VFmulMT( vx_vect, X.vx_vect, MA.vx_m0, MA.u_ht, MA.u_len );
	}

	inline void tVecObj::LUsolve( cmta MLU, cvta B, const uVecObj& Ind )
	{
		if( ui_size != MLU.u_len || ui_size != MLU.u_ht ||
		    ui_size != B.ui_size || ui_size != Ind.ui_size ) throw OptiVec_SizeErr();
		MFLUsolve( vx_vect, MLU.vx_m0, B.vx_vect, Ind.vx_vect, MLU.u_len );
	}

	inline void tVecObj::LUimprove( cvta B, cmta MA, cmta MLU, const uVecObj& Ind )
	{
		if( ui_size != MLU.u_len || ui_size != MLU.u_ht ||
		    ui_size != MA.u_len || ui_size != MA.u_ht ||
		    ui_size != B.ui_size || ui_size != Ind.ui_size ) throw OptiVec_SizeErr();
		MFLUimprove( vx_vect, B.vx_vect, MA.vx_m0, MLU.vx_m0, Ind.vx_vect, MLU.u_len );
	}

	inline void tVecObj::SVsolve( cmta MU, cmta MV, cvta W, cvta B )
	{ // result may overwrite B only for "normal", not for over- or under-determined systems
		if( MU.u_ht < MU.u_len || MU.u_len != MV.u_len || MV.u_len != MV.u_ht ||
		    B.ui_size != ui_size || W.ui_size != MV.u_len || ui_size != MU.u_ht )
		    throw OptiVec_SizeErr();
		MFSVsolve( vx_vect, MU.vx_m0, MV.vx_m0, W.vx_vect, B.vx_vect, (unsigned)B.ui_size, MU.u_len );
	}

	inline int tVecObj::solve( cmta MA, cvta B )
	{
		if( ui_size != MA.u_len || ui_size != MA.u_ht ||
		    ui_size != B.ui_size ) throw OptiVec_SizeErr();
		return MFsolve( vx_vect, MA.vx_m0, B.vx_vect, MA.u_len );
	}

	inline int tVecObj::solveBySVD( cmta MA, cvta B )
	{
		if(ui_size != MA.u_len || (B.ui_size != MA.u_ht && B.vx_vect != vx_vect)) throw OptiVec_SizeErr();
		return MFsolveBySVD( vx_vect, MA.vx_m0, B.vx_vect, MA.u_ht, MA.u_len );
	}

	inline int tVecObj::safeSolve( cmta MA, cvta B )
	{
		if( ui_size != MA.u_len || ui_size != MA.u_ht ||
		    ui_size != B.ui_size ) throw OptiVec_SizeErr();
		return MFsafeSolve( vx_vect, MA.vx_m0, B.vx_vect, MA.u_len );
	}

	inline int tVecObj::sym_eigenvalues( tMatObj EigM, cmta MA, int CalcEigenVec )
	{
		if( ui_size != MA.u_len || ui_size != MA.u_ht ||
		    (CalcEigenVec && (ui_size != EigM.u_len || ui_size != EigM.u_ht)) )
		    throw OptiVec_SizeErr();
		return( MFs_eigenvalues( vx_vect, EigM.vx_m0, MA.vx_m0, MA.u_len, CalcEigenVec));
	}

	inline int tVecObj::sym_eigenvalues( tMatObj* EigM, cmta MA, int CalcEigenVec )
	{
		if( ui_size != MA.u_len || ui_size != MA.u_ht ||
		    (CalcEigenVec && (ui_size != EigM->u_len || ui_size != EigM->u_ht)) )
		throw OptiVec_SizeErr();
		return( MFs_eigenvalues( vx_vect, EigM->vx_m0, MA.vx_m0, MA.u_len, CalcEigenVec));
	}

	inline void tVecObj::polyfit( cvta X, cvta Y )
	{	unsigned deg = ui_size-1;
		if( deg >= X.ui_size || X.ui_size != Y.ui_size) throw OptiVec_SizeErr();
		VF_polyfit( vx_vect, deg, X.vx_vect, Y.vx_vect, X.ui_size );
	}

	inline void tVecObj::polyfitwW( tMatObj Covar, cvta X, cvta Y, cvta InvVar )
	{	unsigned deg = ui_size-1;
		if( deg >= X.ui_size || X.ui_size != Y.ui_size || X.ui_size != InvVar.ui_size ||
		    Covar.u_ht != Covar.u_len || Covar.u_ht != (unsigned)ui_size) throw OptiVec_SizeErr();
		VFpolyfitwW( vx_vect, Covar.vx_m0, deg, X.vx_vect, Y.vx_vect, InvVar.vx_vect, X.ui_size );
	}

	inline void tVecObj::polyfitwW( tMatObj* Covar, cvta X, cvta Y, cvta InvVar )
	{	unsigned deg = ui_size-1;
		if( deg >= X.ui_size || X.ui_size != Y.ui_size || X.ui_size != InvVar.ui_size ||
		    Covar->u_ht != Covar->u_len || Covar->u_ht != (unsigned)ui_size) throw OptiVec_SizeErr();
		VFpolyfitwW( vx_vect, Covar->vx_m0, deg,
	 	            X.vx_vect, Y.vx_vect, InvVar.vx_vect, X.ui_size );
	}

	inline void tVecObj::linfit( const iVecObj& AStatus, cvta X, cvta Y,
									void __vf funcs(tVector BasFuncs, t x, unsigned nfuncs) )
	{
		if( ui_size > X.ui_size || ui_size != AStatus.ui_size ||
		    X.ui_size != Y.ui_size) throw OptiVec_SizeErr();
		VF_linfit( vx_vect, AStatus.vx_vect, ui_size, X.vx_vect, Y.vx_vect, X.ui_size, funcs );
	}

	inline void tVecObj::linfitwW( tMatObj Covar, const iVecObj& AStatus, cvta X, cvta Y, cvta InvVar,
	             void __vf funcs(tVector BasFuncs, t x, unsigned nfuncs) )
	{
		if( ui_size > X.ui_size || ui_size != AStatus.ui_size ||
		    X.ui_size != Y.ui_size || X.ui_size != InvVar.ui_size ||
		    Covar.u_ht != Covar.u_len || Covar.u_ht != (unsigned)ui_size) throw OptiVec_SizeErr();
		VFlinfitwW( vx_vect, Covar.vx_m0, AStatus.vx_vect, ui_size,
		            X.vx_vect, Y.vx_vect, InvVar.vx_vect, X.ui_size, funcs );
	}

	inline void tVecObj::linfitwW( tMatObj* Covar, const iVecObj& AStatus, cvta X, cvta Y, cvta InvVar,
	             void __vf funcs(tVector BasFuncs, t x, unsigned nfuncs) )
	{
		if( ui_size > X.ui_size || ui_size != AStatus.ui_size ||
		    X.ui_size != Y.ui_size || X.ui_size != InvVar.ui_size ||
		    Covar->u_ht != Covar->u_len || Covar->u_ht != (unsigned)ui_size) throw OptiVec_SizeErr();
		VFlinfitwW( vx_vect, Covar->vx_m0, AStatus.vx_vect, ui_size,
		            X.vx_vect, Y.vx_vect, InvVar.vx_vect, X.ui_size, funcs );
	}

	inline void tVecObj::linfit( const iVecObj& AStatus, cvta X, cvta Y, cmta MZ,
	            void __vf funcs(tVector BasFuncs, t x, t y, unsigned nfuncs) )
	{
		if( ui_size > MZ.ui_size || ui_size != AStatus.ui_size ||
		    X.ui_size != MZ.u_len || Y.ui_size != MZ.u_ht) throw OptiVec_SizeErr();
		MFlinfit( vx_vect, AStatus.vx_vect, ui_size,
		          X.vx_vect, Y.vx_vect, MZ.vx_m0, MZ.u_ht, MZ.u_len, funcs );
	}

	inline void tVecObj::linfitwW( tMatObj Covar, const iVecObj& AStatus,
		      cvta X, cvta Y, cmta MZ, cmta MInvVar,
		      void __vf funcs(tVector BasFuncs, t x, t y, unsigned nfuncs) )
	{
		if( ui_size > MZ.ui_size || ui_size != AStatus.ui_size ||
		    X.ui_size != MZ.u_len || Y.ui_size != MZ.u_ht ||
		    MZ.u_ht != MInvVar.u_ht || MZ.u_len != MInvVar.u_len ||
		    Covar.u_ht != Covar.u_len || Covar.u_ht != (unsigned)ui_size)
		    throw OptiVec_SizeErr();
		MFlinfitwW( vx_vect, Covar.vx_m0, AStatus.vx_vect, ui_size,
		    X.vx_vect, Y.vx_vect, MZ.vx_m0, MInvVar.vx_m0, MZ.u_ht, MZ.u_len, funcs );
	}

	inline void tVecObj::linfitwW( tMatObj* Covar, const iVecObj& AStatus,
		           cvta X, cvta Y, cmta MZ, cmta MInvVar,
		           void __vf funcs(tVector BasFuncs, t x, t y, unsigned nfuncs) )
	{
		if( ui_size > MZ.ui_size || ui_size != AStatus.ui_size ||
		    X.ui_size != MZ.u_len || Y.ui_size != MZ.u_ht ||
		    MZ.u_ht != MInvVar.u_ht || MZ.u_len != MInvVar.u_len ||
		    Covar->u_ht != Covar->u_len || Covar->u_ht != (unsigned)ui_size)
		    throw OptiVec_SizeErr();
		MFlinfitwW( vx_vect, Covar->vx_m0, AStatus.vx_vect, ui_size,
		            X.vx_vect, Y.vx_vect, MZ.vx_m0, MInvVar.vx_m0, MZ.u_ht, MZ.u_len, funcs );
	}

	inline t tVecObj::nonlinfit( const iVecObj& AStatus, cvta X, cvta Y,
		           void __vf modelfunc(tVector YModel, tVector XModel, ui size),
		           void __vf derivatives(tVector dYdAi,tVector X, ui size, unsigned iPar) )
	{
		if( ui_size > X.ui_size || ui_size != AStatus.ui_size ||
		    X.ui_size != Y.ui_size) throw OptiVec_SizeErr();
		return VF_nonlinfit( vx_vect, AStatus.vx_vect, ui_size, X.vx_vect, Y.vx_vect, X.ui_size,
		                     modelfunc, derivatives );
	}

	inline void tVecObj::setLinfitNeglect( const t& Thresh ) const
	{ VF_setLinfitNeglect( Thresh ); }

	inline t tVecObj::getLinfitNeglect() const
	{  return VF_getLinfitNeglect(); }

	inline t tVecObj::nonlinfitwW( tMatObj Covar, const iVecObj& AStatus, cvta X, cvta Y, cvta InvVar,
		        void __vf modelfunc(tVector YModel, tVector XModel, ui size),
		        void __vf derivatives(tVector dYdAi,tVector X, ui size, unsigned iPar) )
	{
		if( ui_size > X.ui_size || ui_size != AStatus.ui_size ||
		    X.ui_size != Y.ui_size || X.ui_size != InvVar.ui_size ||
		    Covar.u_ht != Covar.u_len || Covar.u_ht != (unsigned)ui_size) throw OptiVec_SizeErr();
		return VFnonlinfitwW( vx_vect, Covar.vx_m0, AStatus.vx_vect, ui_size,
		              X.vx_vect, Y.vx_vect, InvVar.vx_vect, X.ui_size,
		              modelfunc, derivatives );
	}

	inline t tVecObj::nonlinfitwW( tMatObj* Covar, const iVecObj& AStatus, cvta X, cvta Y, cvta InvVar,
		        void __vf modelfunc(tVector YModel, tVector XModel, ui size),
		        void __vf derivatives(tVector dYdAi,tVector X, ui size, unsigned iPar) )
	{
		if( ui_size > X.ui_size || ui_size != AStatus.ui_size ||
		    X.ui_size != Y.ui_size || X.ui_size != InvVar.ui_size ||
		    Covar->u_ht != Covar->u_len || Covar->u_ht != (unsigned)ui_size)
		    throw OptiVec_SizeErr();
		return VFnonlinfitwW( vx_vect, Covar->vx_m0, AStatus.vx_vect, ui_size,
		               X.vx_vect, Y.vx_vect, InvVar.vx_vect, X.ui_size,
		               modelfunc, derivatives );
	}

	inline t tVecObj::nonlinfit( const iVecObj& AStatus, cvta X, cvta Y, cmta MZ,
		      void __vf modelfunc(tMatrix MZModel, unsigned htZ, unsigned lenZ, tVector X, tVector Y ),
		      void __vf derivatives(tMatrix dZdAi, unsigned htZ, unsigned lenZ, tVector X, tVector Y, unsigned ipar) )
	{
		if( ui_size > MZ.ui_size || ui_size != AStatus.ui_size ||
		    X.ui_size != MZ.u_len || Y.ui_size != MZ.u_ht) throw OptiVec_SizeErr();
		return MFnonlinfit( vx_vect, AStatus.vx_vect, ui_size,
		           X.vx_vect, Y.vx_vect, MZ.vx_m0, MZ.u_ht, MZ.u_len,
		           modelfunc, derivatives );
	}

	inline t tVecObj::nonlinfitwW( tMatObj Covar, const iVecObj& AStatus,
		      cvta X, cvta Y, cmta MZ, cmta MInvVar,
		      void __vf modelfunc(tMatrix MZModel, unsigned htZ, unsigned lenZ, tVector X, tVector Y ),
		      void __vf derivatives(tMatrix dZdAi, unsigned htZ, unsigned lenZ, tVector X, tVector Y, unsigned ipar) )
	{
		if( ui_size > MZ.ui_size || ui_size != AStatus.ui_size ||
		    X.ui_size != MZ.u_len || Y.ui_size != MZ.u_ht ||
		    MZ.u_ht != MInvVar.u_ht || MZ.u_len != MInvVar.u_len ||
		    Covar.u_ht != Covar.u_len || Covar.u_ht != (unsigned)ui_size)
		    throw OptiVec_SizeErr();
		return MFnonlinfitwW( vx_vect, Covar.vx_m0, AStatus.vx_vect, ui_size,
		          X.vx_vect, Y.vx_vect, MZ.vx_m0, MInvVar.vx_m0,
		          MZ.u_ht, MZ.u_len, modelfunc, derivatives );
	}

	inline t tVecObj::nonlinfitwW( tMatObj* Covar, const iVecObj& AStatus,
		      cvta X, cvta Y, cmta MZ, cmta MInvVar,
		      void __vf modelfunc(tMatrix MZModel, unsigned htZ, unsigned lenZ, tVector X, tVector Y ),
		      void __vf derivatives(tMatrix dZdAi, unsigned htZ, unsigned lenZ, tVector X, tVector Y, unsigned ipar) )
	{
		if( ui_size > MZ.ui_size || ui_size != AStatus.ui_size ||
		    X.ui_size != MZ.u_len || Y.ui_size != MZ.u_ht ||
		    MZ.u_ht != MInvVar.u_ht || MZ.u_len != MInvVar.u_len ||
		    Covar->u_ht != Covar->u_len || Covar->u_ht != (unsigned)ui_size)
		    throw OptiVec_SizeErr();
		return MFnonlinfitwW( vx_vect, Covar->vx_m0, AStatus.vx_vect, ui_size,
		              X.vx_vect, Y.vx_vect, MZ.vx_m0, MInvVar.vx_m0,
		              MZ.u_ht, MZ.u_len, modelfunc, derivatives );
	}

// all multifit functions can be called only in their normal C form,
// because MF_EXPERIMENT cannot work with VecObjects, but requires actual pointers

#ifdef DMATOBJ_H
	inline void fMatObj::DtoF(const dMatObj& X)
	{
		enforceSize(X.u_ht, X.u_len);
		V_DtoF(vx_m0,X.vx_m0,ui_size);
	}

	inline void dMatObj::FtoD(const fMatObj& X)
	{
		enforceSize(X.u_ht, X.u_len);
		V_FtoD(vx_m0,X.vx_m0,ui_size);
	}
#endif //DMATOBJ_H

#ifdef EMATOBJ_H
#ifdef __BORLANDC__ //80-bit numbers defined
	inline void fMatObj::EtoF(const eMatObj& X)
	{
		enforceSize(X.u_ht, X.u_len);
		V_EtoF(vx_m0,X.vx_m0,ui_size);
	}

	inline void eMatObj::FtoE(const fMatObj& X)
	{
		enforceSize(X.u_ht, X.u_len);
		V_FtoE(vx_m0,X.vx_m0,ui_size);
	}
#endif //__BORLANDC__
#endif //EMATOBJ_H

#if defined _MSC_VER || (defined __BORLANDC__ && __BORLANDC__ >= 0x500)
}  // end of namespace OptiVec
#endif

#undef cmta
#undef cvta
#undef cua
#undef cta

#undef vta
#undef mta
#undef cuint
#undef ta

#undef vt
#undef mt
#undef t
#undef tVector
#undef tMatrix

#undef vtcplx
#undef tcplx

#undef tMatObj
#undef tVecObj
#ifdef __BORLANDC__
#pragma warn .inl
#endif

#define FMATOBJ_H
#endif //FMATOBJ_H
