spamatrix.hh

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

const int SpaMatrix<Type>::nzRow(int i) const {  
  return (*irow)[i];}

				//row operations
template<class Type>
int SpaMatrix<Type>::rowCount(int ir) const
{
    
   int k = (*irow)[ir+1]-(*irow)[ir]; 
    
   if (ir==nrow-1) k++;   
   return k;
}

 
template<class Type>
Type SpaMatrix<Type>::rowDot(int ir, const Vector<Type>& x) 
{
    
   assert(ncol==x.size()); 
    
   int n = (*irow)[ir]; 
    
   Type sum = 0;
    
   for (int i=n; ( i<(*irow)[ir+1]);  i++)
      sum += (*elem)[i]*x[(*icol)[i]];
    
   return sum; 
}

 
template<class Type>
Vector<Type> SpaMatrix<Type>::rowScale(int ir, Type c) 
{
    
   Vector<Type> v(ncol); 
    
   int n = (*irow)[ir];

    
   for (int i=0;   i<ncol;   i++)
   {
      if ((*icol)[n] == i) 
      {
	 v[i] = c*(*elem)[n]; n++; 
      }
      else v[i] = (Type)0; 
   }
    
   return v;
}

 
template<class Type>
Type SpaMatrix<Type>::rowNorm(int ir, int p) 	//row-norm of the sparse matrix
{
    
   int nelem = rowCount(ir); 
    
   Vector<Type> v(nelem);
    
   for (int i=0;   i<nelem;   i++) v[i] = (*elem)[irow[0][ir]+i]; 
    
   return v.norm(p); 
}
 
template<class Type>
Type SpaMatrix<Type>::rowNorm2(int ir) 	//row-2-norm squared of the sparse matrix
{
    
  Type p = rowNorm(ir,2);  
  return p*p;
}

				// Overloading operators
template<class Type>
SpaMatrix<Type>& SpaMatrix<Type>::operator=(const SpaMatrix<Type>& A){
 
  assert(nz==A.nz);  
  assert(nrow==A.nrow);  
  assert(ncol==A.ncol);

  *elem = *A.elem; 
  *icol = *A.icol; 
  *irow = *A.irow; 
 
  return *this;
}
 
template<class Type>
SpaMatrix<Type>&  SpaMatrix<Type>::operator+=(Type c){

  *elem += c;  
  return *this;
}
 
template<class Type>
SpaMatrix<Type>&  SpaMatrix<Type>::operator-=(Type c){

  *elem -= c;  
  return *this;
}
 
template<class Type>
SpaMatrix<Type>&  SpaMatrix<Type>::operator*=(Type c){

  *elem *= c;  
  return *this;
} 
 
template<class Type>
SpaMatrix<Type>&  SpaMatrix<Type>::operator/=(Type c){

  *elem /= c;  
  return *this;
}


				//friends functions
