matpackII.cc

Go to the documentation of this file.

/* Copyright (C) 2001-2008
   Stefan Buehler <sbuehler@ltu.se>
   Mattias Ekstroem <ekstrom@rss.chalmers.se>
 
   This program is free software; you can redistribute it and/or modify it
   under the terms of the GNU General Public License as published by the
   Free Software Foundation; either version 2, or (at your option) any
   later version.
 
   This program is distributed in the hope that it will be useful,
   but WITHOUT ANY WARRANTY; without even the implied warranty of
   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
   GNU General Public License for more details.
 
   You should have received a copy of the GNU General Public License
   along with this program; if not, write to the Free Software
   Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307,
   USA. */
 
// #include <vector>
#include <algorithm>
#include <set>
#include <iostream>             // For debugging.
#include <cmath>
#include <iterator>
#include "matpackII.h"
 
 
// Simple member Functions
// ----------------
 
Index Sparse::nrows() const
{
  return mrr;
}
 
Index Sparse::ncols() const
{
  return mcr;
}
 
Index Sparse::nnz() const
{
  return *(mcolptr->end()-1);
}
 
// Index Operators
// ---------------
 
 
Numeric& Sparse::rw(Index r, Index c)
{
  // Check if indices are valid:
  assert( 0<=r );
  assert( 0<=c );
  assert( r<mrr );
  assert( c<mcr );
 
  // Get index of first data element of this column:
  Index i = (*mcolptr)[c];
 
  // Get index of first data element of next column. (This always works,
  // because mcolptr has one extra element pointing behind the last
  // column.)
  const Index end = (*mcolptr)[c+1];
 
  // See if we find an element with the right row index in this
  // range. We assume that the elements are sorted by ascending row
  // index. We use the index i as the loop counter, which we have
  // initialized above.
  for ( ; i<end; ++i )
    {
      Index rowi = (*mrowind)[i];
      if ( r <  rowi )
        break;
      if ( r == rowi )
        return (*mdata)[i];
    }
 
  // If we are here, then the requested data element does not yet
  // exist. 
 
  // We have to adjust the array of column pointers. The values
  // in all columns above the current one have to be increased by 1.
  for ( vector<Index>::iterator j = mcolptr->begin() + c + 1;
        j < mcolptr->end();
        ++j )
    ++(*j);
 
  // We have to insert the new element in *mrowind and *mdata just one
  // position before the index i. We can use insert to achieve
  // this. Because they return an iterator to the newly inserted
  // element, we can return directly from the second call.
  mrowind->insert( mrowind->begin()+i, r );
  return *( mdata->insert( mdata->begin()+i, 0 ) );
}
 
 
Numeric Sparse::operator() (Index r, Index c) const
{ return this->ro(r, c); }
 
 
Numeric Sparse::ro (Index r, Index c) const
{
  // Check if indices are valid:
  assert( 0<=r );
  assert( 0<=c );
  assert( r<mrr );
  assert( c<mcr );
 
  // Get index of first data element of this column:
  Index begin = (*mcolptr)[c];
 
  // Get index of first data element of next column. (This always works,
  // because mcolptr has one extra element pointing behind the last
  // column.)
  const Index end = (*mcolptr)[c+1];
 
  // See if we find an element with the right row index in this
  // range. We assume that the elements are sorted by ascending row
  // index. 
  for ( Index i=begin; i<end; ++i )
    {
      Index rowi = (*mrowind)[i];
      if ( r <  rowi )
        return 0;
      if ( r == rowi )
        return (*mdata)[i];
    }
  return 0;
}
 
 
// Constructors
// ------------
 
 
Sparse::Sparse() :
  mdata(NULL),
  mrowind(NULL),
  mcolptr(NULL),
  mrr(0),
  mcr(0)
{
  // Nothing to do here
}
 
 
Sparse::Sparse(Index r, Index c) :
  mdata(new vector<Numeric>),
  mrowind(new vector<Index>),
  mcolptr(new vector<Index>(c+1,0)),
  mrr(r),
  mcr(c)
{
  // Nothing to do here.
}
 
 
 
