Matrix.h

/* -*- c++ -*- */

/* WARNING: don't really know yet how to properly use templates: is it ok to
 * have all the code in this file and not in a .cpp file? If it is in a .cpp
 * file it does not work: i guess the needed functions are generated 
 * at linker time and the linker does not seem to go searching for them
 * in a different file than the file needed?
 */

#ifndef __Matrix_h__
#define __Matrix_h__

#include <cassert>
#include <vector>
using std::vector;

#include <iostream>
using std::cout;
using std::cerr;
using std::endl;

#include <fstream>
using std::ofstream;

#include <cstdlib>

#include <assert.h>

template<typename T>
class Matrix {
public:
  // a default constructor
  Matrix();

  Matrix(unsigned nrows, unsigned ncols);
  // Throws a BadSize object if either size is zero
  class BadSize { };

  Matrix(T * pointer, unsigned nrows, unsigned ncols);
  
  // Based on the Law Of The Big Three:
  ~Matrix();
  Matrix(const Matrix<T>& m);
   Matrix<T>& operator= (const Matrix<T>& m);
 
   // Access methods to get the (i,j) element:
   T&       operator() (unsigned i, unsigned j);
   const T& operator() (unsigned i, unsigned j) const;
   // These throw a BoundsViolation object if i or j is too big
   //class BoundsViolation { };
 
   unsigned get_nrows() const;
   unsigned get_ncols() const;

   void print(const std::string & filename);
   void print(const std::string & filename, 
              const vector<float> & xvals, const vector<float> & yvals);

 private:
   unsigned nrows_, ncols_;
   T* data_;
   bool delete_flag;
 };


 
template<typename T>
inline unsigned Matrix<T>::get_nrows() const
{ return nrows_; }
 
template<typename T>
inline unsigned Matrix<T>::get_ncols() const
{ return ncols_; }

 template<typename T>
 inline T& Matrix<T>::operator() (unsigned row, unsigned col)
 {
   if (row >= nrows_ || col >= ncols_) 
     assert(1==0);
     //throw BoundsViolation();
   return data_[row*ncols_ + col];
 }
 
 template<typename T>
 inline const T& Matrix<T>::operator() (unsigned row, unsigned col) const
 {
   if (row >= nrows_ || col >= ncols_)
     assert(1==0);
     //throw BoundsViolation();
   return data_[row*ncols_ + col];
 }

 template<typename T>
 inline Matrix<T>::Matrix(T * pointer, unsigned nrows, unsigned ncols)
   : nrows_ (nrows)
   , ncols_ (ncols)
{
  data_ = pointer;
  delete_flag = false;
}

template<typename T>
inline Matrix<T>::Matrix()
{
  nrows_ = 1;
  ncols_ = 1;
  data_ = new T[nrows_ * ncols_];
  delete_flag = true;
}
 
 template<typename T>
 inline Matrix<T>::Matrix(unsigned nrows, unsigned ncols)
   : nrows_ (nrows)
   , ncols_ (ncols)
 //, data_ <--initialized below (after the 'if/throw' statement)
 {
   if (nrows == 0 || ncols == 0)
     {
       assert(1==0);
       throw BadSize();
     }
   data_ = new T[nrows * ncols];
   delete_flag = true;
 }
 
 template<typename T>
 inline Matrix<T>::~Matrix()
 {
   if (delete_flag)
     delete[] data_;
 }

template<typename T>
Matrix<T>& Matrix<T>::operator= (const Matrix<T>& m) 
{
  if (this != &m) {  // make sure not same object
    if (delete_flag)
      {
        /* default action */
        delete [] data_;                      // Delete old name's memory.
        data_ = new T[m.get_nrows() * m.get_ncols()]; // Get new space

        // S- added
        nrows_ = m.get_nrows();
        ncols_ = m.get_ncols();
      }
    else
      {
        /* bad code, for pgplot use */
        if (this->get_nrows() != m.get_nrows() || 
            this->get_ncols() != m.get_ncols())
          {
            cout << "Trying to copy a matrix to a matrix of a different size"
                 << ", constructed from a pointer" << endl
                 << "exiting" << endl;
            exit(-1);
          }
      }

    for (unsigned i = 0; i < m.get_nrows(); ++i)
      for (unsigned j = 0; j < m.get_nrows(); ++j)
        {
          (*this)(i,j) = m(i,j);
        }
  }
  return *this;    // Return ref for multiple assignment
}//end operator=

template<typename T>
void Matrix<T>::print (const std::string & filename)
{
  ofstream outdata; // outdata is like ci
  outdata.open(filename.c_str()); // opens the file

  if( !outdata ) { // file couldn't be opened
    cerr << "Error: file could not be opened" << endl;
    exit(1);
  }

  outdata << "# i j (*this)(j,i)" << endl;
  for (unsigned j = 0; j < get_nrows(); ++j)
    {
      for (unsigned i = 0; i < get_ncols(); ++i){
        outdata << i << " " << j << " " << (*this)(j,i) << endl;
      }

      outdata << endl;
    }
  
   outdata.close();
}


template<typename T>
void Matrix<T>::print (const std::string & filename, 
                       const vector<float> & xvals, const vector<float> & yvals)
{
  ofstream outdata;
  outdata.open(filename.c_str()); // opens the file

  if( !outdata ) { // file couldn't be opened
    cerr << "Error: file " << filename << " could not be opened" << endl;
    exit(1);
  }

  outdata << "# x y val (*this)(j (y),i (x))" << endl;
  for (unsigned j = 0; j < get_nrows(); ++j)
    {
      // output: j,i: j: row-index e.g. y, i: col-index e.g. x
      for (unsigned i = 0; i < get_ncols(); ++i){
        outdata << xvals[i] << "\t" << yvals[j] << "\t" << (*this)(j,i) << endl;
      }

      outdata << endl;
    }
  
   outdata.close();
}


#endif // __Matrix_h__

---