densmatrix.hh

COOOL：CWP面向对象最优化库（CWP Object Oriented Optimization Library） COOOL是C++类的一个集合

#ifndef DENSMATRIX_HH
#define DENSMATRIX_HH

#ifdef __GNUC__
#pragma interface
#endif

#include "Matrix.hh"

static const char*  myNamed =  "dense matrix class";

//===========================================
//@Man:
//@Memo: a dense matrix class template
//==========================================
/*@Doc: This is a class of regular dense matrix stored row-wise.

*/

template<class Type>
class DensMatrix : public Matrix<Type> {

  private:
    Vector<Type>  *aa;

  public:
    //@ManMemo: the default constructor
    DensMatrix();
    //@ManMemo: a constructor with specified number of rows and columns
    DensMatrix(///number of rows
	       int m, 
	       ///number of columns
	       int n);
    //@ManMemo: construct a densmatrix which has the contant of a vector array 
    DensMatrix(/// number of rows 
	       int m, 
	       /// pointer to a vector array
	       const Vector<Type>* p);
    //@ManMemo:  construct a densmatrix which non-zero elements only at the diagonal
    DensMatrix(const Vector<Type>& p);
    //@ManMemo:  copy a const dense matrix to a new object
    DensMatrix(const DensMatrix<Type>& A);
    ~DensMatrix(){ delete [] aa;}

    //@ManMemo: type name of the dense matrix
    const char* matrixType() const {return myNamed;}

    //@ManMemo: adjoint of the DensMatrix
    DensMatrix<Type> adjoint();

    //@ManMemo: access the ith row of the matrix, writable
    Vector<Type>& operator[](int i);
    //@ManMemo: access the ith row of the matrix, unwritable
    Vector<Type> operator[](int i) const;
    //@ManMemo: access the ith row of the matrix, unwritable
    Vector<Type> rowVector(int i) const {return aa[i];}
    //@ManMemo: access the jth column of the matrix, unwritable
    Vector<Type> colVector(int j) const
    {
       Vector<Type> u(nrow); 
       for(int i=0; i<nrow; i++) 
	  u[i] = aa[i][j]; 
       return u;  
    }

    ///augmenting the DensMatrix
    DensMatrix<Type>& chaSize(int, int);

    //@ManMemo:  sum across rows
    Type operator() (/// anything, could be empty
		     char*, 
		     /// the row index
		     int j)
    {
       Type sum=0;
       assert (j < ncol);
       for (int i = 0; i < nrow ; i++)
	  sum += aa[i][j];
       return sum;
    }

    //@ManMemo:  sum across columns
    Type operator() (/// the column index
		     int i, 
		     /// anything, could be empty
		     char*)
    {
       Type sum=0;
       assert (i < nrow);
       for (int j = 0; j < ncol; j++)  
	  sum += aa[i][j];
       return sum;
    }
    //@ManMemo:   return the largest element of the ith row
    Type rowMax(int i) const{ return aa[i].max(); }

    //@ManMemo:    return the smallest element of the ith row
    Type rowMin(int i) const{ return aa[i].min(); }
    //@ManMemo:    return the largest element of the jth col
    Type colMax(int j) const
    { 
       Type m = aa[0][j];	
       for(int i=1; i<nrow; i++) 
	  m = Max(m,aa[i][j]); 
       return m; 
    }	

    //@ManMemo:    return the smallest element of the jth col
    Type colMin(int j) const
    { 
       Type m = aa[0][j];	
       for(int i=1; i<nrow; i++) 
	  m = Min(m,aa[i][j]); 
       return m; 
    }	

            // implicit type conversion
    operator DensMatrix<int>() const;
    operator DensMatrix<long>() const;
    operator DensMatrix<float>() const;
    operator DensMatrix<double>() const;

    //@ManMemo:    Negating Operator
    DensMatrix<Type> operator-();
    //@ManMemo:    Assigning values of dense Matrix
    DensMatrix<Type>& operator= (const DensMatrix<Type>&);
    //@ManMemo:    Assigning values of a Vector to the diagonal
    DensMatrix<Type>& operator= (const Vector<Type>&);
    //@ManMemo:    Assigning a constant value
    DensMatrix<Type>& operator= (const Type);
    //@ManMemo:    Add a dense matrix
    DensMatrix<Type>& operator+= (const DensMatrix<Type>&);
    //@ManMemo:    Increment values by a constant
    DensMatrix<Type>& operator+= (const Type);
    //@ManMemo:    Subtract a dense matrix
    DensMatrix<Type>& operator-= (const DensMatrix<Type>&);
    //@ManMemo:    Decrement by a constant
    DensMatrix<Type>& operator-= (const Type);
    //@ManMemo:    Multiply by a constant
    DensMatrix<Type>& operator*= (const Type);
    //@ManMemo:    Divide by a constant
    DensMatrix<Type>& operator/= (const Type);