Sparse::Sparse(const Sparse& m) :
//  mdata(new vector<Numeric>(*m.mdata)),
//  mrowind(new vector<Index>(*m.mrowind)),
//  mcolptr(new vector<Index>(*m.mcolptr)),    
  mdata(NULL),
  mrowind(NULL),
  mcolptr(NULL),    
  mrr(m.mrr),
  mcr(m.mcr)
{
  if (m.mdata)   mdata = new vector<Numeric>(*m.mdata);
  if (m.mrowind) mrowind = new vector<Index>(*m.mrowind);
  if (m.mcolptr) mcolptr = new vector<Index>(*m.mcolptr);
 
//   cout << "Original:\n" << m << "\n";
//   cout << "Copied:\n" << *this << "\n";
}
 
 
 
Sparse::~Sparse()
{
  delete mdata;
  delete mrowind;
  delete mcolptr;
}
 
 
void Sparse::insert_row(Index r, Vector v)
{
  // Check if the row index and the Vector length are valid
  assert( 0<=r );
  assert( r<mrr );
  assert( v.nelem()==mcr );
 
  // Calculate number of non-zero elements in Vector v.
  Index vnnz=0;
  for (Index i=0; i<v.nelem(); i++)
    {
      if (v[i]!=0)
        {
          vnnz++;
        }
    }
 
  // Count number of already existing elements in this row. Create
  // reference mrowind and mdata vector that copies of the real mrowind and
  // mdata. Resize the real mrowind and mdata to the correct output size.
  Index rnnz = vnnz - count(mrowind->begin(),mrowind->end(),r);
  vector<Index> mrowind_ref(mrowind->size());
  copy(mrowind->begin(), mrowind->end(), mrowind_ref.begin());
  vector<Numeric> mdata_ref(mdata->size());
  copy(mdata->begin(), mdata->end(), mdata_ref.begin());
  mrowind->resize(mrowind_ref.size()+rnnz);
  mdata->resize(mdata_ref.size()+rnnz);
 
  // Create iterators to the output vectors to keep track of current
  // positions.
  vector<Index>::iterator mrowind_it = mrowind->begin();
  vector<Numeric>::iterator mdata_it = mdata->begin();
 
  // Create a variable to store the change to mcolptr for each run
  Index colptr_mod = 0;
 
  // Loop through Vector v and insert the non-zero elements.
  for (Index i=0; i<v.nelem(); i++)
    {
      // Get mdata- and mrowind iterators to start and end of this
      // (the i:th) reference column.
      vector<Numeric>::iterator dstart =
        mdata_ref.begin()+(*mcolptr)[i];
      vector<Numeric>::iterator dend =
        mdata_ref.begin()+(*mcolptr)[i+1];
      vector<Index>::iterator rstart =
        mrowind_ref.begin()+(*mcolptr)[i];
      vector<Index>::iterator rend =
        mrowind_ref.begin()+(*mcolptr)[i+1];
 
      // Apply mcolptr change now that we have the iterators to the
      // data and row indices.
      (*mcolptr)[i] = colptr_mod;
 
      if (v[i]!=0)
        {
          // Check if r exist within this column, and get iterator for
          // mdata
          vector<Index>::iterator rpos = find(rstart,rend,r);
          vector<Numeric>::iterator dpos = dstart+(rpos-rstart);
          if (rpos!=rend)
            {
              // The index was found, replace the value in mdata with the
              // value from v.
              *dpos = v[i];
 
              // Copy this column to the ouput vectors.
              copy(rstart,rend,mrowind_it);
              copy(dstart,dend,mdata_it);
 
              // Adjust the position iterators accordingly.
              mrowind_it += rend-rstart;
              mdata_it += dend-dstart;
 
              // Set the mcolptr step, for next loop
              colptr_mod = (*mcolptr)[i]+(rend-rstart);
            }
          else
            {
              // The row index was not found, look for the first index
              // greater than r
              rpos = find_if(rstart,rend,bind2nd(greater<Index>(),r));
              dpos = dstart+(rpos-rstart);
 
              // Copy the first part of the column to the output vector.
              copy(rstart,rpos,mrowind_it);
              copy(dstart,dpos,mdata_it);
 
              // Make sure mrowind_it and mdata_it points at the first
              // 'empty' position.
              mrowind_it += rpos-rstart;
              mdata_it += dpos-dstart;
 
              // Insert the new value from v in mdata and the row index r in
              // mrowind.
              *mrowind_it = r;
              *mdata_it = v[i];
 
              // Again, make sure mrowind_it and mdata_it points at the
              // first 'empty' position.
              mrowind_it++;
              mdata_it++;
 
               // Copy the rest of this column to the output vectors.
              copy(rpos,rend,mrowind_it);
              copy(dpos,dend,mdata_it);
 
              // Adjust the iterators a last time
              mrowind_it += rend-rpos;
              mdata_it += dend-dpos;
 
              // Set the mcolptr step, for next loop
              colptr_mod = (*mcolptr)[i]+(rend-rstart+1);
            }
        }
      else
        {
          // Check if r exist within this column, and get iterator for
          // mdata
          vector<Index>::iterator rpos = find(rstart,rend,r);
          vector<Numeric>::iterator dpos = dstart+(rpos-rstart);
          if (rpos!=rend)
            {
              // The index was found and we use remove_copy to copy the rest
              // of the column to mrowind_new and mdata_new.
              remove_copy(rstart,rend,mrowind_it, *rpos);
              remove_copy(dstart,dend,mdata_it, *dpos);
 
              // Increase the mrowind_it and mdata_it
              mrowind_it += rend-rstart-1;
              mdata_it += dend-dstart-1;
 
              // Set the mcolptr step, for next loop
              colptr_mod = (*mcolptr)[i]+(rend-rstart-1);
            }
          else
            {
              // The row index was not found, all we need is to copy this
              // columns to the output mrowind and mdata vectors.
              copy(rstart,rend,mrowind_it);
              copy(dstart,dend,mdata_it);
 
              // Adjust the mrowind_it and mdata_it iterators
              mrowind_it += rend-rstart;
              mdata_it += dend-dstart;
 
              // Set the mcolptr step, for next loop
              colptr_mod = (*mcolptr)[i]+(rend-rstart);
            }
        }
    }
  // Apply mcolptr change for the one extra mcolptr element.
  *(mcolptr->end()-1) = colptr_mod;
 
  // Clean up?
  //delete mrowind_ref;
  //delete mdata_ref;
}
 
 
void Sparse::make_I( Index r, Index c)
{
  assert( 0<=r );
  assert( 0<=c );
 
  // First get number of ones in the identity matrix
  Index n = min(r,c);
 
  // Remake and assign values to vectors
  delete mcolptr;
  mcolptr = new vector<Index>( c+1, n);
  delete mrowind;
  mrowind = new vector<Index>(n);
  delete mdata;
  mdata = new vector<Numeric>(n,1.0);
 
  // Loop over number of ones and assign values [0:n-1] to mcolptr and
  // mrowind
  vector<Index>::iterator rit = mrowind->begin();
  vector<Index>::iterator cit = mcolptr->begin();
  for( Index i=0; i<n; i++ ) {
    *rit++ = i;
    *cit++ = i;
  }
 
  mrr = r;
  mcr = c;
}
 
 
void Sparse::resize(Index r, Index c)
{
  assert( 0<=r );
  assert( 0<=c );
  if ( mrr!=r || mcr!=c )
    {
      delete mdata;
      mdata = new vector<Numeric>;
      delete mrowind;
      mrowind = new vector<Index>;
      delete mcolptr;
      mcolptr = new vector<Index>(c+1,0);
 
      mrr = r;
      mcr = c;
    }
}
 
 
Sparse& Sparse::operator=(const Sparse& m)
{
  // It is important that we first delete previous content of the
  // target, otherwise we create a memory leak!
  //
  // The scalar delete operator is the correct one to use here, since
  // only one vector object has been allocated by new in each
  // case. (That the object itself is a vector does not matter.)
  //
  // Caveat: We should delete only if there really *is* previous
  // content, otherwise we will generate a core dump.
  if (mdata!=NULL)
    {
      assert (mrowind!=NULL);
      assert (mcolptr!=NULL);
 
      delete mdata;
      delete mrowind;
      delete mcolptr;
 
      mdata   = NULL;
      mrowind = NULL;
      mcolptr = NULL;
    }
 
  mrr = m.mrr;
  mcr = m.mcr;
 
  if (m.mdata!=NULL)
    {
      mdata   = new vector<Numeric>(*m.mdata);
      mrowind = new vector<Index>(*m.mrowind);
      mcolptr = new vector<Index>(*m.mcolptr);
    }
 
  return *this;
}
 
 
ostream& operator<<(ostream& os, const Sparse& v)
{
  for (size_t c=0; c<v.mcolptr->size()-1; ++c)
    {
      // Get index of first data element of this column:
      Index begin = (*v.mcolptr)[c];
 
      // Get index of first data element of next column. (This always works,
      // because mcolptr has one extra element pointing behind the last
      // column.)
      const Index end = (*v.mcolptr)[c+1];
 
      // Loop through the elements in this column:
      for ( Index i=begin; i<end; ++i )
        {
          // Get row index:
          Index r = (*v.mrowind)[i];
 
          // Output everything:
          os << setw(3) << r << " "
             << setw(3) << c << " "
             << setw(3) << (*v.mdata)[i] << "\n";
        }
    }
 
  return os;
}
 
