spamatrix.hh

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

//============================================================
// COOOL           version 1.1           ---     Nov,  1995
//   Center for Wave Phenomena, Colorado School of Mines
//============================================================
//
//   This code is part of a preliminary release of COOOL (CWP
// Object-Oriented Optimization Library) and associated class 
// libraries. 
//
// The COOOL library is a free software. You can do anything you want
// with it including make a fortune.  However, neither the authors,
// the Center for Wave Phenomena, nor anyone else you can think of
// makes any guarantees about anything in this package or any aspect
// of its functionality.
//
// Since you've got the source code you can also modify the
// library to suit your own purposes. We would appreciate it 
// if the headers that identify the authors are kept in the 
// source code.



#ifndef SPAMATRIX_HH
#define SPAMATRIX_HH

//===========================
// Sparse-Matrix class template library for efficient algebraic operations
//
//  H. Lydia Deng,  01/24/1994
//===========================
// .NAME SpaMatrix class template, an inheritance of Matrix class
// .LIBRARY Base
// .HEADER c++ uitility classes
// .INCLUDE Matrix.hh
// .FILE SpaMatrix.hh

// .SECTION Description
// .B SpaMatrix<Type>
// is a derived sparse matrix class from Matrix
// it manages to do algebraic operations of two-dimensional sparse-matrix
// by considering the sparsity of the matrix
// it also include some basic row operations.
// 
// .SECTION Caveats

#ifdef __GNUC__
#pragma interface
#endif

#include "Vector.hh"
#include "Matrix.hh"

static const char*  myName =  "sparse matrix class";

//*****************************************************************************
//definition of the sparse SpaMatrix using row index class and its member functions
//@Man:
//@Memo: a sparse matrix class template
/*@Doc:  This is a class for efficient representation and manipulation of 
  sparse matrix. 
*/

template <class Type>
class SpaMatrix : public Matrix<Type> {
  private:
    int nz;
    //@ManMemo: 
    Vector<Type>* elem;
    //@ManMemo: 
    Vector<int>* icol;
    //@ManMemo: 
    Vector<int>* irow;
    public:			//constructors and destructor 
    //@ManMemo: 
    SpaMatrix(int,int,int); 
    //@ManMemo: 
    SpaMatrix(Vector<Type>&, Vector<int>&, Vector<int>&);
    //@ManMemo: 
    SpaMatrix(Vector<Type>&, int*, int*);
    //@ManMemo: 
    SpaMatrix(const SpaMatrix<Type>&);
    ~SpaMatrix();

    //@ManMemo: access the size description of the matrix
    const int nzSize() const;

    //@ManMemo: 
    const char* matrixType() const;

    //@ManMemo:the ith non-zero element
    Type& operator[](int );
    //@ManMemo:the ith row-vector 
    Vector<Type> rowVector(int i) const;	
    //@ManMemo:the jth coumn-vector 
    Vector<Type> colVector(int j) const;		

    //@ManMemo:sum across rows
    Type operator()(char*, int);
    //@ManMemo:sum across columns 
    Type operator()(int, char*);
    //@ManMemo:return the largest non-zero member of the row
    Type rowMax(int) const;	
    //@ManMemo:return the smallest non-zero member of the row
    Type rowMin(int) const;	
    //@ManMemo:return the largest non-zero member of the col
    Type colMax(int) const;	
    //@ManMemo:return the smallest non-zero member of the col
    Type colMin(int) const;	

    const Type nzElem(int i) const;
    //@ManMemo: 
    const int colPos(int i) const;
    //@ManMemo: 
    const int rowPos(int i) const;
    //@ManMemo: 
    const int nzRow(int i) const;

    //@ManMemo: count number of nonzeros 
    int rowCount(int) const;
    //@ManMemo: row-norm of the sparse matrix
    Type rowNorm(int ir, int p); 	
    //@ManMemo: row-2-norm squared of the sparse matrix
    Type rowNorm2(int ir);
    //@ManMemo: dot product of one row with a Vector
    Type rowDot(int, const Vector<Type>&);
    //@ManMemo: row scaling
    Vector<Type> rowScale(int, Type); 	

    //implicit type converting
    operator SpaMatrix<int>() const
    {
       int j, m = irow->size();
       Vector<int> mm(elem[0]);
       Vector<int> row(elem->size());
       for (int i=0; i<m; i++)
	  for (j=(*irow)[i]; j<(*irow)[i+1]); j++)  row[j] = i;
       SpaMatrix<int> A(mm,row,icol[0]);
       return A;
    }
    operator SpaMatrix<long>() const
    {
       int j, m = irow->size();
       Vector<long> mm(elem[0]);
       Vector<int> row(elem->size());
       for (int i=0; i<m; i++)
	  for (j=(*irow)[i]; (j<(*irow)[i+1]); j++)  row[j] = i;
       SpaMatrix<long> A(mm,row,icol[0]);
       return A;
    }
    operator SpaMatrix<float>() const
    {
       int j, m = irow->size();
       Vector<float> mm(elem[0]);
       Vector<int> row(elem->size());
       for (int i=0; i<m; i++)
	  for (j=(*irow)[i]; (j<(*irow)[i+1]); j++)  row[j] = i;
       SpaMatrix<float> A(mm,row,icol[0]);
       return A;
    }
    operator SpaMatrix<double>() const
    {
       int j, m = irow->size();
       Vector<double> mm(elem[0]);
       Vector<int> row(elem->size());
       for (int i=0; i<m; i++)
	  for (j=(*irow)[i]; (j<(*irow)[i+1]); j++)  row[j] = i;
       SpaMatrix<double> A(mm,row,icol[0]);
       return A;
    }
    //@ManMemo: 
    SpaMatrix<Type>&  	operator=(const SpaMatrix<Type>&);
    //@ManMemo: 
    SpaMatrix<Type>&  	operator+=(Type c);
    //@ManMemo: 
    SpaMatrix<Type>&  	operator-=(Type c);
    //@ManMemo: 
    SpaMatrix<Type>&  	operator*=(Type c);
    //@ManMemo: 
    SpaMatrix<Type>&  	operator/=(Type c);
    //@ManMemo: 
    SpaMatrix<Type> 	operator-();

    //@ManMemo: 
    Vector<Type> adotx(const Vector<Type>& x);
    //@ManMemo: 
    Vector<Type> atdotx(const Vector<Type>& x);

    //@ManMemo: 
