escript/test/DataMathsTestCase.cpp


/*******************************************************
*
* Copyright (c) 2003-2008 by University of Queensland
* Earth Systems Science Computational Center (ESSCC)
* http://www.uq.edu.au/esscc
*
* Primary Business: Queensland, Australia
* Licensed under the Open Software License version 3.0
* http://www.opensource.org/licenses/osl-3.0.php
*
*******************************************************/


#include "escript/DataAlgorithm.h"
#include "escript/DataVector.h"
#include "esysUtils/EsysException.h"

#include "DataMathsTestCase.h"
#include "escript/DataTypes.h"

#include <iostream>

using namespace CppUnitTest;
using namespace esysUtils;
using namespace escript;
using namespace std;
using namespace escript::DataTypes;
using namespace escript::DataMaths;

void DataMathsTestCase::setUp() {
  //
  // This is called before each test is run

}

void DataMathsTestCase::tearDown() {
  //
  // This is called after each test has been run

}


void DataMathsTestCase::testMatMult()
{

  {
    cout << endl;
    cout << "\tTest result shape." << endl;

    DataTypes::ShapeType leftShape;
    leftShape.push_back(1);
    leftShape.push_back(3);
    DataTypes::ValueType leftData(DataTypes::noValues(leftShape),0);
//     DataArrayView leftDataView(leftData,leftShape);

    DataTypes::ShapeType rightShape;
    rightShape.push_back(3);
    rightShape.push_back(2);
    DataTypes::ValueType rightData(DataTypes::noValues(rightShape),0);
//     DataArrayView rightDataView(rightData,rightShape);

    DataTypes::ShapeType resultShape=DataMaths::determineResultShape(leftShape,rightShape);

    assert(resultShape.size()==2);
    assert(resultShape[0]==1);
    assert(resultShape[1]==2);

    DataTypes::ValueType resultData(DataTypes::noValues(resultShape),0);

    cout << "\tTest matrix multiplication.";
    double aValue=0.0;
    for (int i=0;i<leftShape[0];i++) {
      for (int j=0;j<leftShape[1];j++) {
    leftData[getRelIndex(leftShape,i,j)]=++aValue;
      }
    }
    aValue=0.0;
    for (int i=0;i<rightShape[0];i++) {
      for (int j=0;j<rightShape[1];j++) {
    rightData[getRelIndex(rightShape,i,j)]=++aValue;
      }
    }

    DataMaths::matMult(leftData,leftShape,0,rightData,rightShape,0,resultData, resultShape);
    assert((resultData[0]==22) && (resultData[1]==28));
  }

  cout << endl;

}

void DataMathsTestCase::testUnaryOp()
{

  // This typedef allows function names to be cast to pointers
  // to unary functions of the appropriate type.
  typedef double (*UnaryDFunPtr)(double);

  {
    cout << endl;
    cout << "\tTest unaryOp on scalar Data.";

    // define the shape for the Data
    DataTypes::ShapeType shape;

    // allocate the data for the Data
    int npoints=4;
    DataTypes::ValueType data(DataTypes::noValues(shape)*npoints,0);

    double tmp;
    int offset=0;
    // step the view along each data point in the underlying data
    for (int p=0;p<npoints;p++) {

      // assign values to the data point
      data[offset]=p;

      // apply a unary operation to this data point
      unaryOp(data,shape,offset,(UnaryDFunPtr)std::sin);

      // check the results
      tmp = std::sin((double)p);
      assert(std::abs(data[offset]-tmp)<=REL_TOL*std::abs(tmp));

      if (p<npoints-1) {
        offset++;
      }

    }

  }

  {
    cout << endl;
    cout << "\tTest unaryOp on shape (2,3) Data.";

    // define the shape for the DataArrayView
    DataTypes::ShapeType shape;
    shape.push_back(2);
    shape.push_back(3);

    // allocate the data for the DataArrayView
    int npoints=4;
    DataTypes::ValueType data(DataTypes::noValues(shape)*npoints,0);


    int offset=0;
    // step the view along each data point in the underlying data
    for (int p=0;p<npoints;p++) {

      // assign values to the data point
      for (int i=0;i<shape[0];i++) {
        for (int j=0;j<shape[1];j++) {
          data[offset+getRelIndex(shape,i,j)]=offset+getRelIndex(shape,i,j);
        }
      }

      // apply a unary operation to this data point
//      dataView.unaryOp((UnaryDFunPtr)std::sqrt);
      unaryOp(data,shape,offset,(UnaryDFunPtr)std::sqrt);

      // check the results
      for (int i=0;i<shape[0];i++) {
        for (int j=0;j<shape[1];j++) {
//           assert(std::abs(dataView(i,j)-std::sqrt((double)dataView.index(i,j)))<=REL_TOL*std::sqrt((double)dataView.index(i,j)));
          assert(std::abs(data[offset+getRelIndex(shape,i,j)]-std::sqrt((double)offset+getRelIndex(shape,i,j)))<=REL_TOL*std::sqrt((double)offset+getRelIndex(shape,i,j)));
        }
      }

      if (p<npoints-1) {
        offset++;
      }

    }

  }

  {
    cout << endl;
    cout << "\tTest unaryOp on shape (9,8,5,11) Data.";

    // define the shape for the DataArrayView
    DataTypes::ShapeType shape;
    shape.push_back(9);
    shape.push_back(8);
    shape.push_back(5);
    shape.push_back(11);

    // allocate the data for the DataArrayView
    int npoints=4;
    DataTypes::ValueType data(DataTypes::noValues(shape)*npoints,0);

    int offset=0;
    // step the view along each data point in the underlying data
    for (int p=0;p<npoints;p++) {

      // assign values to the data point
      for (int i=0;i<shape[0];i++) {
        for (int j=0;j<shape[1];j++) {
          for (int k=0;k<shape[2];k++) {
            for (int l=0;l<shape[3];l++) {
              data[offset+getRelIndex(shape,i,j,k,l)]=data[offset+getRelIndex(shape,i,j,k,l)]+1;
            }
          }
        }
      }

      // apply a unary operation to this data point
//       dataView.unaryOp((UnaryDFunPtr)std::log);
      unaryOp(data,shape,offset,(UnaryDFunPtr)std::log);

      // check the results
      for (int i=0;i<shape[0];i++) {
        for (int j=0;j<shape[1];j++) {
          for (int k=0;k<shape[2];k++) {
            for (int l=0;l<shape[3];l++) {
              assert(std::abs(data[offset+getRelIndex(shape,i,j,k,l)]-std::log(1+(double)data[offset+getRelIndex(shape,i,j,k,l)]))<=REL_TOL*std::abs(std::log(1+(double)data[offset+getRelIndex(shape,i,j,k,l)])));
            }
          }
        }
      }

      if (p<npoints-1) {
        offset++;
      }

    }

  }

  cout << endl;

}

