/* ceMatObj.H

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

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

#ifndef CEMATOBJ_H

#include <VecObj.h>
#include <ceVecObj.h>
#include <eMatObj.h>
#include <uVecObj.h>
#include <iVecObj.h>
#ifdef __BORLANDC__   // 80-bit reals supported ?

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

#define vta  vector<eComplex>&
#define mta  matrix<eComplex>&
#define ta   eComplex&

#define vt   vector<eComplex>
#define mt   matrix<eComplex>
#define t    eComplex
#define tVector ceVector
#define tMatrix ceMatrix
#define tMatObj ceMatObj
#define tVecObj ceVecObj

#define tbase  extended

#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=MCE_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.Re && !fill.Im) //fill == 0.0
		 	 	{	mx_matr=MCE_matrix0(ht,len); vx_m0=mx_matr[0];}
				else
				{	mx_matr=MCE_matrix(ht,len);  vx_m0=mx_matr[0];
					VCE_equC(vx_m0,ui_size,fill);}
	 		}
		}

		matrix(cmta init)
		{	if(init.ui_size)
			{	mx_matr= MCE_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;
				VCE_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) free(); //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= MCE_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) free(); //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= MCE_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; }

		void 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)
		{
			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();
			VCE_equV(vx_m0, matr.vx_m0, ui_size);
			return *this;}

		mta operator=(cta scalar)
		{	enforceSize();
			VCE_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);
			VCE_accV(vx_m0, X.vx_m0, ui_size);
			return *this;}

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

		mta operator*=(cta scalar)
		{	enforceSize();
			VCE_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);
			VCE_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();
			VCE_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);
			VCE_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();
			VCE_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();
			MCEmulM(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();
			MCEmulV(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 ceVecObj, see below

		mt operator*(cta C)
		{	tMatObj tmp(this->u_ht, this->u_len);
			tmp.enforceSize();
			VCE_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();
			VCE_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);
			VCE_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);
			VCE_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);
			VCE_mulC(tmp.vx_m0,X.vx_m0,tmp.ui_size,C);
			return tmp;}

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

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

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


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

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


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

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

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

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

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

		void outerprod(cvta X, cvta Y)
		{
			enforceSize( X.ui_size, Y.ui_size);
			MCEouterprod(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();
			MCERow_equ0(vx_m0, u_ht, u_len, iRow);
		}

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

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

		void Row_equC(cuint iRow, cta C)
		{
			if( !u_len || (u_ht <= iRow) )throw OptiVec_SmallSizeErr();
			MCERow_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();
			MCECol_equC(vx_m0, u_ht, u_len, iCol, C);
		}

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

		void Row_equV(cuint iRow, cvta X)
		{
			if( !u_len || (u_ht <= iRow) )throw OptiVec_SmallSizeErr();
			MCERow_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();
			MCECol_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();
			MCEDia_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();
			MCETrd_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();
			MCETrd_extract(vx_m0, MA.vx_m0, u_len);
		}

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

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

		void conj(cmta X)
		{
			enforceSize(X.u_ht, X.u_len);
			MCEconj(vx_m0, X.vx_m0, u_ht, u_len);
		}

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

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

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

		//defined by including cdMatObj.h
		void CDtoCE(const cdMatObj& 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();
			MCEsubmatrix(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();
			MCEsubmatrix_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();
			MCEequMblock(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();
			MCEequMblockT(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();
			MCEblock_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();
			MCEblock_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
				MCERow_delete(vx_m0, vx_m0, u_ht--, u_len, iRow);
			}
			else
			{	enforceSize(MA.u_ht-1, MA.u_len);
				MCERow_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 = MCE_matrix( u_ht, u_len-1 );
				MCECol_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);
				MCECol_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 = MCE_matrix( u_ht, u_len );
				MCERow_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, MA.u_len+1);
				MCERow_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 = MCE_matrix( u_ht, u_len );
				MCECol_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);
				MCECol_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();
			MCERow_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();
			MCECol_addC(vx_m0, u_ht, u_len, iCol, C);
		}

	 	 void Dia_addC(cta C)
		{
			if( !u_len || (u_len != u_ht) )throw OptiVec_SmallSizeErr();
			MCEDia_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();
			MCERow_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();
			MCECol_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();
			MCEDia_addV(vx_m0, u_len, X.vx_vect);
		}

		void Row_subC(cuint iRow, cta C)
		{
			if( !u_len || (u_ht <= iRow) )throw OptiVec_SmallSizeErr();
			MCERow_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();
			MCECol_subC(vx_m0, u_ht, u_len, iCol, C);
		}

	 	 void Dia_subC(cta C)
		{
			if( !u_len || (u_len != u_ht) )throw OptiVec_SmallSizeErr();
			MCEDia_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();
			MCERow_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();
			MCECol_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();
			MCEDia_subV(vx_m0, u_len, X.vx_vect);
		}

		void Row_subrC(cuint iRow, cta C)
		{
			if( !u_len || (u_ht <= iRow) )throw OptiVec_SmallSizeErr();
			MCERow_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();
			MCECol_subrC(vx_m0, u_ht, u_len, iCol, C);
		}

	 	void Dia_subrC(cta C)
		{
			if( !u_len || (u_len != u_ht) )throw OptiVec_SmallSizeErr();
			MCEDia_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();
			MCERow_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();
			MCECol_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();
			MCEDia_subrV(vx_m0, u_len, X.vx_vect);
		}

		void Row_mulC(cuint iRow, cta C)
		{
			if( !u_len || (u_ht <= iRow) )throw OptiVec_SmallSizeErr();
			MCERow_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();
			MCECol_mulC(vx_m0, u_ht, u_len, iCol, C);
		}

	 	void Dia_mulC(cta C)
		{
			if( !u_len || (u_len != u_ht) )throw OptiVec_SmallSizeErr();
			MCEDia_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();
			MCERow_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();
			MCECol_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();
			MCEDia_mulV(vx_m0, u_len, X.vx_vect);
		}

		void Row_divC(cuint iRow, cta C)
		{
			if( !u_len || (u_ht <= iRow) )throw OptiVec_SmallSizeErr();
			MCERow_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();
			MCECol_divC(vx_m0, u_ht, u_len, iCol, C);
		}

	 	void Dia_divC(cta C)
		{
			if( !u_len || (u_len != u_ht) )throw OptiVec_SmallSizeErr();
			MCEDia_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();
			MCERow_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();
			MCECol_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();
			MCEDia_divV(vx_m0, u_len, X.vx_vect);
		}

		void Row_divrC(cuint iRow, cta C)
		{
			if( !u_len || (u_ht <= iRow) )throw OptiVec_SmallSizeErr();
			MCERow_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();
			MCECol_divrC(vx_m0, u_ht, u_len, iCol, C);
		}

	 	void Dia_divrC(cta C)
		{
			if( !u_len || (u_len != u_ht) )throw OptiVec_SmallSizeErr();
			MCEDia_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();
			MCERow_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();
			MCECol_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();
			MCEDia_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 eVecObj.h)

		t Dia_maxReIm()
		{
			if( !u_len || (u_len != u_ht) )throw OptiVec_SizeErr();
			return MCEDia_maxReIm(vx_m0, u_len);
		}

		t Dia_minReIm()
		{
			if( !u_len || (u_len != u_ht) )throw OptiVec_SizeErr();
			return MCEDia_minReIm(vx_m0, u_len);
		}

		tbase Dia_absmax()
		{
			if( !u_len || (u_len != u_ht) )throw OptiVec_SizeErr();
			return MCEDia_absmax(vx_m0, u_len);
		}

		tbase Dia_absmin()
		{
			if( !u_len || (u_len != u_ht) )throw OptiVec_SizeErr();
			return MCEDia_absmin(vx_m0, u_len);
		}

		t Dia_absmaxReIm()
		{
			if( !u_len || (u_len != u_ht) )throw OptiVec_SizeErr();
			return MCEDia_absmaxReIm(vx_m0, u_len);
		}

		t Dia_absminReIm()
		{
			if( !u_len || (u_len != u_ht) )throw OptiVec_SizeErr();
			return MCEDia_absminReIm(vx_m0, u_len);
		}

		t Dia_sabsmax()
		{
			if( !u_len || (u_len != u_ht) )throw OptiVec_SizeErr();
			return MCEDia_sabsmax(vx_m0, u_len);
		}

		t Dia_sabsmin()
		{
			if( !u_len || (u_len != u_ht) )throw OptiVec_SizeErr();
			return MCEDia_sabsmin(vx_m0, u_len);
		}

		t Dia_cabsmax()
		{
			if( !u_len || (u_len != u_ht) )throw OptiVec_SizeErr();
			return MCEDia_cabsmax(vx_m0, u_len);
		}

		t Dia_cabsmin()
		{
			if( !u_len || (u_len != u_ht) )throw OptiVec_SizeErr();
			return MCEDia_cabsmin(vx_m0, u_len);
		}

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

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

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

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

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

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

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

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

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

		void Cols_reflect()
		{
			enforceSize();
			MCECols_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();
			MCERows_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();
			MCECols_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();
			MCERows_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();
			MCECols_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();
			MCERows_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();
			MCECols_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();
			MCERows_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();
			MCECols_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();
			MCERows_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();
			MCECols_lincomb(vx_m0, u_ht, u_len, destCol, destC, srceCol, srceC);
		}

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

		void hermconj(cmta MA)
		{
			enforceSize(MA.u_len, MA.u_ht);
			MCEhermconj(vx_m0, MA.vx_m0, u_ht, 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);
			MCEaddM(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);
			MCEaddMT(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);
			MCEsubM(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);
			MCEsubMT(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);
			MCEsubM(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);
			MCEsubrMT(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);
			MCEmulC(vx_m0, MA.vx_m0, u_ht, u_len, C);
		}

		void divC(cmta MA, cta C)
		{
			enforceSize(MA.u_ht, MA.u_len);
			MCEdivC(vx_m0, MA.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);
			MCElincomb(vx_m0, MA.vx_m0, MB.vx_m0, u_ht, u_len, CA, CB);
		}

		//mixed multiplications matrix * vector belong to eVecObj
		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();
	 		MCEmulM( 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();
	 		MCEmulMT( vx_m0, MA.vx_m0, MB.vx_m0, MA.u_ht, MA.u_len, MB.u_ht);
		}

		void mulMH( cmta MA, cmta MB )
		{
	 		if( u_ht != MA.u_ht || u_len != MB.u_ht || MB.u_len != MA.u_len )throw OptiVec_SizeErr();
	 		MCEmulMH( 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();
	 		MCETmulM( vx_m0, MA.vx_m0, MB.vx_m0, MA.u_ht, MA.u_len, MB.u_len);
		}

		void HmulM( cmta MA, cmta MB )
		{
	 		if( u_ht != MA.u_len || u_len != MB.u_len || MB.u_ht != MA.u_ht )throw OptiVec_SizeErr();
	 		MCEHmulM( 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();
	 		MCETmulMT( 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();
	 		MCEdiamulM( 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();
	 		MCEdiamulMT( 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();
			MCEmulMdia( 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();
			MCETmulMdia( 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 MCELUdecompose( 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 MCELUdecompose( vx_m0, Ind->vx_vect, MA.vx_m0, u_len );
		}

		int  LUDresult()
		{
			return MCE_LUDresult();
		}

		void LUDsetEdit(const tbase& Thresh)
		{
			MCE_LUDsetEdit(Thresh);
		}

		tbase LUDgetEdit()
		{
			return MCE_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();
			MCELUinv( vx_m0, MLU.vx_m0, Ind.vx_vect, u_len );
		}

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

/*  SVD not yet implemented:
		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 MCESVdecompose( 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 MCESVdecompose( vx_m0, MV.vx_m0, W.vx_vect,
		                           MA.vx_m0, MA.u_ht, MA.u_len );
		}

		void SVDsetEdit(cta Thresh)
		{
			MCE_SVDsetEdit(Thresh);
		}

		t SVDgetEdit()
		{
			return MCE_SVDgetEdit();
		}
*/

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

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


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

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

		void p_FFTtoC(const matrix<tbase>& X)
		{
			enforcePow2();
			MEpFFTtoC(vx_m0,X.vx_m0,u_ht,u_len);
		}

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

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

		void l_FFTtoC(const matrix<tbase>& X)
		{
			enforcePow2();
			MElFFTtoC(vx_m0,X.vx_m0,u_ht,u_len);
		}

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

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

		void s_FFTtoC(const matrix<tbase>& X)
		{
			enforcePow2();
			MEsFFTtoC(vx_m0,X.vx_m0,u_ht,u_len);
		}

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

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

		void FFTtoC(const matrix<tbase>& X)
		{
			enforcePow2();
			MEFFTtoC(vx_m0,X.vx_m0,u_ht,u_len);
		}

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


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

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

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

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

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

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

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

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


	inline vt tVecObj::operator*(cmta MY)
	{
		tVecObj tmp(this->ui_size);
		tmp.enforceSize();
		VCEmulM(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();
		MCERow_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();
		MCECol_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();
		MCEDia_extract( vx_vect, MA.vx_m0, MA.u_len );
	}

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

	inline void vector<tbase>::Cols_absmin(cmta MA)
	{
		if( !MA.u_ht || !MA.u_len || ui_size < MA.u_len ) throw OptiVec_SizeErr();
		MCECols_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();
		MCERows_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();
		MCECols_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();
		MCERows_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();
		MCECols_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();
		MCEmulV( 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();
		MCETmulV( 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();
		VCEmulM( 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();
		VCEmulMT( 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();
		MCELUsolve( 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();
		MCELUimprove( vx_vect, B.vx_vect, MA.vx_m0, MLU.vx_m0, Ind.vx_vect, MLU.u_len );
	}

/* SVD not yet implemented:
	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();
		MCESVsolve( 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 MCEsolve( 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 MCEsolveBySVD( 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 MCEsafeSolve( vx_vect, MA.vx_m0, B.vx_vect, MA.u_len );
	}
*/

#ifdef CFMATOBJ_H
	inline void cfMatObj::CEtoCF(const ceMatObj& X)
	{
		enforceSize(X.u_ht, X.u_len);
		V_CEtoCF(vx_m0,X.vx_m0,ui_size);
	}

	inline void ceMatObj::CFtoCE(const cfMatObj& X)
	{
		enforceSize(X.u_ht, X.u_len);
		V_CFtoCE(vx_m0,X.vx_m0,ui_size);
	}
#endif //CEMATOBJ_H

#ifdef CDMATOBJ_H
	inline void cdMatObj::CEtoCD(const ceMatObj& X)
	{
		enforceSize(X.u_ht, X.u_len);
		V_CEtoCD(vx_m0,X.vx_m0,ui_size);
	}

	inline void ceMatObj::CDtoCE(const cdMatObj& X)
	{
		enforceSize(X.u_ht, X.u_len);
		V_CDtoCE(vx_m0,X.vx_m0,ui_size);
	}
#endif //CDMATOBJ_H

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

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

#undef vta
#undef mta
#undef ta

#undef vt
#undef mt
#undef t
#undef tVector
#undef tMatrix
#undef tVecObj
#undef tMatObj

#undef vtbase
#undef tbase

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

#define CEMATOBJ_H
#endif //CEMATOBJ_H