template<class Type>
Vector<Type> operator*(const SpaMatrix<Type>& A, const Vector<Type>& v) 
{
    
   assert(A.numOfCols()==v.size());
    
   int n = A.numOfRows();
    
   Vector<Type> u(n);

    
   int i, j;
    
   for(i=0;   i<n;   i++)
     for(j = A.nzRow(i); (j<A.nzRow(i+1));  j++)
       u[i] +=  A.nzElem(j)*v[A.colPos(j)];
    
   return u;
}

 
template<class Type>
Vector<Type> SpaMatrix<Type>::adotx(const Vector<Type>& x)
{
    
   SpaMatrix<Type> A(*this);
    
   return A*x;
}

 
template<class Type>
Vector<Type> SpaMatrix<Type>::atdotx(const Vector<Type>& x)
{
    
   if (x.size() != nrow) inValidSize();
    
   Vector<Type> v(ncol);

    
   int i, j = 0;
    
   for(i = 0;  i < nz;   i++)
   {
     if((i>=(*irow)[j+1])&&(i!=nz)) j++;
     v[(*icol)[i]] += (*elem)[i] * x[j];
   }
    
   return v;
}

 
template<class Type>			//sum across columns
Type SpaMatrix<Type>::operator() (int ir, char*)
{
     
    Type sum = 0;
     
    for (int j=(*irow)[ir]; ( j<(*irow)[ir+1]);   j++)
	sum += (*elem)[j];
     
    return sum;
}

 
template<class Type>			//sum across rows
Type SpaMatrix<Type>::operator() (char*, int ic)
{
     
    Type sum = 0;
     
    int k, iend = (*irow)[0];
	
     
    for (int j=0;   j<nrow;   j++){
	k = iend;
	iend = (*irow)[j+1];
	for(; (k<iend); k++)
	  if ((*icol)[k] == ic) { 
	    sum += (*elem)[k]; break;
	  }
    }
     
    return sum;
}

 
template<class Type>
Type SpaMatrix<Type>::rowMax(int ir) const //return the largest non-zero member of the row
{
     
    int j = (*irow)[ir]+1;
     
    Type m = (*elem)[j-1];
     
    for (; (j<(*irow)[ir+1]);  j++)
	m = Max(m, (*elem)[j]);
     
    return m;
}

 
template<class Type>
Type SpaMatrix<Type>::rowMin(int ir) const	//return the smallest non-zero member of the row
{
     
    int j = (*irow)[ir]+1;
     
    Type m = (*elem)[j-1];
     
    for (; (j<(*irow)[ir+1]);   j++)
	m = Min(m, (*elem)[j]);
     
    return m;
}

 
template<class Type>
Type SpaMatrix<Type>::colMax(int ic) const	//return the largest non-zero member of the col
{
     
    int k = 0;
     
    while (ic!=(*icol)[k]) k++;
     
    Type m = (*elem)[k];
     
    int ir = 0;
     
    while ((*irow)[ir] < k) ir++;
     
    int iend = (*irow)[ir+1];
     
    for (;   ir<nrow;   ir++) {
    	iend = (*irow)[ir+1];
	while (ic!=(*icol)[k] && k<iend) k++;
	if (k != iend) {
	    m = Max(m, (*elem)[k]);
	    k = iend;
	}
    }
     
    return m;
}
	
 
template<class Type>
Type SpaMatrix<Type>::colMin(int ic) const	//return the smallest non-zero member of the col
{
     
    int k = 0;
     
    while ( ic!=(*icol)[k]) k++;
     
    Type m = (*elem)[k];
     
    int ir = 0;
     
    while ((*irow)[ir] < k) ir++;
     
    int iend;
     
    for (;   ir<nrow;   ir++) {
	iend = (*irow)[ir+1];
	while (ic!=(*icol)[k] && k<iend) k++;
	if (k != iend ) {
	    m = Min(m, (*elem)[k]);
	    k = iend;
	    iend = (*irow)[ir+2];
	}
    }
     
    return m;
}


 
template<class Type>
SpaMatrix<Type> SpaMatrix<Type>::operator-(){
     
    SpaMatrix<Type> A(nz,nrow,ncol); 
    *(A.icol) = *icol; 
    *(A.irow) = *irow; 
    *(A.elem) = -(*elem);
     
    return A;
}
					//Overloading operators

template<class Type>
ostream& operator<<(ostream& osp, const SpaMatrix<Type>& A)
{
 
	int j, k = 0;
 
	for (int i=0;   i<A.nrow;   i++){
	  for(j=0; j<A.ncol; j++){
	    if((i<A.nrow-1 && k==(A.irow[0])[i+1]) || (A.icol[0])[k]!=j) osp<<"0  ";
	    else {osp<<(*A.elem)[k]<<"  "; k++; }
	  } osp<<endl;
	}
 
	return osp;
}

 
template <class Type>
inline SpaMatrix<Type>  operator+(const SpaMatrix<Type>& A, Type c){
     
  SpaMatrix<Type> B(A);  
  B += c;  
  return B;
}
 
template <class Type>
inline SpaMatrix<Type>  operator+(Type c, const SpaMatrix<Type>& A)
{ 
   SpaMatrix<Type> B(A);  
   B += c;  
   return B;
}
 
template <class Type>
inline SpaMatrix<Type>  operator-(const SpaMatrix<Type>& A, Type c){
     
  SpaMatrix<Type> B(A);  
  B -= c;  
  return B;
}
 
template <class Type>
inline SpaMatrix<Type>  operator-(Type c, const SpaMatrix<Type>& A){
     
  SpaMatrix<Type> B(A);  
  B -= c;  
  return -B;
}
 
template <class Type>
inline SpaMatrix<Type>  operator*(const SpaMatrix<Type>& A, Type c)
{
   SpaMatrix<Type> B(A);  
   B *= c;  
   return B;
}
 
template <class Type>
inline SpaMatrix<Type>  operator*(Type c, const SpaMatrix<Type>& A)
{ 
   SpaMatrix<Type> B(A);  
   B *= c;  
   return B;
}
 
template <class Type>
inline SpaMatrix<Type>  operator/(const SpaMatrix<Type>& A, Type c)
{
   SpaMatrix<Type> B(A);  
   B /= c;  
   return B;
}

#endif