// General matrix functions
 
 
void abs(       Sparse& A,
          const Sparse& B )
{
  // Check dimensions
  assert( A.nrows() == B.nrows() );
  assert( A.ncols() == B.ncols() );
 
  // Here we allocate memory for the A.mdata vector so that it matches the
  // input matrix, and then store the absolute values of B.mdata in it.
  A.mdata->resize( B.mdata->size() );
  Index end = B.mdata->size();
  for (Index i=0; i<end; i++)
    {
      (*A.mdata)[i] = fabs((*B.mdata)[i]);
    }
  
  // The column pointer and row index vectors are copies of the input  
  A.mcolptr->resize( B.mcolptr->size() );
  copy(B.mcolptr->begin(), B.mcolptr->end(), A.mcolptr->begin());
  A.mrowind->resize( B.mrowind->size() );
  copy(B.mrowind->begin(), B.mrowind->end(), A.mrowind->begin());
  
}
 
 
void mult( VectorView y,
           const Sparse& M,
           ConstVectorView x )
{
  // Check dimensions:
  assert( y.nelem() == M.nrows() );
  assert( M.ncols() == x.nelem() );
 
  // Initialize y to all zeros:
  y = 0.0;
 
  // Looping through every element of M, multiplying it with the
  // appropriate elements of x.
  for (size_t c=0; c<M.mcolptr->size()-1; ++c)
    {
      // Get index of first data element of this column:
      Index begin = (*M.mcolptr)[c];
 
      // Get index of first data element of next column. (This always works,
      // because mcolptr has one extra element pointing behind the last
      // column.)
      const Index end = (*M.mcolptr)[c+1];
 
      // Loop through the elements in this column:
      for ( Index i=begin; i<end; ++i )
        {
          // Get row index:
          Index r = (*M.mrowind)[i];
 
          // Compute this element:
          y[r] += (*M.mdata)[i] * x[c];
        }
    }
}
 
 
void mult( MatrixView A,
           const Sparse& B,
           ConstMatrixView C )
{
  // Check dimensions:
  assert( A.nrows() == B.nrows() );
  assert( A.ncols() == C.ncols() );
  assert( B.ncols() == C.nrows() );
 
  // Set all elements of A to zero:
  A = 0.0;
 
  // Loop through the elements of B:
  for (size_t c=0; c<B.mcolptr->size()-1; ++c)
    {
      // Get index of first data element of this column:
      Index begin = (*B.mcolptr)[c];
 
      // Get index of first data element of next column. (This always works,
      // because mcolptr has one extra element pointing behind the last
      // column.)
      const Index end = (*B.mcolptr)[c+1];
 
      // Loop through the elements in this column:
      for ( Index i=begin; i<end; ++i )
        {
          // Get row index:
          Index r = (*B.mrowind)[i];
 
          // Multiply this element with the corresponding row of C and
          // add the product to the right row of A
          for (Index j=0; j<C.ncols(); j++)
            {
              /* Conceptually:
                 A(r,j) += B.ro(r,c) * C(c,j);
                 But we don't need to use the index operator, because
                 we have the right element right here: */
              A(r,j) += (*B.mdata)[i] * C(c,j);
            }
        }
    }
}
 
 
 
