/* cdMatObj.H

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

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

#ifndef CDMATOBJ_H

#include <VecObj.h>
#include <cdVecObj.h>
#include <dMatObj.h>
#include <uVecObj.h>
#include <iVecObj.h>

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

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

#define vt   vector<dComplex>
#define mt   matrix<dComplex>
#define t    dComplex
#define tVector cdVector
#define tMatrix cdMatrix
#define tMatObj cdMatObj
#define tVecObj cdVecObj

#define tbase  double

#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=MCD_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=MCD_matrix0(ht,len); vx_m0=mx_matr[0];}
				else
				{	mx_matr=MCD_matrix(ht,len);  vx_m0=mx_matr[0];
					VCD_equC(vx_m0,ui_size,fill);}
            }
		}

		matrix(cmta init)
		{	if(init.ui_size)
			{	mx_matr= MCD_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;
				VCD_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= MCD_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= MCD_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) 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();
			VCD_equV(vx_m0, matr.vx_m0, ui_size);
			return *this;}

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

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

		mta operator*=(cta scalar)
		{	enforceSize();
			VCD_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);
			VCD_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();
			VCD_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);
			VCD_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();
			VCD_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();
			MCDmulM(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();
			MCDmulV(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 cdVecObj, see below

		mt operator*(cta C)
		{	tMatObj tmp(this->u_ht, this->u_len);
			tmp.enforceSize();
			VCD_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();
			VCD_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);
			VCD_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);
			VCD_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);
			VCD_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= MCD_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= MCD_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 MCDPelement(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 MCDelement(vx_m0,u_ht,u_len,m,n);
		}


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

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

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

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

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

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

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

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

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

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

		void Row_equV(cuint iRow, cvta X)
		{
			if( !u_len || (u_ht <= iRow) )throw OptiVec_SmallSizeErr();
			MCDRow_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();
			MCDCol_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();
			MCDDia_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();
			MCDTrd_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();
			MCDTrd_extract(vx_m0, MA.vx_m0, u_len);
		}

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

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

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

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

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

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

		//defined by including ceMatObj.h
		void CEtoCD(const ceMatObj& 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();
			MCDsubmatrix(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();
			MCDsubmatrix_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();
			MCDequMblock(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();
			MCDequMblockT(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();
			MCDblock_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();
			MCDblock_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
				MCDRow_delete(vx_m0, vx_m0, u_ht--, u_len, iRow);
			}
			else
			{	enforceSize(MA.u_ht-1, MA.u_len);
				MCDRow_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 = MCD_matrix( u_ht, u_len-1 );
				MCDCol_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);
				MCDCol_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 = MCD_matrix( u_ht, u_len );
				MCDRow_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);
				MCDRow_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 = MCD_matrix( u_ht, u_len );
				MCDCol_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);
				MCDCol_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();
			MCDRow_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();
			MCDCol_addC(vx_m0, u_ht, u_len, iCol, C);
		}

		void Dia_addC(cta C)
		{
			if( !u_len || (u_len != u_ht) )throw OptiVec_SmallSizeErr();
			MCDDia_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();
			MCDRow_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();
			MCDCol_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();
			MCDDia_addV(vx_m0, u_len, X.vx_vect);
		}

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

		void Dia_subC(cta C)
		{
			if( !u_len || (u_len != u_ht) )throw OptiVec_SmallSizeErr();
			MCDDia_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();
			MCDRow_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();
			MCDCol_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();
			MCDDia_subV(vx_m0, u_len, X.vx_vect);
		}

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

		void Dia_subrC(cta C)
		{
			if( !u_len || (u_len != u_ht) )throw OptiVec_SmallSizeErr();
			MCDDia_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();
			MCDRow_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();
			MCDCol_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();
			MCDDia_subrV(vx_m0, u_len, X.vx_vect);
		}

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

		void Dia_mulC(cta C)
		{
			if( !u_len || (u_len != u_ht) )throw OptiVec_SmallSizeErr();
			MCDDia_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();
			MCDRow_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();
			MCDCol_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();
			MCDDia_mulV(vx_m0, u_len, X.vx_vect);
		}

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

		void Dia_divC(cta C)
		{
			if( !u_len || (u_len != u_ht) )throw OptiVec_SmallSizeErr();
			MCDDia_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();
			MCDRow_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();
			MCDCol_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();
			MCDDia_divV(vx_m0, u_len, X.vx_vect);
		}

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

		void Dia_divrC(cta C)
		{
			if( !u_len || (u_len != u_ht) )throw OptiVec_SmallSizeErr();
			MCDDia_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();
			MCDRow_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();
			MCDCol_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();
			MCDDia_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 dVecObj.h)

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

		void hermconj(cmta MA)
		{
			enforceSize(MA.u_len, MA.u_ht);
			MCDhermconj(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);
			MCDaddM(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);
			MCDaddMT(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);
			MCDsubM(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);
			MCDsubMT(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);
			MCDsubM(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);
			MCDsubrMT(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);
			MCDmulC(vx_m0, MA.vx_m0, u_ht, u_len, C);
		}

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

		//mixed multiplications matrix * vector belong to dVecObj
        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();
            MCDmulM( 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();
            MCDmulMT( 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();
            MCDmulMH( 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();
            MCDTmulM( 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();
            MCDHmulM( 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();
            MCDTmulMT( 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();
            MCDdiamulM( 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();
            MCDdiamulMT( 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();
            MCDmulMdia( 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();
            MCDTmulMdia( 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 MCDLUdecompose( 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 MCDLUdecompose( vx_m0, Ind->vx_vect, MA.vx_m0, u_len );
        }

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

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

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

        t LUdet( int permut ) const
        {
            if( u_ht != u_len ) throw OptiVec_SizeErr();
            return MCDLUdet( 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 MCDSVdecompose( 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 MCDSVdecompose( vx_m0, MV.vx_m0, W.vx_vect,
                                   MA.vx_m0, MA.u_ht, MA.u_len );
        }

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

		t SVDgetEdit() const
		{
			return MCD_SVDgetEdit();
		}
*/

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

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


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

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

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

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

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

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

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

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

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

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

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

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

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


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

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

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

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

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

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

		void store(FILE* stream) const
		{
			enforceSize();
			fwrite(&u_ht,sizeof(unsigned),1,stream);
			fwrite(&u_len,sizeof(unsigned),1,stream);
			VCD_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);
			VCD_recall(vx_m0,ui_size,stream);
		}
	};


	inline vt tVecObj::operator*(cmta MY)
	{
		tVecObj tmp(this->ui_size);
		tmp.enforceSize();
		VCDmulM(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();
		MCDRow_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();
		MCDCol_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();
		MCDDia_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();
		MCDRows_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();
        MCDCols_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();
        MCDRows_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();
        MCDCols_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();
        MCDRows_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();
        MCDCols_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();
        MCDRows_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();
        MCDCols_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();
        MCDRows_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();
        MCDCols_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();
        MCDRows_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();
        MCDCols_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();
        MCDRows_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();
        MCDCols_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();
        MCDRows_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();
        MCDCols_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();
        MCDRows_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();
        MCDCols_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();
        MCDRows_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();
        MCDCols_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();
        MCDRows_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();
        MCDCols_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();
        MCDRows_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();
        MCDCols_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();
        MCDmulV( 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();
        MCDTmulV( 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();
        VCDmulM( 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();
        VCDmulMT( 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();
        MCDLUsolve( 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();
        MCDLUimprove( 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();
        MCDSVsolve( 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 MCDsolve( 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 MCDsolveBySVD( 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 MCDsafeSolve( vx_vect, MA.vx_m0, B.vx_vect, MA.u_len );
    }
*/

#ifdef CFMATOBJ_H
	inline void cfMatObj::CDtoCF(const cdMatObj& X)
	{
		enforceSize(X.u_ht, X.u_len);
		V_CDtoCF(vx_m0,X.vx_m0,ui_size);
	}

	inline void cdMatObj::CFtoCD(const cfMatObj& X)
	{
		enforceSize(X.u_ht, X.u_len);
		V_CFtoCD(vx_m0,X.vx_m0,ui_size);
	}
#endif //CDMATOBJ_H

#ifdef CEMATOBJ_H
#ifdef __BORLANDC__ //80-bit numbers defined
	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 //__BORLANDC__
#endif //CEMATOBJ_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 CDMATOBJ_H
#endif //CDMATOBJ_H