void DataMathsTestCase::testBinaryOp()
{

  // This typedef allows function names to be cast to pointers
  // to binary functions of the appropriate type.
  typedef double (*BinaryDFunPtr)(double,double);

  {
    cout << endl;
    cout << "\tTest binaryOp on scalar Data.";

    // define the shape for the DataArrayViews
    DataTypes::ShapeType shape;

    // allocate the data for the DataArrayViews
    int npoints=4;
    DataTypes::ValueType data1(DataTypes::noValues(shape)*npoints,0);
    DataTypes::ValueType data2(DataTypes::noValues(shape)*npoints,0);

    int offset=0;
    // step the views along each data point in the underlying data
    for (int p=0;p<npoints;p++) {

      // assign values to the data points
      data1[offset]=p;
      data2[offset]=p;

      // apply a binary operation to these data points
      binaryOp(data1,scalarShape,offset,data2,scalarShape,offset, plus<double>());

      // check the results
      assert(data1[offset]==p+p);

      if (p<npoints-1) {
    ++offset;
      }

    }

  }

  {
    cout << endl;
    cout << "\tTest binaryOp on shape (2,3) Data.";

    // define the shape for the DataArrayViews
    DataTypes::ShapeType shape;
    shape.push_back(2);
    shape.push_back(3);

    // allocate the data for the DataArrayViews
    int npoints=4;
    DataTypes::ValueType data1(DataTypes::noValues(shape)*npoints,0);
    DataTypes::ValueType data2(DataTypes::noValues(shape)*npoints,0);

    int offset=0;
    // step the views along each data point in the underlying data
    for (int p=0;p<npoints;p++) {

      // assign values to the data points
      for (int i=0;i<shape[0];i++) {
        for (int j=0;j<shape[1];j++) {
          data1[offset+getRelIndex(shape,i,j)]=offset+getRelIndex(shape,i,j);
          data2[offset+getRelIndex(shape,i,j)]=offset+getRelIndex(shape,i,j);
        }
      }

      // apply a binary operation to these data points
/*      dataView1.binaryOp(dataView2,multiplies<double>());*/
      binaryOp(data1,shape,offset,data2,shape,offset,multiplies<double>());

      // check the results
      for (int i=0;i<shape[0];i++) {
        for (int j=0;j<shape[1];j++) {
          assert(data1[offset+getRelIndex(shape,i,j)]==(offset+getRelIndex(shape,i,j))*(offset+getRelIndex(shape,i,j)));
        }
      }

      if (p<npoints-1) {
    offset+=noValues(shape);
      }

    }

  }

  {
    cout << endl;
    cout << "\tTest binaryOp on shape (9,8,5,11) Data.";

    // define the shape for the DataArrayViews
    DataTypes::ShapeType shape;
    shape.push_back(9);
    shape.push_back(8);
    shape.push_back(5);
    shape.push_back(11);

    // allocate the data for the DataArrayViews
    int npoints=4;
    DataTypes::ValueType data1(DataTypes::noValues(shape)*npoints,0);
    DataTypes::ValueType data2(DataTypes::noValues(shape)*npoints,0);

    int offset=0;
    // step the views along each data point in the underlying data
    for (int p=0;p<npoints;p++) {

      // assign values to the data points
      for (int i=0;i<shape[0];i++) {
        for (int j=0;j<shape[1];j++) {
          for (int k=0;k<shape[2];k++) {
            for (int l=0;l<shape[3];l++) {
              data1[offset+getRelIndex(shape,i,j,k,l)]=offset+getRelIndex(shape,i,j,k,l);
              data2[offset+getRelIndex(shape,i,j,k,l)]=offset+getRelIndex(shape,i,j,k,l);
            }
          }
        }
      }

      // apply a binary operation to these data points
//      dataView1.binaryOp(dataView2,multiplies<double>());
      binaryOp(data1,shape,offset,data2,shape,offset,multiplies<double>());

      // check the results
      for (int i=0;i<shape[0];i++) {
        for (int j=0;j<shape[1];j++) {
          for (int k=0;k<shape[2];k++) {
            for (int l=0;l<shape[3];l++) {
              assert(data1[offset+getRelIndex(shape,i,j,k,l)]==(offset+getRelIndex(shape,i,j,k,l))*(offset+getRelIndex(shape,i,j,k,l)));
            }
          }
        }
      }

      if (p<npoints-1) {
    offset+=noValues(shape);
      }

    }

  }

  {
    cout << endl;
    cout << "\tTest binaryOp on scalar Data and single value.";

    // define the shape for the DataArrayView
    DataTypes::ShapeType shape;

    // allocate the data for the DataArrayView
    int npoints=4;
    DataTypes::ValueType data(DataTypes::noValues(shape)*npoints,0);

    int offset=0;
    // step the view along each data point in the underlying data
    for (int p=0;p<npoints;p++) {

      // assign values to the data point
      data[offset]=p;

      // apply a binary operation to this data point
//       dataView.binaryOp(4.9,plus<double>());
      binaryOp(data,shape,offset,4.9,plus<double>());

      // check the results
      assert(data[offset]==4.9+p);

      if (p<npoints-1) {
        ++offset;
      }

    }

  }

  {
    cout << endl;
    cout << "\tTest binaryOp on shape (2,3) Data and single value.";

    // define the shape for the DataArrayView
    DataTypes::ShapeType shape;
    shape.push_back(2);
    shape.push_back(3);

    // allocate the data for the DataArrayView
    int npoints=4;
    DataTypes::ValueType data(DataTypes::noValues(shape)*npoints,0);

    int offset=0;
    // step the view along each data point in the underlying data
    for (int p=0;p<npoints;p++) {

      // assign values to the data point
      for (int i=0;i<shape[0];i++) {
        for (int j=0;j<shape[1];j++) {
          data[offset+getRelIndex(shape,i,j)]=offset+getRelIndex(shape,i,j);
        }
      }

      // apply a binary operation to the data point
//       dataView.binaryOp(5.8,multiplies<double>());
      binaryOp(data,shape,offset,5.8,multiplies<double>());

      double tmp;
      // check the results
      for (int i=0;i<shape[0];i++) {
        for (int j=0;j<shape[1];j++) {
          tmp=5.8*(offset+getRelIndex(shape,i,j));
          assert(std::abs(data[offset+getRelIndex(shape,i,j)]-tmp)<=REL_TOL*std::abs(tmp));
        }
      }

      if (p<npoints-1) {
        offset+=noValues(shape);
      }

    }

  }

  {
    cout << endl;
    cout << "\tTest binaryOp on shape (9,8,5,11) Data and single value.";

    // define the shape for the DataArrayView
    DataTypes::ShapeType shape;
    shape.push_back(9);
    shape.push_back(8);
    shape.push_back(5);
    shape.push_back(11);

    // allocate the data for the DataArrayView
    int npoints=4;
    DataTypes::ValueType data(DataTypes::noValues(shape)*npoints,0);

    int offset=0;
    // step the view along each data point in the underlying data
    for (int p=0;p<npoints;p++) {

      // assign values to the data point
      for (int i=0;i<shape[0];i++) {
        for (int j=0;j<shape[1];j++) {
          for (int k=0;k<shape[2];k++) {
            for (int l=0;l<shape[3];l++) {
              data[offset+getRelIndex(shape,i,j,k,l)]=offset+getRelIndex(shape,i,j,k,l);
            }
          }
        }
      }

      // apply a binary operation to the data point
//       dataView.binaryOp(5.4,multiplies<double>());
      binaryOp(data,shape,offset,5.4,multiplies<double>());

      double tmp;
      // check the results
      for (int i=0;i<shape[0];i++) {
        for (int j=0;j<shape[1];j++) {
          for (int k=0;k<shape[2];k++) {
            for (int l=0;l<shape[3];l++) {
              tmp=5.4*(offset+getRelIndex(shape,i,j,k,l));
              assert(std::abs(data[offset+getRelIndex(shape,i,j,k,l)]-tmp)<=REL_TOL*std::abs(tmp));
            }
          }
        }
      }

      if (p<npoints-1) {
        offset+=noValues(shape);
      }

    }

  }

  cout << endl;

}