void transpose( Sparse& A,
                const Sparse& B )
{
  // Check dimensions
  assert( A.nrows() == B.ncols() );
  assert( A.ncols() == B.nrows() );
  if ( B.mdata->size() == 0) return;
 
  // Here we allocate memory for the A.mdata and A.mrowind vectors so that
  // it matches the input matrix. (NB: that mdata and mrowind already has
  // the same size.)
  A.mdata->resize( B.mdata->size() );
  A.mrowind->resize( B.mrowind->size() );
 
  // Find the minimum and maximum existing row number in B.
  // (This maybe unnecessary in our case since we mostly treat diagonally
  // banded matrices where the minimum and maximum existing row numbers
  // coincide with the border numbers of the matrix.)
  vector<Index>::iterator startrow, stoprow;
  startrow = min_element( B.mrowind->begin(), B.mrowind->end() );
  stoprow = max_element( B.mrowind->begin(), B.mrowind->end() );
 
  // To create the A.mcolptr vector we loop through the existing row
  // numbers of B and count how many there are in B.mrowind.
  // First we make sure it is initialized to all zeros.
  A.mcolptr->assign( A.mcr+1, 0 );
  for (Index i=*startrow; i<=*stoprow; i++)
    {
      Index n = count( B.mrowind->begin(), B.mrowind->end(), i);
 
      // The column pointer for the column above the current one in A
      // are then set to the current one plus this amount.
      (*A.mcolptr)[i+1] = (*A.mcolptr)[i] + n;
    }
 
  // Next we loop through the columns of A and search for the corresponding
  // row index within each column of B (i.e. two loops). The elements are
  // then stored in A. We know that for every column in B we will have
  // the corresponding rowindex in A.
  vector<Numeric>::iterator dataA_it = A.mdata->begin();
  vector<Index>::iterator rowindA_it = A.mrowind->begin();
 
  // Looping over columns in A to keep track of what row number we are
  // looking for.
  for (size_t c=0; c<A.mcolptr->size(); ++c)
    {
      // Looping over columns in B to get the elements in the right order,
      // this will turn into row index in A.
      for (size_t i=0; i<B.mcolptr->size()-1; i++)
        {
          // Since only one element can occupy one row in a column, we
          // only need to call find() once per column in B.
          vector<Index>::iterator elem_it =
            find(B.mrowind->begin()+*(B.mcolptr->begin()+i),
                 B.mrowind->begin()+*(B.mcolptr->begin()+i+1), Index (c));
 
          // If we found the element, store it in A.mrowind and A.mdata
          if (elem_it != B.mrowind->begin()+*(B.mcolptr->begin()+i+1))
            {
              *rowindA_it = i;
              rowindA_it++;
 
              // To get the corresponding element in B.mdata we subtract
              // the initial value of the row index iterator from elem_it
              // and add the difference to a mdata iterator.
              Index diff = elem_it - B.mrowind->begin();
              vector<Numeric>::iterator elemdata_it=B.mdata->begin() + diff;
              *dataA_it = *elemdata_it;
              dataA_it++;
            }
        }
    }
}
 
 
 