    //@ManMemo: insert a vector at the ith row
    DensMatrix<Type>& addVal(const int i, const Vector<Type>&);
    //@ManMemo: insert a dense matrix starting from the ith row
    DensMatrix<Type>& addVal(const int i, const DensMatrix<Type>&);
    //@ManMemo: insert a vector at the jth column
    DensMatrix<Type>& addVal(const Vector<Type>&, const int j);
    //@ManMemo: insert a dense matrix starting from the jth column
    DensMatrix<Type>& addVal(const DensMatrix<Type>&, const int j);
    //@ManMemo: delete ith row
    DensMatrix<Type>& delVal(const int i);
    
    //@ManMemo: saxpy, the current matrix $\times$ vector $x$
    Vector<Type> adotx(const Vector<Type>& x)
    {
       DensMatrix<Type>  a(*this);
       return a*x;
    }
    //@ManMemo: saxpy, the transpose of current matrix $\times$ vector $x$
    Vector<Type> atdotx(const Vector<Type>& x)
    {
       DensMatrix<Type>  at(this->adjoint());
       return at*x;
    }
    
    //@ManMemo: a friend, write matrix $A$ to stream $os$
friend ostream& operator<<(ostream& os, const DensMatrix<Type>& A);
    //@ManMemo: Matrix and Vector products, $A \times x$
friend Vector<Type> operator*(const DensMatrix<Type>& A, const Vector<Type>& x)
{
   int n = A.numOfRows();
   if (x.size() != n) inValidSize();
   Vector<Type> w(n);
   for (int i = 0; i < n; i++)  w[i] = A.rowVector(i)*x;
   return w;
}
    //@ManMemo: Matrix and Vector products, $A^T \times x$
friend Vector<Type> operator*(const Vector<Type>& x, const DensMatrix<Type>& A)
{
   int n = A.numOfRows();
   if (x.size() != n) inValidSize();
   Vector<Type> w(n);
   for (int i = 0; i < n; i++)  w[i] = A.rowVector(i)*x;
   return w;
}


};

#ifdef __GNUC__
#pragma implementation
#endif

  // Constructors for DensMatrix
template<class Type>
DensMatrix<Type>::DensMatrix()
{
   nrow = 1; 
   ncol = 1;
   aa = new Vector<Type>[1]; 
   aa[0].chaSize(1);
}

template<class Type>
DensMatrix<Type>::DensMatrix(int m, int n) 
{
   nrow = m; 
   ncol = n;
   aa = new Vector<Type>[m];
   for(int i=0; i<m; i++)
      aa[i].chaSize(n);
}

template<class Type>
DensMatrix<Type>::DensMatrix(int nn, const Vector<Type>* vp)
{ 
   nrow = nn; 
   ncol = vp->size();	
   aa = new Vector<Type>[nn]; 
   for (int i=0; i<nrow; i++) 
   { 
      aa[i].chaSize(ncol); 
      aa[i] = vp[i];
   } 
}

template<class Type>
DensMatrix<Type>::DensMatrix(const Vector<Type>& p)
{
   nrow = p.size();
   ncol = nrow;
   
   aa = new Vector<Type>[nrow];
   for(int i=0; i<nrow; i++)
   {
      aa[i].chaSize(ncol);
      aa[i][i] = p[i];
   }
}

template<class Type>
DensMatrix<Type>::DensMatrix(const DensMatrix<Type>& A)
{
   nrow = A.nrow; 
   ncol = A.ncol; 
   aa = new Vector<Type>[nrow];
   for (int i=0; i<nrow; i++) 
   {
      aa[i].chaSize(ncol); 
      aa[i] = (A.aa)[i];
   }
}

	//adjoint of the DensMatrix
template<class Type>
DensMatrix<Type> DensMatrix<Type>::adjoint()
{
   DensMatrix<Type> at(ncol,nrow);
   for (int i=0; i<ncol; i++)  
      at[i] = colVector(i);
   return at;
}

	//Fetch dimensions and overloading operators
template<class Type>
Vector<Type>& DensMatrix<Type>::operator[](int i) 
{return aa[i];}
template<class Type>
Vector<Type> DensMatrix<Type>::operator[](int i) const 
{return aa[i];}

			//augmenting the DensMatrix