void DataMathsTestCase::testReductionOp()
{

  {
    cout << endl;
    cout << "\tTest reductionOp on scalar Data.";

    // define the shape for the DataArrayView
    DataTypes::ShapeType shape;

    // allocate the data for the DataArrayView
    int npoints=4;
    DataTypes::ValueType data(DataTypes::noValues(shape)*npoints,0);

    int offset=0;
    // step the view along each data point in the underlying data
    for (int p=0;p<npoints;p++) {

      // assign values to the data point
      data[offset]=p;

      // apply a reduction operation to this data point and check the results
      FMax fmax_func;
//       assert(std::abs(dataView.reductionOp(fmax_func,numeric_limits<double>::max()*-1)-p)<=REL_TOL*p);
      assert(std::abs(reductionOp(data,shape,offset,fmax_func,numeric_limits<double>::max()*-1)-p)<=REL_TOL*p);


      if (p<npoints-1) {
        ++offset;
      }

    }

  }

  {
    cout << endl;
    cout << "\tTest reductionOp on shape (2,3) Data.";

    // define the shape for the DataArrayView
    DataTypes::ShapeType shape;
    shape.push_back(2);
    shape.push_back(3);

    // allocate the data for the DataArrayView
    int npoints=4;
    DataTypes::ValueType data(DataTypes::noValues(shape)*npoints,0);

    int offset=0;
    // step the view along each data point in the underlying data
    for (int p=0;p<npoints;p++) {

      // assign values to the data point
      for (int i=0;i<shape[0];i++) {
        for (int j=0;j<shape[1];j++) {
          data[offset+getRelIndex(shape,i,j)]=offset+getRelIndex(shape,i,j);
        }
      }

      // apply a reduction operation to this data point and check the results
      FMin fmin_func;
      assert(std::abs(reductionOp(data,shape,offset,fmin_func,numeric_limits<double>::max())-offset)<=REL_TOL*std::abs(offset));

      if (p<npoints-1) {
        offset+=noValues(shape);
      }

    }

  }

  {
    cout << endl;
    cout << "\tTest reductionOp on shape (9,8,5,11) Data.";

    // define the shape for the DataArrayView
    DataTypes::ShapeType shape;
    shape.push_back(9);
    shape.push_back(8);
    shape.push_back(5);
    shape.push_back(11);

    // allocate the data for the DataArrayView
    int npoints=4;
    DataTypes::ValueType data(DataTypes::noValues(shape)*npoints,0);

    int offset=0;
    // step the view along each data point in the underlying data
    for (int p=0;p<npoints;p++) {

      // assign values to the data point
      for (int i=0;i<shape[0];i++) {
        for (int j=0;j<shape[1];j++) {
          for (int k=0;k<shape[2];k++) {
            for (int l=0;l<shape[3];l++) {
              data[offset+getRelIndex(shape,i,j,k,l)]=offset+getRelIndex(shape,i,j,k,l);
            }
          }
        }
      }

      // apply a reduction operation to this data point and check the results
      AbsMax absmax_func;
      assert(reductionOp(data,shape,offset,absmax_func,0)==offset+getRelIndex(shape,8,7,4,10));

      if (p<npoints-1) {
        offset+=noValues(shape);
      }

    }

  }

  cout << endl;

}