void mult( Sparse& A,
           const Sparse& B,
           const Sparse& C )
{
  // Check dimensions and make sure that A is empty
  assert( A.nrows() == B.nrows() );
  assert( A.ncols() == C.ncols() );
  assert( B.ncols() == C.nrows() );
  A.mcolptr->assign( A.mcr+1, 0);
  A.mrowind->clear();
  A.mdata->clear();
 
  // Transpose B to simplify multiplication algorithm, after transposing we
  // can extract columns form the two matrices and multiply them, (which is
  // easier than extacting rows.)
  Sparse Bt(B.ncols(), B.nrows());
  //  cout << "B.nrows, B.ncols = " << B.nrows() << ", " << B.ncols() << "\n";
  //  cout << "B = \n" << B << "\n";
  transpose(Bt,B);
 
  // By looping over columns in C and multiply them with every column in Bt
  // (instead of the conventional loooping over Bt), we get the output
  // elements in the right order for storing them.
  for (size_t c=0; c<C.mcolptr->size()-1; ++c)
    {
      // Get row indices of this column
      //Index beginC = (*C.mcolptr)[c];
      //Index endC = (*C.mcolptr)[c+1];
 
      for (size_t b=0; b<Bt.mcolptr->size()-1; ++b)
        {
          // Get the intersection between the elements in the two columns and
          // and store them in a temporary vector
          set<Index> colintersec;
          set_intersection(C.mrowind->begin()+*(C.mcolptr->begin()+c),
            C.mrowind->begin()+*(C.mcolptr->begin()+c+1),
            Bt.mrowind->begin()+*(Bt.mcolptr->begin()+b),
            Bt.mrowind->begin()+*(Bt.mcolptr->begin()+b+1),
            inserter(colintersec, colintersec.begin()));
 
          // If we got an intersection, loop through it and multiply the
          // element pairs from C and Bt and store result in A
          if (!colintersec.empty())
            {
              Numeric tempA = 0.0;
              for (set<Index>::iterator i=colintersec.begin();
                i!=colintersec.end(); ++i)
                {
                  // To get iterators to the data elements in Bt and C, we
                  // subtract the iterator that points to the start of the
                  // mrowind vector from the iterators that points to the
                  // values in colintersec. We then get the number of steps 
                  // from the start of the both vectors mrowind and mdata to
                  // the intersecting values and can add them to the
                  // iterator that points to the beginning of mdata.
                  vector<Index>::iterator rowindBt_it =
                    find(Bt.mrowind->begin()+*(Bt.mcolptr->begin()+b),
                      Bt.mrowind->begin()+*(Bt.mcolptr->begin()+b+1), *i);
                  vector<Index>::iterator rowindC_it =
                    find(C.mrowind->begin()+*(C.mcolptr->begin()+c),
                      C.mrowind->begin()+*(C.mcolptr->begin()+c+1), *i);
 
                  vector<Numeric>::iterator dataBt_it =
                    Bt.mdata->begin()+(rowindBt_it-Bt.mrowind->begin());
                  vector<Numeric>::iterator dataC_it =
                    C.mdata->begin()+(rowindC_it-C.mrowind->begin());
 
                  tempA += *dataBt_it * *dataC_it;
                }
              A.rw(b,c) = tempA;
            }
        }
    }
}
 
 
 