template<class Type>
DensMatrix<Type>& DensMatrix<Type>::chaSize(int m, int n)
{
    DensMatrix<Type> A(*this);
    delete [] aa; 
    aa = new Vector<Type>[m];
    for(int i=0; i<m; i++) 
    {
       aa[i].chaSize(n); 
       if(i < nrow) 
       {
	  (A.aa)->chaSize(n); 
	  aa[i] = (A.aa)[i]; 
       }
    }
    nrow = m; 
    ncol = n;

    return *this;
}


template<class Type>   //overloading the operator -
DensMatrix<Type> DensMatrix<Type>::operator-()
{ 
   DensMatrix<Type> A(nrow,ncol);
   for (int i=0; i<nrow; i++) (A.aa)[i] = -aa[i];
   return A;
}


template<class Type>
DensMatrix<Type>& DensMatrix<Type>::operator=(const DensMatrix<Type>& B)
{
    if (nrow != B.nrow)
    {
       nrow = B.nrow;
       if (aa != NULL) delete aa;
       
       aa = new Vector<Type>[nrow];
       if (ncol != B.ncol)
       {
	  ncol = B.ncol;
	  for (int i=0; i<nrow; i++)
	     aa[i].chaSize(ncol);
       }
    }
    
    for (int i = 0; i < nrow; i++)  aa[i] = (B.aa)[i];
    return *this;
}

template<class Type>
DensMatrix<Type>& DensMatrix<Type>::operator= (const Vector<Type>& a)
{
    for (int i = 0; i < nrow; i++) aa[i][i] = a[i];
    return *this;
}
template<class Type>
DensMatrix<Type>& DensMatrix<Type>::operator= (const Type c)
{
    for (int i = 0; i < nrow; i++) aa[i] = c;
    return *this;
}
template<class Type>
DensMatrix<Type>& DensMatrix<Type>::operator+= (const DensMatrix<Type>& B)
{
    if (nrow != B.nrow) inValidSize();
    for (int i = 0; i < nrow; i++) aa[i] += B.aa[i];
    return *this;
}
template<class Type>
DensMatrix<Type>& DensMatrix<Type>::operator+= (const Type c)
{
    for (int i = 0; i < nrow; i++) aa[i] += c;
    return *this;
}
template<class Type>
DensMatrix<Type>& DensMatrix<Type>::operator-= (const DensMatrix<Type>& B)
{
    if (nrow != B.nrow) inValidSize();
    for (int i = 0; i < nrow; i++) aa[i] -= B.aa[i];
    return *this;
}
template<class Type>
DensMatrix<Type>& DensMatrix<Type>::operator-= (const Type c)
{
    for (int i = 0; i < nrow; i++) aa[i] -= c;
    return *this;
}
template<class Type>
DensMatrix<Type>& DensMatrix<Type>::operator*= (const Type c)
{
    for (int i = 0; i < nrow; i++) aa[i] *= c;
    return *this;
}
template<class Type>
DensMatrix<Type>& DensMatrix<Type>::operator/= (const Type c)
{
    assert (c!= 0);
    for (int i = 0; i < nrow; i++) aa[i] /= c;
    return *this;
}


//augment DensMatrix 
template<class Type>
DensMatrix<Type>& DensMatrix<Type>::addVal(const int, const Vector<Type>& v)
{
    if (v.size() != ncol) inValidSize();
    chaSize(nrow+1,ncol);
    aa[nrow-1] = v;
    return *this;
}

template<class Type>
DensMatrix<Type>& DensMatrix<Type>::delVal(const int irow)
{
   if (irow >= nrow) inValidSize();
   Vector<Type>* temp;
   temp = new Vector<Type>[nrow];
   
   int i;
   for (i=0; i<nrow; i++)
      temp[i] = aa[i];
   
   nrow --;
   delete [] aa;
   aa = new Vector<Type>[nrow];

   for (i=0; i<irow; i++)
      aa[i] = temp[i];
   for (i=irow; i<nrow; i++)
      aa[i] = temp[i+1];
   
   return *this;
}

template<class Type>
DensMatrix<Type>& DensMatrix<Type>::addVal(const int, const DensMatrix<Type>& A)
{
    if (ncol!=A.ncol) inValidSize();
    int nrowOld = nrow;
    chaSize(nrowOld+A.nrow,ncol);
    for (int i=nrowOld; i<nrow; i++) aa[i] = A.aa[i-nrowOld];
    return *this;
}