TestSuite* DataMathsTestCase::suite ()
{
  //
  // create the suite of tests to perform.
  TestSuite *testSuite = new TestSuite ("DataMathsTestCase");

  testSuite->addTest (new TestCaller< DataMathsTestCase>("testUnaryOp",&DataMathsTestCase::testUnaryOp));
  testSuite->addTest (new TestCaller< DataMathsTestCase>("testBinaryOp",&DataMathsTestCase::testBinaryOp));
  testSuite->addTest (new TestCaller< DataMathsTestCase>("testReductionOp",&DataMathsTestCase::testReductionOp));
  testSuite->addTest (new TestCaller< DataMathsTestCase>("testMatMult",&DataMathsTestCase::testMatMult));
  return testSuite;
}
1
2	/*******************************************************
3	*
4	* Copyright (c) 2003-2008 by University of Queensland
5	* Earth Systems Science Computational Center (ESSCC)
6	* http://www.uq.edu.au/esscc
7	*
8	* Primary Business: Queensland, Australia
9	* Licensed under the Open Software License version 3.0
10	* http://www.opensource.org/licenses/osl-3.0.php
11	*
12	*******************************************************/
13
14
15	#include "escript/DataAlgorithm.h"
16	#include "escript/DataVector.h"
17	#include "esysUtils/EsysException.h"
18
19	#include "DataMathsTestCase.h"
20	#include "escript/DataTypes.h"
21
22	#include <iostream>
23
24	using namespace CppUnitTest;
25	using namespace esysUtils;
26	using namespace escript;
27	using namespace std;
28	using namespace escript::DataTypes;
29	using namespace escript::DataMaths;
30
31	void DataMathsTestCase::setUp() {
32	//
33	// This is called before each test is run
34
35	}
36
37	void DataMathsTestCase::tearDown() {
38	//
39	// This is called after each test has been run
40
41	}
42
43
44	void DataMathsTestCase::testMatMult()
45	{
46
47	{
48	cout << endl;
49	cout << "\tTest result shape." << endl;
50
51	DataTypes::ShapeType leftShape;
52	leftShape.push_back(1);
53	leftShape.push_back(3);
54	DataTypes::ValueType leftData(DataTypes::noValues(leftShape),0);
55	// DataArrayView leftDataView(leftData,leftShape);
56
57	DataTypes::ShapeType rightShape;
58	rightShape.push_back(3);
59	rightShape.push_back(2);
60	DataTypes::ValueType rightData(DataTypes::noValues(rightShape),0);
61	// DataArrayView rightDataView(rightData,rightShape);
62
63	DataTypes::ShapeType resultShape=DataMaths::determineResultShape(leftShape,rightShape);
64
65	assert(resultShape.size()==2);
66	assert(resultShape[0]==1);
67	assert(resultShape[1]==2);
68
69	DataTypes::ValueType resultData(DataTypes::noValues(resultShape),0);
70
71	cout << "\tTest matrix multiplication.";
72	double aValue=0.0;
73	for (int i=0;i<leftShape[0];i++) {
74	for (int j=0;j<leftShape[1];j++) {
75	leftData[getRelIndex(leftShape,i,j)]=++aValue;
76	}
77	}
78	aValue=0.0;
79	for (int i=0;i<rightShape[0];i++) {
80	for (int j=0;j<rightShape[1];j++) {
81	rightData[getRelIndex(rightShape,i,j)]=++aValue;
82	}
83	}
84
85	DataMaths::matMult(leftData,leftShape,0,rightData,rightShape,0,resultData, resultShape);
86	assert((resultData[0]==22) && (resultData[1]==28));
87	}
88
89	cout << endl;
90
91	}
92
93	void DataMathsTestCase::testUnaryOp()
94	{
95
96	// This typedef allows function names to be cast to pointers
97	// to unary functions of the appropriate type.
98	typedef double (*UnaryDFunPtr)(double);
99
100	{
101	cout << endl;
102	cout << "\tTest unaryOp on scalar Data.";
103
104	// define the shape for the Data
105	DataTypes::ShapeType shape;
106
107	// allocate the data for the Data
108	int npoints=4;
109	DataTypes::ValueType data(DataTypes::noValues(shape)*npoints,0);
110
111	double tmp;
112	int offset=0;
113	// step the view along each data point in the underlying data
114	for (int p=0;p<npoints;p++) {
115
116	// assign values to the data point
117	data[offset]=p;
118
119	// apply a unary operation to this data point
120	unaryOp(data,shape,offset,(UnaryDFunPtr)std::sin);
121
122	// check the results
123	tmp = std::sin((double)p);
124	assert(std::abs(data[offset]-tmp)<=REL_TOL*std::abs(tmp));
125
126	if (p<npoints-1) {
127	offset++;
128	}
129
130	}
131
132	}
133
134	{
135	cout << endl;
136	cout << "\tTest unaryOp on shape (2,3) Data.";
137
138	// define the shape for the DataArrayView
139	DataTypes::ShapeType shape;
140	shape.push_back(2);
141	shape.push_back(3);
142
143	// allocate the data for the DataArrayView
144	int npoints=4;
145	DataTypes::ValueType data(DataTypes::noValues(shape)*npoints,0);
146
147
148	int offset=0;
149	// step the view along each data point in the underlying data
150	for (int p=0;p<npoints;p++) {
151
152	// assign values to the data point
153	for (int i=0;i<shape[0];i++) {
154	for (int j=0;j<shape[1];j++) {
155	data[offset+getRelIndex(shape,i,j)]=offset+getRelIndex(shape,i,j);
156	}
157	}
158
159	// apply a unary operation to this data point
160	// dataView.unaryOp((UnaryDFunPtr)std::sqrt);
161	unaryOp(data,shape,offset,(UnaryDFunPtr)std::sqrt);
162
163	// check the results
164	for (int i=0;i<shape[0];i++) {
165	for (int j=0;j<shape[1];j++) {
166	// assert(std::abs(dataView(i,j)-std::sqrt((double)dataView.index(i,j)))<=REL_TOL*std::sqrt((double)dataView.index(i,j)));
167	assert(std::abs(data[offset+getRelIndex(shape,i,j)]-std::sqrt((double)offset+getRelIndex(shape,i,j)))<=REL_TOL*std::sqrt((double)offset+getRelIndex(shape,i,j)));
168	}
169	}
170
171	if (p<npoints-1) {
172	offset++;
173	}
174
175	}
176
177	}
178
179	{
180	cout << endl;
181	cout << "\tTest unaryOp on shape (9,8,5,11) Data.";
182
183	// define the shape for the DataArrayView
184	DataTypes::ShapeType shape;
185	shape.push_back(9);
186	shape.push_back(8);
187	shape.push_back(5);
188	shape.push_back(11);
189
190	// allocate the data for the DataArrayView
191	int npoints=4;
192	DataTypes::ValueType data(DataTypes::noValues(shape)*npoints,0);
193
194	int offset=0;
195	// step the view along each data point in the underlying data
196	for (int p=0;p<npoints;p++) {
197
198	// assign values to the data point
199	for (int i=0;i<shape[0];i++) {
200	for (int j=0;j<shape[1];j++) {
201	for (int k=0;k<shape[2];k++) {
202	for (int l=0;l<shape[3];l++) {
203	data[offset+getRelIndex(shape,i,j,k,l)]=data[offset+getRelIndex(shape,i,j,k,l)]+1;
204	}
205	}
206	}
207	}
208
209	// apply a unary operation to this data point
210	// dataView.unaryOp((UnaryDFunPtr)std::log);
211	unaryOp(data,shape,offset,(UnaryDFunPtr)std::log);
212
213	// check the results
214	for (int i=0;i<shape[0];i++) {
215	for (int j=0;j<shape[1];j++) {
216	for (int k=0;k<shape[2];k++) {
217	for (int l=0;l<shape[3];l++) {
218	assert(std::abs(data[offset+getRelIndex(shape,i,j,k,l)]-std::log(1+(double)data[offset+getRelIndex(shape,i,j,k,l)]))<=REL_TOL*std::abs(std::log(1+(double)data[offset+getRelIndex(shape,i,j,k,l)])));
219	}
220	}
221	}
222	}
223
224	if (p<npoints-1) {
225	offset++;
226	}
227
228	}
229
230	}
231
232	cout << endl;
233
234	}
235
236	void DataMathsTestCase::testBinaryOp()
237	{
238
239	// This typedef allows function names to be cast to pointers
240	// to binary functions of the appropriate type.
241	typedef double (*BinaryDFunPtr)(double,double);
242
243	{
244	cout << endl;
245	cout << "\tTest binaryOp on scalar Data.";
246
247	// define the shape for the DataArrayViews
248	DataTypes::ShapeType shape;
249
250	// allocate the data for the DataArrayViews
251	int npoints=4;
252	DataTypes::ValueType data1(DataTypes::noValues(shape)*npoints,0);
253	DataTypes::ValueType data2(DataTypes::noValues(shape)*npoints,0);
254
255	int offset=0;
256	// step the views along each data point in the underlying data
257	for (int p=0;p<npoints;p++) {
258
259	// assign values to the data points
260	data1[offset]=p;
261	data2[offset]=p;
262
263	// apply a binary operation to these data points
264	binaryOp(data1,scalarShape,offset,data2,scalarShape,offset, plus<double>());
265
266	// check the results
267	assert(data1[offset]==p+p);
268
269	if (p<npoints-1) {
270	++offset;
271	}
272
273	}
274
275	}
276
277	{
278	cout << endl;
279	cout << "\tTest binaryOp on shape (2,3) Data.";
280
281	// define the shape for the DataArrayViews
282	DataTypes::ShapeType shape;
283	shape.push_back(2);
284	shape.push_back(3);
285
286	// allocate the data for the DataArrayViews
287	int npoints=4;
288	DataTypes::ValueType data1(DataTypes::noValues(shape)*npoints,0);
289	DataTypes::ValueType data2(DataTypes::noValues(shape)*npoints,0);
290
291	int offset=0;
292	// step the views along each data point in the underlying data
293	for (int p=0;p<npoints;p++) {
294
295	// assign values to the data points
296	for (int i=0;i<shape[0];i++) {
297	for (int j=0;j<shape[1];j++) {
298	data1[offset+getRelIndex(shape,i,j)]=offset+getRelIndex(shape,i,j);
299	data2[offset+getRelIndex(shape,i,j)]=offset+getRelIndex(shape,i,j);
300	}
301	}
302
303	// apply a binary operation to these data points
304	/* dataView1.binaryOp(dataView2,multiplies<double>());*/
305	binaryOp(data1,shape,offset,data2,shape,offset,multiplies<double>());
306
307	// check the results
308	for (int i=0;i<shape[0];i++) {
309	for (int j=0;j<shape[1];j++) {
310	assert(data1[offset+getRelIndex(shape,i,j)]==(offset+getRelIndex(shape,i,j))*(offset+getRelIndex(shape,i,j)));
311	}
312	}
313
314	if (p<npoints-1) {
315	offset+=noValues(shape);
316	}
317
318	}
319
320	}
321
322	{
323	cout << endl;
324	cout << "\tTest binaryOp on shape (9,8,5,11) Data.";
325
326	// define the shape for the DataArrayViews
327	DataTypes::ShapeType shape;
328	shape.push_back(9);
329	shape.push_back(8);
330	shape.push_back(5);
331	shape.push_back(11);
332
333	// allocate the data for the DataArrayViews
334	int npoints=4;
335	DataTypes::ValueType data1(DataTypes::noValues(shape)*npoints,0);
336	DataTypes::ValueType data2(DataTypes::noValues(shape)*npoints,0);
337
338	int offset=0;
339	// step the views along each data point in the underlying data
340	for (int p=0;p<npoints;p++) {
341
342	// assign values to the data points
343	for (int i=0;i<shape[0];i++) {
344	for (int j=0;j<shape[1];j++) {
345	for (int k=0;k<shape[2];k++) {
346	for (int l=0;l<shape[3];l++) {
347	data1[offset+getRelIndex(shape,i,j,k,l)]=offset+getRelIndex(shape,i,j,k,l);
348	data2[offset+getRelIndex(shape,i,j,k,l)]=offset+getRelIndex(shape,i,j,k,l);
349	}
350	}
351	}
352	}
353
354	// apply a binary operation to these data points
355	// dataView1.binaryOp(dataView2,multiplies<double>());
356	binaryOp(data1,shape,offset,data2,shape,offset,multiplies<double>());
357
358	// check the results
359	for (int i=0;i<shape[0];i++) {
360	for (int j=0;j<shape[1];j++) {
361	for (int k=0;k<shape[2];k++) {
362	for (int l=0;l<shape[3];l++) {
363	assert(data1[offset+getRelIndex(shape,i,j,k,l)]==(offset+getRelIndex(shape,i,j,k,l))*(offset+getRelIndex(shape,i,j,k,l)));
364	}
365	}
366	}
367	}
368
369	if (p<npoints-1) {
370	offset+=noValues(shape);
371	}
372
373	}
374
375	}
376
377	{
378	cout << endl;
379	cout << "\tTest binaryOp on scalar Data and single value.";
380
381	// define the shape for the DataArrayView
382	DataTypes::ShapeType shape;
383
384	// allocate the data for the DataArrayView
385	int npoints=4;
386	DataTypes::ValueType data(DataTypes::noValues(shape)*npoints,0);
387
388	int offset=0;
389	// step the view along each data point in the underlying data
390	for (int p=0;p<npoints;p++) {
391
392	// assign values to the data point
393	data[offset]=p;
394
395	// apply a binary operation to this data point
396	// dataView.binaryOp(4.9,plus<double>());
397	binaryOp(data,shape,offset,4.9,plus<double>());
398
399	// check the results
400	assert(data[offset]==4.9+p);
401
402	if (p<npoints-1) {
403	++offset;
404	}
405
406	}
407
408	}
409
410	{
411	cout << endl;
412	cout << "\tTest binaryOp on shape (2,3) Data and single value.";
413
414	// define the shape for the DataArrayView
415	DataTypes::ShapeType shape;
416	shape.push_back(2);
417	shape.push_back(3);
418
419	// allocate the data for the DataArrayView
420	int npoints=4;
421	DataTypes::ValueType data(DataTypes::noValues(shape)*npoints,0);
422
423	int offset=0;
424	// step the view along each data point in the underlying data
425	for (int p=0;p<npoints;p++) {
426
427	// assign values to the data point
428	for (int i=0;i<shape[0];i++) {
429	for (int j=0;j<shape[1];j++) {
430	data[offset+getRelIndex(shape,i,j)]=offset+getRelIndex(shape,i,j);
431	}
432	}
433
434	// apply a binary operation to the data point
435	// dataView.binaryOp(5.8,multiplies<double>());
436	binaryOp(data,shape,offset,5.8,multiplies<double>());
437
438	double tmp;
439	// check the results
440	for (int i=0;i<shape[0];i++) {
441	for (int j=0;j<shape[1];j++) {
442	tmp=5.8*(offset+getRelIndex(shape,i,j));
443	assert(std::abs(data[offset+getRelIndex(shape,i,j)]-tmp)<=REL_TOL*std::abs(tmp));
444	}
445	}
446
447	if (p<npoints-1) {
448	offset+=noValues(shape);
449	}
450
451	}
452
453	}
454
455	{
456	cout << endl;
457	cout << "\tTest binaryOp on shape (9,8,5,11) Data and single value.";
458
459	// define the shape for the DataArrayView
460	DataTypes::ShapeType shape;
461	shape.push_back(9);
462	shape.push_back(8);
463	shape.push_back(5);
464	shape.push_back(11);
465
466	// allocate the data for the DataArrayView
467	int npoints=4;
468	DataTypes::ValueType data(DataTypes::noValues(shape)*npoints,0);
469
470	int offset=0;
471	// step the view along each data point in the underlying data
472	for (int p=0;p<npoints;p++) {
473
474	// assign values to the data point
475	for (int i=0;i<shape[0];i++) {
476	for (int j=0;j<shape[1];j++) {
477	for (int k=0;k<shape[2];k++) {
478	for (int l=0;l<shape[3];l++) {
479	data[offset+getRelIndex(shape,i,j,k,l)]=offset+getRelIndex(shape,i,j,k,l);
480	}
481	}
482	}
483	}
484
485	// apply a binary operation to the data point
486	// dataView.binaryOp(5.4,multiplies<double>());
487	binaryOp(data,shape,offset,5.4,multiplies<double>());
488
489	double tmp;
490	// check the results
491	for (int i=0;i<shape[0];i++) {
492	for (int j=0;j<shape[1];j++) {
493	for (int k=0;k<shape[2];k++) {
494	for (int l=0;l<shape[3];l++) {
495	tmp=5.4*(offset+getRelIndex(shape,i,j,k,l));
496	assert(std::abs(data[offset+getRelIndex(shape,i,j,k,l)]-tmp)<=REL_TOL*std::abs(tmp));
497	}
498	}
499	}
500	}
501
502	if (p<npoints-1) {
503	offset+=noValues(shape);
504	}
505
506	}
507
508	}
509
510	cout << endl;
511
512	}
513
514	void DataMathsTestCase::testReductionOp()
515	{
516
517	{
518	cout << endl;
519	cout << "\tTest reductionOp on scalar Data.";
520
521	// define the shape for the DataArrayView
522	DataTypes::ShapeType shape;
523
524	// allocate the data for the DataArrayView
525	int npoints=4;
526	DataTypes::ValueType data(DataTypes::noValues(shape)*npoints,0);
527
528	int offset=0;
529	// step the view along each data point in the underlying data
530	for (int p=0;p<npoints;p++) {
531
532	// assign values to the data point
533	data[offset]=p;
534
535	// apply a reduction operation to this data point and check the results
536	FMax fmax_func;
537	// assert(std::abs(dataView.reductionOp(fmax_func,numeric_limits<double>::max()-1)-p)<=REL_TOLp);
538	assert(std::abs(reductionOp(data,shape,offset,fmax_func,numeric_limits<double>::max()-1)-p)<=REL_TOLp);
539
540
541
542	if (p<npoints-1) {
543	++offset;
544	}
545
546	}
547
548	}
549
550	{
551	cout << endl;
552	cout << "\tTest reductionOp on shape (2,3) Data.";
553
554	// define the shape for the DataArrayView
555	DataTypes::ShapeType shape;
556	shape.push_back(2);
557	shape.push_back(3);
558
559	// allocate the data for the DataArrayView
560	int npoints=4;
561	DataTypes::ValueType data(DataTypes::noValues(shape)*npoints,0);
562
563	int offset=0;
564	// step the view along each data point in the underlying data
565	for (int p=0;p<npoints;p++) {
566
567	// assign values to the data point
568	for (int i=0;i<shape[0];i++) {
569	for (int j=0;j<shape[1];j++) {
570	data[offset+getRelIndex(shape,i,j)]=offset+getRelIndex(shape,i,j);
571	}
572	}
573
574	// apply a reduction operation to this data point and check the results
575	FMin fmin_func;
576	assert(std::abs(reductionOp(data,shape,offset,fmin_func,numeric_limits<double>::max())-offset)<=REL_TOL*std::abs(offset));
577
578	if (p<npoints-1) {
579	offset+=noValues(shape);
580	}
581
582	}
583
584	}
585
586	{
587	cout << endl;
588	cout << "\tTest reductionOp on shape (9,8,5,11) Data.";
589
590	// define the shape for the DataArrayView
591	DataTypes::ShapeType shape;
592	shape.push_back(9);
593	shape.push_back(8);
594	shape.push_back(5);
595	shape.push_back(11);
596
597	// allocate the data for the DataArrayView
598	int npoints=4;
599	DataTypes::ValueType data(DataTypes::noValues(shape)*npoints,0);
600
601	int offset=0;
602	// step the view along each data point in the underlying data
603	for (int p=0;p<npoints;p++) {
604
605	// assign values to the data point
606	for (int i=0;i<shape[0];i++) {
607	for (int j=0;j<shape[1];j++) {
608	for (int k=0;k<shape[2];k++) {
609	for (int l=0;l<shape[3];l++) {
610	data[offset+getRelIndex(shape,i,j,k,l)]=offset+getRelIndex(shape,i,j,k,l);
611	}
612	}
613	}
614	}
615
616	// apply a reduction operation to this data point and check the results
617	AbsMax absmax_func;
618	assert(reductionOp(data,shape,offset,absmax_func,0)==offset+getRelIndex(shape,8,7,4,10));
619
620	if (p<npoints-1) {
621	offset+=noValues(shape);
622	}
623
624	}
625
626	}
627
628	cout << endl;
629
630	}
631
632	TestSuite* DataMathsTestCase::suite ()
633	{
634	//
635	// create the suite of tests to perform.
636	TestSuite *testSuite = new TestSuite ("DataMathsTestCase");
637
638	testSuite->addTest (new TestCaller< DataMathsTestCase>("testUnaryOp",&DataMathsTestCase::testUnaryOp));
639	testSuite->addTest (new TestCaller< DataMathsTestCase>("testBinaryOp",&DataMathsTestCase::testBinaryOp));
640	testSuite->addTest (new TestCaller< DataMathsTestCase>("testReductionOp",&DataMathsTestCase::testReductionOp));
641	testSuite->addTest (new TestCaller< DataMathsTestCase>("testMatMult",&DataMathsTestCase::testMatMult));
642	return testSuite;
643	}