void add( Sparse& A,
          const Sparse& B,
          const Sparse& C )
{
  // Check dimensions
  assert( B.ncols() == C.ncols() );
  assert( B.nrows() == C.nrows() );
 
  // We copy the smaller matrix of B or C to A. This way we can
  // loop over the matrix with the fewer number of elements to perform
  // the actual addition later. 
  const Sparse* D;
  if (B.data()->size() < C.data()->size())
    {
      A=C;
      D=&B;
    }
  else
    {
      A=B;
      D=&C;
    }
 
  for (size_t c = 0; c < D->mcolptr->size()-1; ++c)
    {
      // Loop through the elements in this column:
      for (Index i = (*D->mcolptr)[c]; i < (*D->mcolptr)[c+1]; ++i)
        {
          A.rw((*D->mrowind)[i], c) += (*D->mdata)[i];
        }
    }
}
 
 
 
void sub( Sparse& A,
          const Sparse& B,
          const Sparse& C )
{
  // Check dimensions
  assert( B.ncols() == C.ncols() );
  assert( B.nrows() == C.nrows() );
 
  A=B;
 
  for (size_t c = 0; c < C.mcolptr->size()-1; ++c)
    {
      // Loop through the elements in this column:
      for (Index i = (*C.mcolptr)[c]; i < (*C.mcolptr)[c+1]; ++i)
        {
          A.rw((*C.mrowind)[i], c) -= (*C.mdata)[i];
        }
    }
}