template<class Type>
DensMatrix<Type>& DensMatrix<Type>::addVal(const Vector<Type>& v, const int)
{	
    if (nrow!=v.size()) inValidSize();
    int ncolOld = ncol;
    chaSize(nrow,ncolOld+1);
    for (int i=0; i<nrow; i++)  aa[i][ncolOld] = v[i];
    return *this;
}

template<class Type>
DensMatrix<Type>& DensMatrix<Type>::addVal(const DensMatrix<Type>& A, const int)
{
    if (nrow!=A.nrow) inValidSize();
    for (int i=0; i<A.ncol; i++) addVal(A.colVector(i),1);
    return *this;
}

               // implicit type conversion
template<class Type>
DensMatrix<Type>::operator DensMatrix<int>() const
{
   DensMatrix<int> A(nrow,ncol);
   for (int i=0; i<nrow; i++) A[i] = aa[i];
   return A;
}
template<class Type>
DensMatrix<Type>::operator DensMatrix<long>() const
{
   DensMatrix<long> A(nrow,ncol);
   for (int i=0; i<nrow; i++) A[i] = aa[i]; 
   return A;
}
template<class Type>
DensMatrix<Type>::operator DensMatrix<float>() const
{
   DensMatrix<float> A(nrow,ncol);
   for (int i=0; i<nrow; i++) A[i] = aa[i];
   return A;
}
template<class Type>
DensMatrix<Type>::operator DensMatrix<double>() const
{
   DensMatrix<double> A(nrow,ncol);
   for (int i=0; i<nrow; i++) A[i] = aa[i]; 
   return A;
}

		// I/O for the DensMatrix  class
template<class Type>
ostream& operator<<(ostream& ofp, const DensMatrix<Type>& A)
{
   for (int i = 0; i < A.nrow; i++)
   {
      for (int j = 0; j < A.ncol; j++)
	 ofp << (A.aa)[i][j] << " ";
      ofp << endl;
   }
   return ofp;
}

//inline functions
template<class Type>
inline DensMatrix<Type> operator+(const DensMatrix<Type>& A, const DensMatrix<Type>& B)
{ 	
   DensMatrix<Type> C(A); 
   C += B; 
   return C; 
}
template<class Type>
inline 	DensMatrix<Type> operator+(const Type b, const DensMatrix<Type>& A)
{
   DensMatrix<Type> C(A);
   C += b;
   return C;
}
template<class Type>
inline 	DensMatrix<Type> operator+(const DensMatrix<Type>& A, const Type b)
{
   DensMatrix<Type> C(A);    
   C += b;    
   return C;
}
template<class Type>
inline 	DensMatrix<Type> operator-(const DensMatrix<Type>& A, const DensMatrix<Type>& B)
{
    DensMatrix<Type> C(A);    
    C -= B;    
    return C;
}
template<class Type>
inline 	DensMatrix<Type> operator-(const Type b, const DensMatrix<Type>& A)
{
    DensMatrix<Type> C(A);    
    C -= b;    
    return (-C);
}
template<class Type>
inline 	DensMatrix<Type> operator-(const DensMatrix<Type>& A, const Type b)
{
    DensMatrix<Type> C(A);    
    C -= b;    
    return C;
}

template<class Type>
inline DensMatrix<Type> operator* (const DensMatrix<Type>& A, const DensMatrix<Type>& B)
{
	if (A.numOfCols() != B.numOfRows()) inValidSize();
	int m = A.numOfRows(); n = B.numOfCols();
	DensMatrix<Type> C(m, n);
	int j;
	for (int i = 1; i < m; i++) 
	    for (j = 1; j < n; j++) 
		C[i][j] += A.rowVector(i)*A.colVector(j);
	return C;
}
template<class Type>
inline	DensMatrix<Type> operator* (const DensMatrix<Type>& A, const Type p)
{
	DensMatrix<Type> B(A);
	B *= p;
	return B;
}
template<class Type>
inline	DensMatrix<Type> operator* (const Type p, const DensMatrix<Type>& A)
{
	DensMatrix<Type> B(A);
	B *= p;
	return B;
}

template<class Type>
inline	DensMatrix<Type> operator/ (const DensMatrix<Type>& A, const Type p)
{
	DensMatrix<Type> B(A);
	B /= p;
	return B;
}

template<class Type>
inline DensMatrix<Type> outProduct(const Vector<Type> u, const Vector<Type> v)
{
   int n = u.size();
   assert(n == v.size());
   
   DensMatrix<Type> A(n,n);

   for (int i=0; i<n; i++)
      for (int j=0; j<n; j++)
	 A[i][j] = u[i]*v[j];

   return A;
}

#endif