/*
 * ArrayClass.hxx
 * --------------
 *    file created:   January   16, 2013
 *    last modified:  March     24, 2016
 *    author:         Matthew MacLean
 *
 *    purpose:        define classes to simulate multi-dimensional dynamic arrays
 *
 * notable revisions:
 * -----------------
 *    1.4.0     - minor cleanup for web posting
 *    1.3.0     - base class for inheritance of operator() and indexing accessors
 *    1.2.0     - changed =operator to copy function to avoid massive
 *                memory allocations for mismatched arrays
 *    1.1.0     - add exception support & namespace confinement
 *    1.0.0     - first "official" version
 *
 *  ------------------------------------------------------------------------------
 *  ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
 *
 *  The contents of this file are released under the terms of either the Mozilla
 *  Public License or the GNU Lesser Public License as defined below.  You may
 *  indicate your preference by deleting the block of text below that applies to
 *  the other license.  If you do not delete either license, any recipient may
 *  use your versions of this file under either license also.
 *
 *  ------------------------------------------------------------------------------
 *  ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
 *
 *  
 *  This Source Code Form is subject to the terms of the Mozilla Public License,
 *  v. 2.0. If a copy of the MPL was not distributed with this file, 
 *  You can obtain one at http://mozilla.org/MPL/2.0/.
 *
 *  Software distributed under the License is distributed on an "AS IS"
 *  basis, WITHOUT WARRANTY OF ANY KIND, either express or implied. See the
 *  License for the specific language governing rights and limitations
 *  under the License.
 *
 *  The Original Code is ArrayClass header.
 *
 *  The Initial Developer of the Original Code is Matthew MacLean.
 *  Portions created by Matthew MacLean are Copyright (C) 2013-2016 by Matthew MacLean.
 *  All Rights Reserved.
 *
 *  Contributor(s):
 *    original code - Matthew MacLean
 *
 *  ------------------------------------------------------------------------------
 *  ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
 *
 *  This file is part of ArrayClass header.
 *  Copyright (C) 2013-2016 by Matthew MacLean.
 *
 *  This program is free software: you can redistribute it and/or modify
 *  it under the terms of the GNU Lesser General Public License as
 *  published by the Free Software Foundation, either version 3 of the
 *  License, or 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 Lesser General Public License for more details.
 *
 *  You should have received a copy of the GNU Lesser General Public License
 *  along with this program.  If not, see <http://www.gnu.org/licenses/>.
 *
 *  ------------------------------------------------------------------------------
 *  ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
 *
 */

#if !defined(_ARRAY_CLASS_DOT_H)
#define _ARRAY_CLASS_DOT_H

/* global includes */
#include <new>
#include <cstdlib>
#if defined (ARRAY_CLASS_BOUNDS_CHECK)
#  include <stdexcept>
#endif