friend Vector<Type> operator*(const SpaMatrix<Type>& A, const Vector<Type>& v);
    //@ManMemo: write the sparse matrix to a stream
friend ostream& operator<<(ostream&, const SpaMatrix<Type>&);
};

#ifdef __GNUC__
#pragma implementation
#endif

 
template<class Type>
SpaMatrix<Type>::SpaMatrix(int mm, int nn, int mn)
:Matrix<Type>(nn, mn) 
{
    
   nz = mm;
    
   elem = new Vector<Type>(nz);
    
   icol = new Vector<int>(nz);
    
   irow = new Vector<int>(nrow+1);
}

 
template<class Type>
SpaMatrix<Type>::SpaMatrix(Vector<Type>& x, Vector<int>& mm, Vector<int>& nn)
:Matrix<Type>(mm.max()+1, nn.max()+1) {
       
      nz = x.size();
       
      elem = new Vector<Type>(nz);
       
      icol = new Vector<int>(nz);
      
       
      for (int i=0;   i<nz;   i++)
      {
	 (*elem)[i] = x[i];
	 (*icol)[i] = nn[i];
      }

       
      irow = new Vector<int>(nrow+1);
       
      i = 0; (* irow)[i]=0;
       
      for (int j=1;   j<nz;   j++)
	 if (mm[j] != i) { assert(mm[j]==i+1); i++; (*irow)[i]=j; }
       
      assert(i==nrow-1);
      (* irow)[nrow] = nz;
}

 
template<class Type>
SpaMatrix<Type>::SpaMatrix(Vector<Type>& x, int* mm, int* nn)
{
    
   nz = x.size();
    
   Vector<int> vm(nz), vn(nz);
    
   vm = mm;
    
   vn = nn;

    
   SpaMatrix<Type> A(x,vm,vn);
    
   ncol = A.ncol;  
 nrow = A.nrow;
    
   elem = new Vector<Type>(nz);
    
   icol = new Vector<int>(nz);
    
   irow = new Vector<int>(nrow+1);

   * elem = *(A.elem);
   * icol = *(A.icol);
   * irow = *(A.irow);
}

 
template<class Type>
SpaMatrix<Type>::SpaMatrix(const SpaMatrix<Type>& A)
:Matrix<Type>(A.nrow, A.ncol) 
{
    
   nz = A.nz;
    
   elem = new Vector<Type>(nz);
    
   icol = new Vector<int>(nz);
    
   irow = new Vector<int>(nrow+1);
   * elem = *(A.elem);
   * icol = *(A.icol);
   * irow = *(A.irow);
}

 
template<class Type>
SpaMatrix<Type>::~SpaMatrix(){  
 delete elem;  
 delete icol;  
 delete irow;}

    				// access the size description of the matrix
template<class Type>   		// return  number of nonzero elements
const int SpaMatrix<Type>::nzSize() const { 
return nz;}

 
template<class Type>
const char* SpaMatrix<Type>::matrixType() const 
{ return myName;}

				//access elements of the matrix
template<class Type>            //the ith non-zero element
Type& SpaMatrix<Type>::operator[](int i) { 
return elem[0][i];}

 
template<class Type>
Vector<Type> SpaMatrix<Type>::colVector(int ic)	const	//the ith col-vector
{
     
    Vector<Type> v(nrow);
     
    int j, iend;
     
    for(int k=0;   k<nrow;   k++) {
        v[k] = 0;
        j = (*irow)[k];
        iend = (*irow)[k+1];
        if (k==nrow-1) iend ++;
        while ((*icol)[j] < ic && j<iend) j++;
        if ((*icol)[j]==ic && (j<iend)) v[k] = (*elem)[j];
    }
     
    return v;
}
 
 
template<class Type>
Vector<Type> SpaMatrix<Type>::rowVector(int ir) const	//the ith row-vector
{
     
    int j = (*irow)[ir];
     
    int iend = (*irow)[ir+1]-1;
     
    if (ir == nrow-1)  iend = (*irow)[ir];
     
    Vector<Type> v(ncol);
     
    for(int k=0;  k<(*icol)[j];  k++)
	v[k] = 0;
     
    for(; ( k<=(*icol)[iend]);   k++) {
	if ((*icol)[j] == k) {
	    v[k] = (*elem)[j];
	    j++;
	}
	else v[k] = 0;
    }
     
    for(;   k<ncol;   k++)
	v[k] = 0;
    
     
    return v;
}

 
template<class Type>
const Type SpaMatrix<Type>::nzElem(int i) const { 
return (*elem)[i];}
 
template<class Type>
const int SpaMatrix<Type>::colPos(int i) const { 
return (*icol)[i];}
 
template<class Type>
const int SpaMatrix<Type>::rowPos(int i) const 
{
    
   int j=1; 
    
   while((*irow)[j] <= i) j++;
   if (j>nrow) j--;  
   return (--j);
}
 
template<class Type>