Workspace.h

Go to the documentation of this file.

// ========================================================================= //
// This file is part of MyraMath, copyright (c) 2014-2019 by Ryan A Chilton  //
// and distributed by MyraCore, LLC. See LICENSE.txt for license terms.      //
// ========================================================================= //

#ifndef MYRAMATH_UTILITY_WORKSPACE_H
#define MYRAMATH_UTILITY_WORKSPACE_H

// Serialization interface.
#include <myramath/io/Streams.h>
#include <myramath/io/ReaderWriter.h>

// Initialize with diabolical values? (just for testing)
#include <myramath/utility/random.h>

#include <cstdlib>
#include <iostream>

namespace myra {

// Implementation detail of member swap. Why not use std::swap? Trying to avoid dependence on #include<algorithm>
template<class POD> void swap_pod(POD& pod1, POD& pod2)
  {
  POD temp = pod1;
  pod1 = pod2;
  pod2 = temp;
  }

template<class T> class Workspace
  {
  public:

    Workspace() : p(0), S(0)
      { }
  
    explicit Workspace(size_t in_S) : p(0), S(0)
      { allocate(in_S); }

    Workspace(const Workspace& that) : p(0), S(0)
      {
      allocate(that.S);
      for (size_t s = 0; s < S; ++s)
        (*this)[s] = that[s];
      }

    void swap(Workspace& that)
      {
      swap_pod(this->S, that.S);
      swap_pod(this->p, that.p);
      }

    Workspace& operator = (Workspace that)
      { this->swap(that); return *this; }

    explicit Workspace(InputStream& in) : p(0), S(0)
      {
      allocate(in.read<size_t>());
      for (size_t s = 0; s < S; ++s)
        in >> p[s];
      }

    void write(OutputStream& out) const
      {
      out << S;
      for (size_t s = 0; s < S; ++s)
        out << p[s];
      }
   
    const T& operator () (size_t i) const { return p[i]; }
          T& operator () (size_t i)       { return p[i]; }

    const T& operator [] (size_t i) const { return p[i]; }
          T& operator [] (size_t i)       { return p[i]; }
    
    T* begin() { return p; }
    T* end()   { return p+S; }

    const T* begin() const { return p; }
    const T* end()   const { return p+S; }

    size_t size() const
      { return S; }

    void resize(size_t R)
      {
      deallocate(); 
      allocate(R);
      }

    void resize(size_t R, T t)
      {
      resize(R);
      for (size_t s = 0; s < S; ++s)
        (*this)[s] = t;
      }

    void touch()
      { 
      for (size_t s = 0; s < S; s += 1024)
        (*this)[s] = 0;
      }
    
   ~Workspace()
      { deallocate(); } 
     
  private:

    // Allocates with size S
    void allocate(size_t in_S)
      {    
      S = in_S;
      p = 0;
      if (S) 
        {
        p = (T*)malloc(sizeof(T)*S);
        if (!p)
          throw std::bad_alloc();
        }

      // To test/verify that clients don't rely on workspace 
      // being initialized, can populate p[] with junk data..
      // for (int s = 0; s < S; ++s) p[s] = random<T>();

      // ..however some BLAS/LAPACK implementations may actually 
      // require workspace to be initialized with valid numbers 
      // before first touch ***, so consider zeroing it.
      // for (int s = 0; s < S; ++s) p[s] = T(0);
      }

    // *** Note ***
    //
    // This remark is more subtle than it first appears. An unitialized
    // workspace might have any random bit pattern in it. When updating
    // a workspace that is known to be uninitialized, the first touch will
    // always be a careful BLAS write (like gemm, syrk, or herk) with 
    // alpha != 0 and beta = 0. However, if there is a Number c within an
    // unitialized workspace that unluckily has (for instance) some sort 
    // of NaN bit pattern, it cannot actually be cleared by an operation 
    // of the form c = alpha*a*b + beta*c! The pre-existing NaN persists, 
    // despite the fact that we tried to "clear" it by using beta=0.
    // 
    // Some BLAS's can tolerate input NaN's, some cannot. For instance,
    // MKL explicitly guarantees "when beta is equal to zero, then c 
    // need not be set on input" in the documentation for gemm(). MyraKL
    // is also tolerant of input NaN's. The reference BLAS is not.
    //
    // In fact, there are many bit patterns that all can represent a NaN, 
    // so it's not at all unusual to encounter one in unitialized memory 
    // carved from the heap.

    // Deallocates.
    void deallocate()
      { 
      free(p); 
      p = 0;
      }

    // Pointer to internal contents.
    T* p;

    // Total size.    
    size_t S;

  };

template<class T> std::ostream& operator << (std::ostream& out, const Workspace<T>& w)
  {
  out << "[ ";
  for (const T* i = w.begin(); i != w.end(); ++i)
    out << *i << " ";
  out << "] (" << w.size() << ")";
  return out;
  }

} // namespace myra

#endif