/* namespace */
namespace ArrayClass	{


/*
 *  ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
 *  ------------------------------------------------------------------------------
 *  <<<                     PREPROCESSOR DEFINITIONS                           >>>
 *  ------------------------------------------------------------------------------
 *  ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
 */


//---------------------- DEFINED TYPES ----------------------

typedef unsigned int array_index;

//-------------------- END DEFINED TYPES --------------------


/*
 *  ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
 *  ------------------------------------------------------------------------------
 *  <<<                      PUBLIC CLASS DECLARATIONS                         >>>
 *  ------------------------------------------------------------------------------
 *  ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
 */



template <class data_type> class VectorBase
{
public:

	//-----must define operator for vector manipulation-----
	virtual data_type& operator() (array_index row) = 0;
	
	virtual data_type  operator() (array_index row) const = 0;
	
	//-----indexing accessor-----
	virtual array_index getRows() const = 0;

};




template <class data_type> class MatrixBase
{
public:

	//-----must define two-parameter operators for matrix manipulation-----
	virtual data_type& operator() (array_index row, array_index col) = 0;
	
	virtual data_type  operator() (array_index row, array_index col) const = 0;
	
	//-----indexing for these two indices-----
	virtual array_index getRows() const = 0;
	
	virtual array_index getCols() const = 0;

};




/*  ------------------------ <CLASS> Array1D <CLASS> -----------------------------
 *
 *  description:  1D Array class with operators, indexed as (R)
 *
 *  ------------------------------------------------------------------------------
 */
template <class data_type> class Array1D : public VectorBase<data_type>
{
public:
    //-----constructor-----
    Array1D(array_index row=1)
    {
        data = NULL;
        resize(row);
    }
    
    //-----copy constructor-----
    Array1D(const Array1D<data_type>& a)
    {
        data = NULL;
        this->resize(a.getRows());
        copy(a);
    }
    
    //-----destructor-----
    ~Array1D()
    {   delete [] data;     }
    
    //-----accessors-----
    inline array_index getRows() const
    {   return rows_val;    }
    
    inline size_t      getSizeOf() const
    {   return ( rows_val*sizeof(data_type) + sizeof(array_index) );    }
	
    //-----mutators-----
    inline void initialize(data_type val)
    {  (*this) = val;	}
    
    inline void resize(array_index row)
    {
        if (data != NULL) delete [] data;
        rows_val = row;
        if (rows_val <= 0) rows_val = 1;
        data = new data_type[rows_val];
    }
    
    inline bool copy(const Array1D<data_type>& a, bool allow_resize=false)
    {
        if (this->getRows() != a.getRows()) {
            if (allow_resize) this->resize(a.getRows());
            else return false;
        }
        
        for (array_index i=0; i<rows_val; i++)
            (*this)(i) = a(i);
        
        return true;
    }
    
    //-----subscript operators-----
#if defined (ARRAY_CLASS_BOUNDS_CHECK)
    inline data_type& operator() (array_index row)
    {     if (row>=rows_val) throw std::out_of_range("out of bounds access in Array1D:operator(row)"); else return *(data + row);  }
    
    inline data_type  operator() (array_index row) const
    {     if (row>=rows_val) throw std::out_of_range("out of bounds access in Array1D:operator(row) const"); else return *(data + row);  }
#else
    inline data_type& operator() (array_index row)
    {     return *(data + row);  }

    inline data_type  operator() (array_index row) const
    {     return *(data + row);  }
#endif

    //-----assignment operator-----    
    inline Array1D<data_type>& operator= (const data_type& a)
    {
        for (array_index i=0; i<rows_val; i++)
            (*this)(i) = a;
        
        return *this;
    }
   
private:
    //-----set private assignment operator so it cannot happen-----
    Array1D<data_type>& operator= (const Array1D<data_type>& a);
    
protected:
    array_index rows_val;
    data_type   *data;

};
/* ---------------------------<END> Array1D <END>--------------------------------- */



/*  ------------------------ <CLASS> Array2D <CLASS> -----------------------------
 *
 *  description:  2D Array class with operators, indexed as (R, C)
 *
 *  ------------------------------------------------------------------------------
 */
template <class data_type> class Array2D : public MatrixBase<data_type>
{
public:
    //-----constructor-----
    Array2D(array_index row=1, array_index col=1)
    {
        data = NULL;
        resize(row, col);
    }
    
    //-----copy constructor-----
    Array2D(const Array2D<data_type>& a)
    {
        data = NULL;
        this->resize(a.getRows(), a.getCols());
        copy(a);
    }
    
    //-----destructor-----
    ~Array2D()
    {   delete [] data; }
      
    //-----accessors-----
    inline array_index getRows() const
    {     return rows_val;   }
    
    inline array_index getCols() const
    {     return cols_val;   }
    
    inline size_t      getSizeOf() const
    {     return ( rows_val*cols_val*sizeof(data_type) + 2*sizeof(array_index) );    }
	
    //-----mutators-----
    inline void initialize(data_type val)
    {  (*this) = val;	}
    
    inline void resize(array_index row, array_index col)
    {
        if (data != NULL) delete [] data;
        rows_val = row;
        cols_val = col;
        if (rows_val <= 0) rows_val = 1;
        if (cols_val <= 0) cols_val = 1;
        data = new data_type[rows_val * cols_val];
    }
    
    inline bool copy(const Array2D<data_type>& a, bool allow_resize=false)
    {
        if ( (this->getRows() != a.getRows()) || (this->getCols() != a.getCols()) ) {
        	if (allow_resize) this->resize(a.getRows(), a.getCols());
            else return false;
        }
        
        for (array_index i=0; i<rows_val; i++)	{
            for (array_index j=0; j<cols_val; j++)
                (*this)(i,j) = a(i,j);
        }
        
        return true;
    }
    

    //-----subscript operators-----
#if defined (ARRAY_CLASS_BOUNDS_CHECK)
    inline data_type& operator() (array_index row, array_index col)
    {
        if ((row>=rows_val) || (col>=cols_val)) 
            throw std::out_of_range("out of bounds access in Array2D:operator(row,col)");
        else
            return *(data+row*cols_val + col);
    }
    
    inline data_type  operator() (array_index row, array_index col) const
    {
        if ((row>=rows_val) || (col>=cols_val))
            throw std::out_of_range("out of bounds access in Array2D:operator(row,col) const");
        else
            return *(data+row*cols_val + col);
    }
#else
    inline data_type& operator() (array_index row, array_index col)
    {     return *(data+row*cols_val + col);  }
    
    inline data_type  operator() (array_index row, array_index col) const
    {     return *(data+row*cols_val + col);  }
#endif


    //-----assignment operator-----
    inline Array2D<data_type>& operator= (const data_type& a)
    {
        for (array_index i=0; i<rows_val; i++)	{
            for (array_index j=0; j<cols_val; j++)
                (*this)(i,j) = a;
        }
        
        return *this;
    }

private:
    //-----set private assignment operator so it cannot happen-----
    Array2D<data_type>& operator= (const Array2D<data_type>& a);
    
protected:
    array_index rows_val, cols_val;
    data_type   *data;

};
/* ---------------------------<END> Array2D <END>--------------------------------- */



/*  ------------------------ <CLASS> Array3D <CLASS> -----------------------------
 *
 *  description:  3D Array class with operators, indexed as (I, J, K)
 *
 *  ------------------------------------------------------------------------------
 */
template <class data_type> class Array3D
{
public:
    //-----constructor-----
    Array3D(array_index I=1, array_index J=1, array_index K=1)
    {
        data = NULL;
        resize(I, J, K);
    }
    //-----copy constructor-----
    Array3D(const Array3D<data_type>& a)
    {
        data = NULL;
        this->resize(a.get_I_index(), a.get_J_Index(), a.get_K_index());
        copy(a);
    }
    
    //-----destructor-----
    ~Array3D()
    {
        delete [] data;
    }
      
    //-----accessors-----
    inline array_index get_I_index() const
    {     return I_val;   }
    
    inline array_index get_J_index() const
    {     return J_val;   }
    
    inline array_index get_K_index() const
    {     return K_val;   }
    
    inline size_t      getSizeOf() const
    {     return ( K_val*J_val*I_val*sizeof(data_type) + 5*sizeof(array_index) );    }
    
      
    //-----mutators-----
    inline void initialize(data_type val)
    {  (*this) = val;	}
    
    inline void resize(array_index I, array_index J, array_index K)
    {
        if (data != NULL) delete [] data;
        I_val = I;
        J_val = J;
        K_val = K;
        if (I_val <= 0) I_val = 1;
        if (J_val <= 0) J_val = 1;
        if (K_val <= 0)  K_val  = 1;
        K_mult = I_val*J_val;
        J_mult = I_val;
        data = new data_type[I_val * J_val * K_val];
    }
    
    inline bool copy (const Array3D<data_type>& a, bool allow_resize=false)
    {
        if ( (this->get_I_index() != a.get_I_index()) || (this->get_J_index() != a.get_J_index())  || (this->get_K_index() != a.get_K_index()) ) {
        	if (allow_resize) this->resize(a.get_I_index(), a.get_J_index(), a.get_K_index());
            else return false;
        }
        
        for (array_index i=0; i<I_val; i++)	{
            for (array_index j=0; j<J_val; j++)	{
                for (array_index k=0; k<K_val; k++)
                    (*this)(i,j,k) = a(i,j,k);
            }
        }
        
        return true;
    }


    //-----subscript operators-----
#if defined (ARRAY_CLASS_BOUNDS_CHECK)
    inline data_type& operator() (array_index i, array_index j, array_index k)
    {
        if ( (i>=I_val) || (j>=J_val) || (k>=K_val))
            throw std::out_of_range("out of bounds access in Array3D:operator(i,j,k)");
        else
            return *(data+k*K_mult + j*J_mult + i);
    }
    
    inline data_type  operator() (array_index i, array_index j, array_index k) const
    {
        if ( (i>=I_val) || (j>=J_val) || (k>=K_val)) 
            throw std::out_of_range("out of bounds access in Array3D:operator(i,j,k)");
        else
            return *(data+k*K_mult + j*J_mult + i);
    }
#else
    inline data_type& operator() (array_index i, array_index j, array_index k)
    {     return *(data+k*K_mult + j*J_mult + i);  }
    
    inline data_type  operator() (array_index i, array_index j, array_index k) const
    {     return *(data+k*K_mult + j*J_mult + i);  }
#endif


    //-----assignment operators-----
    inline Array3D<data_type>& operator= (const data_type& a)
    {
        for (array_index i=0; i<I_val; i++)	{
            for (array_index j=0; j<J_val; j++)	{
                for (array_index k=0; k<K_val; k++)
                    (*this)(i,j,k) = a;
            }
        }
        
        return *this;
    }

private:
    //-----set private assignment operator so it cannot happen-----
    Array3D<data_type>& operator= (const Array3D<data_type>& a);
    
protected:
    array_index I_val, J_val, K_val;
    array_index K_mult, J_mult;
    data_type   *data;

};
/* ---------------------------<END> Array3D <END>--------------------------------- */



/*  ------------------------ <CLASS> Array4D <CLASS> -----------------------------
 *
 *  description:  4D Array class with operators
 *
 *  ------------------------------------------------------------------------------
 */
template <class data_type> class Array4D : public VectorBase<data_type>
{
public:
    //-----constructor-----
    Array4D(array_index I=1, array_index J=1, array_index K=1, array_index ROW=1)
    {
        data = NULL;
        resize(I, J, K, ROW);
    }
    //-----copy constructor-----
    Array4D(const Array4D<data_type>& a)
    {
        data = NULL;
        this->resize(a.get_I_index(), a.get_J_index(), a.get_K_index(), a.getRows());
        copy(a);
    }
    
    //-----destructor-----
    ~Array4D()
    {
        location = NULL;
        delete [] data;
    }
      
    //-----accessors-----
    inline array_index getRows() const
    {     return rows_val;   }
    
    inline array_index get_I_index() const
    {     return I_val;   }
    
    inline array_index get_J_index() const
    {     return J_val;   }
    
    inline array_index get_K_index() const
    {     return K_val;   }
    
    inline size_t      getSizeOf() const
    {     return ( K_val*J_val*I_val*rows_val*sizeof(data_type) + 7*sizeof(array_index) );    }
    
      
    //-----mutators-----
    inline void setLocation(array_index i = 0, array_index j = 0, array_index k = 0)
    {
#if defined (ARRAY_CLASS_BOUNDS_CHECK)
        if ( (i>=I_val) || (j>=J_val) || (k>=K_val)) 
            throw std::out_of_range("out of bounds access in Array4D::setLocation");
        else
            location = data + k*K_mult + j*J_mult + i*I_mult;
#else
        location = data + k*K_mult + j*J_mult + i*I_mult;
#endif
    }


    inline void initialize(data_type val)
    {   (*this) = val;   }


    inline void resize(array_index I, array_index J, array_index K, array_index Row)
    {
        if (data != NULL) delete [] data;
        I_val    = I;
        J_val    = J;
        K_val    = K;
        rows_val = Row;
        if (I_val <= 0)     I_val = 1;
        if (J_val <= 0)     J_val = 1;
        if (K_val  <= 0)    K_val = 1;
        if (rows_val  <= 0) rows_val = 1;
        K_mult = I_val*J_val*rows_val;
        J_mult = I_val*rows_val;
	  I_mult = rows_val;
        data = new data_type[I_val * J_val * K_val * rows_val];
        location = data;
    }
    
    
    inline bool copy(const Array4D<data_type>& a, bool allow_resize=false)
    {
        if ( (this->get_I_index() != a.get_I_index()) || (this->get_J_index() != a.get_J_index())  || 
	       (this->get_K_index() != a.get_K_index()) || (this->getRows() != a.getRows()) ) {
        	if (allow_resize) this->resize(a.get_I_index(), a.get_J_index(), a.get_K_index(), a.getRows());
            else return false;
        }
        
        for (array_index i=0; i<I_val; i++) {
            for (array_index j=0; j<J_val; j++) {
                for (array_index k=0; k<K_val; k++)	{
                    for (array_index n=0; n<rows_val; n++)
                        (*this)(i,j,k,n) = a(i,j,k,n);
                }
            }
        }
        
        return true;
    }


    //-----subscript operators-----
#if defined (ARRAY_CLASS_BOUNDS_CHECK)
    inline data_type& operator() (array_index i, array_index j, array_index k, array_index row)
    {
        if ((i>=I_val) || (j>=J_val) || (k>=K_val) || (row>=rows_val))
            throw std::out_of_range("out of bounds access in Array4D:operator(i,j,k,row)");
        else
            return *(data+k*K_mult + j*J_mult + i*I_mult + row);
    }
    
    inline data_type  operator() (array_index i, array_index j, array_index k, array_index row) const
    {
        if ((i>=I_val) || (j>=J_val) || (k>=K_val) || (row>=rows_val))
            throw std::out_of_range("out of bounds access in Array4D:operator(i,j,k,row) const");
        else
            return *(data+k*K_mult + j*J_mult + i*I_mult + row);
    }

    inline data_type& operator() (array_index row)
    {
        if ((row>=rows_val))
            throw std::out_of_range("out of bounds access in Array4D:operator(row)");
        else
            return *(location + row);
    }

    inline data_type  operator() (array_index row) const
    {
        if ((row>=rows_val))
            throw std::out_of_range("out of bounds access in Array4D:operator(row) const");
        else
            return *(location + row);
    }
#else
    inline data_type& operator() (array_index i, array_index j, array_index k, array_index row)
    {     return *(data+k*K_mult + j*J_mult + i*I_mult + row);  }

    inline data_type  operator() (array_index i, array_index j, array_index k, array_index row) const
    {     return *(data+k*K_mult + j*J_mult + i*I_mult + row);  }

    inline data_type& operator() (array_index row)
    {     return *(location + row);  }

    inline data_type  operator() (array_index row) const
    {     return *(location + row);  }
#endif


    //-----assignment operators-----
    inline Array4D<data_type>& operator= (const data_type& a)
    {
        for (array_index i=0; i<I_val; i++) {
            for (array_index j=0; j<J_val; j++) {
                for (array_index k=0; k<K_val; k++)	{
                    for (array_index n=0; n<rows_val; n++)
                        (*this)(i,j,k,n) = a;
                }
            }
        }
        
        return *this;
    }
    
private:
    //-----set private assignment operator so it cannot happen-----
    Array4D<data_type>& operator= (const Array4D<data_type>& a);

protected:
    array_index I_val, J_val, K_val, rows_val;
    array_index K_mult, J_mult, I_mult;
    data_type   *data, *location;

};
/* ---------------------------<END> Array4D <END>--------------------------------- */



/*  ------------------------ <CLASS> Array5D <CLASS> -----------------------------
 *
 *  description:  5D Array class with operators
 *
 *  ------------------------------------------------------------------------------
 */
template <class data_type> class Array5D : public MatrixBase<data_type>
{
public:
    //-----constructor-----
    Array5D(array_index I=1, array_index J=1, array_index K=1, array_index ROW=1, array_index COL=1)
    {
        data = NULL;
        resize(I, J, K, ROW, COL);
    }
    //-----copy constructor-----
    Array5D(const Array5D<data_type>& a)
    {
        data = NULL;
        this->resize(a.get_I_index(), a.get_J_index(), a.get_K_index(), a.getRows(), a.getCols());
        copy(a);
    }
    
    //-----destructor-----
    ~Array5D()
    {
        location = NULL;
        delete [] data;
    }
      
    //-----accessors-----
    inline array_index getRows() const
    {     return rows_val;   }
    
    inline array_index getCols() const
    {     return cols_val;   }
    
    inline array_index get_I_index() const
    {     return I_val;   }
    
    inline array_index get_J_index() const
    {     return J_val;   }
    
    inline array_index get_K_index() const
    {     return K_val;   }
    
    inline size_t      getSizeOf() const
    {     return ( K_val*J_val*I_val*rows_val*cols_val*sizeof(data_type) + 9*sizeof(array_index) );    }
    
      
    //-----mutators-----
    inline void setLocation(array_index i = 0, array_index j = 0, array_index k = 0)
    {
#if defined (ARRAY_CLASS_BOUNDS_CHECK)
        if ( (i>=I_val) || (j>=J_val)  || (k>=K_val) ) 
            throw std::out_of_range("out of bounds access in Array5D::setLocation");
        else
            location = data + k*K_mult + j*J_mult + i*I_mult;
#else
        location = data + k*K_mult + j*J_mult + i*I_mult;
#endif
    }


    inline void initialize(data_type val)
    {   (*this) = val;   }


    inline void resize(array_index I, array_index J, array_index K, array_index Row, array_index Col)
    {
        if (data != NULL) delete [] data;
        I_val    = I;
        J_val    = J;
        K_val    = K;
        rows_val = Row;
	  cols_val = Col;
        if (I_val <= 0)     I_val = 1;
        if (J_val <= 0)     J_val = 1;
        if (K_val  <= 0)    K_val = 1;
        if (rows_val  <= 0) rows_val = 1;
	  if (cols_val  <= 0) cols_val = 1;
        K_mult   = I_val*J_val*rows_val*cols_val;
        J_mult   = I_val*rows_val*cols_val;
	  I_mult   = rows_val*cols_val;
	  row_mult = cols_val;
        data = new data_type[I_val * J_val * K_val * rows_val * cols_val];
        location = data;
    }
    
    
    inline bool copy(const Array5D<data_type>& a, bool allow_resize=false)
    {
        if ( (this->get_I_index() != a.get_I_index()) || (this->get_J_index() != a.get_J_index())  || 
	       (this->get_K_index() != a.get_K_index()) || (this->getRows() != a.getRows()) ||
		 (this->getCols() != a.getCols()) ) {
        	if (allow_resize) this->resize(a.get_I_index(), a.get_J_index(), a.get_K_index(), a.getRows(), a.getCols());
            else return false;
        }
        
        for (array_index i=0; i<I_val; i++) {
            for (array_index j=0; j<J_val; j++) {
                for (array_index k=0; k<K_val; k++)	{
                    for (array_index m=0; m<rows_val; m++)	{
                        for (array_index n=0; n<cols_val; n++)
				    (*this)(i,j,k,m,n) = a(i,j,k,m,n);
                    }
                }
            }
        }
        
        return true;
    }


    //-----subscript operators-----
#if defined (ARRAY_CLASS_BOUNDS_CHECK)
    inline data_type& operator() (array_index i, array_index j, array_index k, array_index row, array_index col)
    {
        if ( (i>=I_val) || (j>=J_val) || (k>=K_val) || (row>=rows_val) || (col>=cols_val))
            throw std::out_of_range("out of bounds access in Array5D:operator(i,j,k,row,col)");
        else
            return *(data+k*K_mult + j*J_mult + i*I_mult + row*row_mult + col);
    }
    
    inline data_type  operator() (array_index i, array_index j, array_index k, array_index row, array_index col) const
    {
        if ( (i>=I_val) || (j>=J_val) || (k>=K_val) || (row>=rows_val) || (col>=cols_val)) 
            throw std::out_of_range("out of bounds access in Array5D:operator(i,j,k,row,col) const");
        else
            return *(data+k*K_mult + j*J_mult + i*I_mult + row*row_mult + col);
    }

    inline data_type& operator() (array_index row, array_index col)
    {
        if ((row>=rows_val) || (col>=cols_val))
            throw std::out_of_range("out of bounds access in Array5D:operator(row,col)");
        else
            return *(location + row*row_mult + col);
    }

    inline data_type  operator() (array_index row, array_index col) const
    {
        if ((row>=rows_val) || (col>=cols_val))
            throw std::out_of_range("out of bounds access in Array5D:operator(row,col) const");
        else
            return *(location + row*row_mult + col);
    }
#else
    inline data_type& operator() (array_index i, array_index j, array_index k, array_index row, array_index col)
    {     return *(data+k*K_mult + j*J_mult + i*I_mult + row*row_mult + col);  }

    inline data_type  operator() (array_index i, array_index j, array_index k, array_index row, array_index col) const
    {     return *(data+k*K_mult + j*J_mult + i*I_mult + row*row_mult + col);  }

    inline data_type& operator() (array_index row, array_index col)
    {     return *(location + row*row_mult + col);  }

    inline data_type  operator() (array_index row, array_index col) const
    {     return *(location + row*row_mult + col);  }
#endif


    //-----assignment operators-----
    inline Array5D<data_type>& operator= (const data_type& a)
    {
        for (array_index i=0; i<I_val; i++) {
            for (array_index j=0; j<J_val; j++) {
                for (array_index k=0; k<K_val; k++)	{
                    for (array_index m=0; m<rows_val; m++)	{
                        for (array_index n=0; n<cols_val; n++)
				    (*this)(i,j,k,m,n) = a;
                    }
                }
            }
        }
        
        return *this;
    }
    
private:
    //-----set private assignment operator so it cannot happen-----
    Array5D<data_type>& operator= (const Array5D<data_type>& a);

protected:
    array_index I_val, J_val, K_val, rows_val, cols_val;
    array_index K_mult, J_mult, I_mult, row_mult;
    data_type   *data, *location;

};
/* ---------------------------<END> Array5D <END>--------------------------------- */


} //close namespace

#endif   //  _ARRAY_CLASS_DOT_H