/[escript]/trunk/escript/src/Data.cpp
ViewVC logotype

Diff of /trunk/escript/src/Data.cpp

Parent Directory Parent Directory | Revision Log Revision Log | View Patch Patch

trunk/esys2/escript/src/Data/Data.cpp revision 97 by jgs, Tue Dec 14 05:39:33 2004 UTC trunk/escript/src/Data.cpp revision 1092 by gross, Fri Apr 13 03:39:49 2007 UTC
# Line 1  Line 1 
1  // $Id$  // $Id$
 /*=============================================================================  
2    
3   ******************************************************************************  /*
4   *                                                                            *   ************************************************************
5   *       COPYRIGHT ACcESS 2004 -  All Rights Reserved                         *   *          Copyright 2006 by ACcESS MNRF                   *
6   *                                                                            *   *                                                          *
7   * This software is the property of ACcESS.  No part of this code             *   *              http://www.access.edu.au                    *
8   * may be copied in any form or by any means without the expressed written    *   *       Primary Business: Queensland, Australia            *
9   * consent of ACcESS.  Copying, use or modification of this software          *   *  Licensed under the Open Software License version 3.0    *
10   * by any unauthorised person is illegal unless that                          *   *     http://www.opensource.org/licenses/osl-3.0.php       *
11   * person has a software license agreement with ACcESS.                       *   *                                                          *
12   *                                                                            *   ************************************************************
13   ******************************************************************************  */
14    #include "Data.h"
 ******************************************************************************/  
15    
16  #include "escript/Data/Data.h"  #include "DataExpanded.h"
17    #include "DataConstant.h"
18    #include "DataTagged.h"
19    #include "DataEmpty.h"
20    #include "DataArray.h"
21    #include "DataArrayView.h"
22    #include "DataProf.h"
23    #include "FunctionSpaceFactory.h"
24    #include "AbstractContinuousDomain.h"
25    #include "UnaryFuncs.h"
26    
27  #include <iostream>  #include <fstream>
28  #include <algorithm>  #include <algorithm>
29  #include <vector>  #include <vector>
 #include <exception>  
30  #include <functional>  #include <functional>
 #include <math.h>  
31    
32  #include <boost/python/str.hpp>  #include <boost/python/dict.hpp>
33  #include <boost/python/extract.hpp>  #include <boost/python/extract.hpp>
34  #include <boost/python/long.hpp>  #include <boost/python/long.hpp>
35    
 #include "escript/Data/DataException.h"  
   
 #include "escript/Data/DataExpanded.h"  
 #include "escript/Data/DataConstant.h"  
 #include "escript/Data/DataTagged.h"  
 #include "escript/Data/DataEmpty.h"  
 #include "escript/Data/DataArray.h"  
 #include "escript/Data/DataAlgorithm.h"  
 #include "escript/Data/FunctionSpaceFactory.h"  
 #include "escript/Data/AbstractContinuousDomain.h"  
 #include "escript/Data/UnaryFuncs.h"  
   
36  using namespace std;  using namespace std;
37  using namespace boost::python;  using namespace boost::python;
38  using namespace boost;  using namespace boost;
# Line 52  Data::Data() Line 45  Data::Data()
45    DataAbstract* temp=new DataEmpty();    DataAbstract* temp=new DataEmpty();
46    shared_ptr<DataAbstract> temp_data(temp);    shared_ptr<DataAbstract> temp_data(temp);
47    m_data=temp_data;    m_data=temp_data;
48      m_protected=false;
49  }  }
50    
51  Data::Data(double value,  Data::Data(double value,
# Line 65  Data::Data(double value, Line 59  Data::Data(double value,
59    }    }
60    DataArray temp(dataPointShape,value);    DataArray temp(dataPointShape,value);
61    initialise(temp.getView(),what,expanded);    initialise(temp.getView(),what,expanded);
62      m_protected=false;
63  }  }
64    
65  Data::Data(double value,  Data::Data(double value,
# Line 75  Data::Data(double value, Line 70  Data::Data(double value,
70    DataArray temp(dataPointShape,value);    DataArray temp(dataPointShape,value);
71    pair<int,int> dataShape=what.getDataShape();    pair<int,int> dataShape=what.getDataShape();
72    initialise(temp.getView(),what,expanded);    initialise(temp.getView(),what,expanded);
73      m_protected=false;
74  }  }
75    
76  Data::Data(const Data& inData)  Data::Data(const Data& inData)
77  {  {
78    m_data=inData.m_data;    m_data=inData.m_data;
79      m_protected=inData.isProtected();
80  }  }
81    
82  Data::Data(const Data& inData,  Data::Data(const Data& inData,
# Line 90  Data::Data(const Data& inData, Line 87  Data::Data(const Data& inData,
87    DataAbstract* tmp = inData.m_data->getSlice(region);    DataAbstract* tmp = inData.m_data->getSlice(region);
88    shared_ptr<DataAbstract> temp_data(tmp);    shared_ptr<DataAbstract> temp_data(tmp);
89    m_data=temp_data;    m_data=temp_data;
90      m_protected=false;
91  }  }
92    
93  Data::Data(const Data& inData,  Data::Data(const Data& inData,
# Line 99  Data::Data(const Data& inData, Line 97  Data::Data(const Data& inData,
97      m_data=inData.m_data;      m_data=inData.m_data;
98    } else {    } else {
99      Data tmp(0,inData.getPointDataView().getShape(),functionspace,true);      Data tmp(0,inData.getPointDataView().getShape(),functionspace,true);
100      // Note for Lutz, Must use a reference or pointer to a derived object      // Note: Must use a reference or pointer to a derived object
101      // in order to get polymorphic behaviour. Shouldn't really      // in order to get polymorphic behaviour. Shouldn't really
102      // be able to create an instance of AbstractDomain but that was done      // be able to create an instance of AbstractDomain but that was done
103      // as a boost python work around which may no longer be required.      // as a boost:python work around which may no longer be required.
104      const AbstractDomain& inDataDomain=inData.getDomain();      const AbstractDomain& inDataDomain=inData.getDomain();
105      if  (inDataDomain==functionspace.getDomain()) {      if  (inDataDomain==functionspace.getDomain()) {
106        inDataDomain.interpolateOnDomain(tmp,inData);        inDataDomain.interpolateOnDomain(tmp,inData);
# Line 111  Data::Data(const Data& inData, Line 109  Data::Data(const Data& inData,
109      }      }
110      m_data=tmp.m_data;      m_data=tmp.m_data;
111    }    }
112      m_protected=false;
113  }  }
114    
115  Data::Data(const DataTagged::TagListType& tagKeys,  Data::Data(const DataTagged::TagListType& tagKeys,
# Line 122  Data::Data(const DataTagged::TagListType Line 121  Data::Data(const DataTagged::TagListType
121    DataAbstract* temp=new DataTagged(tagKeys,values,defaultValue,what);    DataAbstract* temp=new DataTagged(tagKeys,values,defaultValue,what);
122    shared_ptr<DataAbstract> temp_data(temp);    shared_ptr<DataAbstract> temp_data(temp);
123    m_data=temp_data;    m_data=temp_data;
124      m_protected=false;
125    if (expanded) {    if (expanded) {
126      expand();      expand();
127    }    }
# Line 132  Data::Data(const numeric::array& value, Line 132  Data::Data(const numeric::array& value,
132             bool expanded)             bool expanded)
133  {  {
134    initialise(value,what,expanded);    initialise(value,what,expanded);
135      m_protected=false;
136  }  }
137    
138  Data::Data(const DataArrayView& value,  Data::Data(const DataArrayView& value,
# Line 139  Data::Data(const DataArrayView& value, Line 140  Data::Data(const DataArrayView& value,
140             bool expanded)             bool expanded)
141  {  {
142    initialise(value,what,expanded);    initialise(value,what,expanded);
143      m_protected=false;
144  }  }
145    
146  Data::Data(const object& value,  Data::Data(const object& value,
# Line 147  Data::Data(const object& value, Line 149  Data::Data(const object& value,
149  {  {
150    numeric::array asNumArray(value);    numeric::array asNumArray(value);
151    initialise(asNumArray,what,expanded);    initialise(asNumArray,what,expanded);
152      m_protected=false;
153  }  }
154    
155  Data::Data(const object& value,  Data::Data(const object& value,
# Line 167  Data::Data(const object& value, Line 170  Data::Data(const object& value,
170      // Create a DataConstant with the same sample shape as other      // Create a DataConstant with the same sample shape as other
171      initialise(temp.getView(),other.getFunctionSpace(),false);      initialise(temp.getView(),other.getFunctionSpace(),false);
172    }    }
173      m_protected=false;
174    }
175    
176    Data::~Data()
177    {
178    
179  }  }
180    
181  escriptDataC  escriptDataC
# Line 185  Data::getDataC() const Line 194  Data::getDataC() const
194    return temp;    return temp;
195  }  }
196    
197  tuple  const boost::python::tuple
198  Data::getShapeTuple() const  Data::getShapeTuple() const
199  {  {
200    const DataArrayView::ShapeType& shape=getDataPointShape();    const DataArrayView::ShapeType& shape=getDataPointShape();
# Line 303  Data::isConstant() const Line 312  Data::isConstant() const
312  }  }
313    
314  void  void
315    Data::setProtection()
316    {
317       m_protected=true;
318    }
319    
320    bool
321    Data::isProtected() const
322    {
323       return m_protected;
324    }
325    
326    
327    
328    void
329  Data::expand()  Data::expand()
330  {  {
331    if (isConstant()) {    if (isConstant()) {
# Line 344  Data::tag() Line 367  Data::tag()
367    }    }
368  }  }
369    
370  void  Data
371  Data::reshapeDataPoint(const DataArrayView::ShapeType& shape)  Data::oneOver() const
372  {  {
373    m_data->reshapeDataPoint(shape);    return escript::unaryOp(*this,bind1st(divides<double>(),1.));
374  }  }
375    
376  Data  Data
# Line 375  Data::whereNonPositive() const Line 398  Data::whereNonPositive() const
398  }  }
399    
400  Data  Data
401  Data::whereZero() const  Data::whereZero(double tol) const
402  {  {
403    return escript::unaryOp(*this,bind2nd(equal_to<double>(),0.0));    Data dataAbs=abs();
404      return escript::unaryOp(dataAbs,bind2nd(less_equal<double>(),tol));
405  }  }
406    
407  Data  Data
408  Data::whereNonZero() const  Data::whereNonZero(double tol) const
409  {  {
410    return escript::unaryOp(*this,bind2nd(not_equal_to<double>(),0.0));    Data dataAbs=abs();
411      return escript::unaryOp(dataAbs,bind2nd(greater<double>(),tol));
412  }  }
413    
414  Data  Data
# Line 443  Data::getDataPointShape() const Line 468  Data::getDataPointShape() const
468    return getPointDataView().getShape();    return getPointDataView().getShape();
469  }  }
470    
471    
472    
473    const
474    boost::python::numeric::array
475    Data:: getValueOfDataPoint(int dataPointNo)
476    {
477      size_t length=0;
478      int i, j, k, l;
479      //
480      // determine the rank and shape of each data point
481      int dataPointRank = getDataPointRank();
482      DataArrayView::ShapeType dataPointShape = getDataPointShape();
483    
484      //
485      // create the numeric array to be returned
486      boost::python::numeric::array numArray(0.0);
487    
488      //
489      // the shape of the returned numeric array will be the same
490      // as that of the data point
491      int arrayRank = dataPointRank;
492      DataArrayView::ShapeType arrayShape = dataPointShape;
493    
494      //
495      // resize the numeric array to the shape just calculated
496      if (arrayRank==0) {
497        numArray.resize(1);
498      }
499      if (arrayRank==1) {
500        numArray.resize(arrayShape[0]);
501      }
502      if (arrayRank==2) {
503        numArray.resize(arrayShape[0],arrayShape[1]);
504      }
505      if (arrayRank==3) {
506        numArray.resize(arrayShape[0],arrayShape[1],arrayShape[2]);
507      }
508      if (arrayRank==4) {
509        numArray.resize(arrayShape[0],arrayShape[1],arrayShape[2],arrayShape[3]);
510      }
511    
512      if (getNumDataPointsPerSample()>0) {
513           int sampleNo = dataPointNo/getNumDataPointsPerSample();
514           int dataPointNoInSample = dataPointNo - sampleNo * getNumDataPointsPerSample();
515           //
516           // Check a valid sample number has been supplied
517           if ((sampleNo >= getNumSamples()) || (sampleNo < 0 )) {
518               throw DataException("Error - Data::convertToNumArray: invalid sampleNo.");
519           }
520                  
521           //
522           // Check a valid data point number has been supplied
523           if ((dataPointNoInSample >= getNumDataPointsPerSample()) || (dataPointNoInSample < 0)) {
524               throw DataException("Error - Data::convertToNumArray: invalid dataPointNoInSample.");
525           }
526           // TODO: global error handling
527           // create a view of the data if it is stored locally
528           DataArrayView dataPointView = getDataPoint(sampleNo, dataPointNoInSample);
529            
530           switch( dataPointRank ){
531                case 0 :
532                    numArray[0] = dataPointView();
533                    break;
534                case 1 :        
535                    for( i=0; i<dataPointShape[0]; i++ )
536                        numArray[i]=dataPointView(i);
537                    break;
538                case 2 :        
539                    for( i=0; i<dataPointShape[0]; i++ )
540                        for( j=0; j<dataPointShape[1]; j++)
541                            numArray[make_tuple(i,j)]=dataPointView(i,j);
542                    break;
543                case 3 :        
544                    for( i=0; i<dataPointShape[0]; i++ )
545                        for( j=0; j<dataPointShape[1]; j++ )
546                            for( k=0; k<dataPointShape[2]; k++)
547                                numArray[make_tuple(i,j,k)]=dataPointView(i,j,k);
548                    break;
549                case 4 :
550                    for( i=0; i<dataPointShape[0]; i++ )
551                        for( j=0; j<dataPointShape[1]; j++ )
552                            for( k=0; k<dataPointShape[2]; k++ )
553                                for( l=0; l<dataPointShape[3]; l++)
554                                    numArray[make_tuple(i,j,k,l)]=dataPointView(i,j,k,l);
555                    break;
556        }
557      }
558      //
559      // return the array
560      return numArray;
561    
562    }
563    void
564    Data::setValueOfDataPointToPyObject(int dataPointNo, const boost::python::object& py_object)
565    {
566        // this will throw if the value cannot be represented
567        boost::python::numeric::array num_array(py_object);
568        setValueOfDataPointToArray(dataPointNo,num_array);
569    
570    
571    }
572    
573    void
574    Data::setValueOfDataPointToArray(int dataPointNo, const boost::python::numeric::array& num_array)
575    {
576      if (isProtected()) {
577            throw DataException("Error - attempt to update protected Data object.");
578      }
579      //
580      // check rank
581      if (num_array.getrank()<getDataPointRank())
582          throw DataException("Rank of numarray does not match Data object rank");
583    
584      //
585      // check shape of num_array
586      for (int i=0; i<getDataPointRank(); i++) {
587        if (extract<int>(num_array.getshape()[i])!=getDataPointShape()[i])
588           throw DataException("Shape of numarray does not match Data object rank");
589      }
590      //
591      // make sure data is expanded:
592      if (!isExpanded()) {
593        expand();
594      }
595      if (getNumDataPointsPerSample()>0) {
596           int sampleNo = dataPointNo/getNumDataPointsPerSample();
597           int dataPointNoInSample = dataPointNo - sampleNo * getNumDataPointsPerSample();
598           m_data->copyToDataPoint(sampleNo, dataPointNoInSample,num_array);
599      } else {
600           m_data->copyToDataPoint(-1, 0,num_array);
601      }
602    }
603    
604    void
605    Data::setValueOfDataPoint(int dataPointNo, const double value)
606    {
607      if (isProtected()) {
608            throw DataException("Error - attempt to update protected Data object.");
609      }
610      //
611      // make sure data is expanded:
612      if (!isExpanded()) {
613        expand();
614      }
615      if (getNumDataPointsPerSample()>0) {
616           int sampleNo = dataPointNo/getNumDataPointsPerSample();
617           int dataPointNoInSample = dataPointNo - sampleNo * getNumDataPointsPerSample();
618           m_data->copyToDataPoint(sampleNo, dataPointNoInSample,value);
619      } else {
620           m_data->copyToDataPoint(-1, 0,value);
621      }
622    }
623    
624    const
625    boost::python::numeric::array
626    Data::getValueOfGlobalDataPoint(int procNo, int dataPointNo)
627    {
628      size_t length=0;
629      int i, j, k, l, pos;
630      //
631      // determine the rank and shape of each data point
632      int dataPointRank = getDataPointRank();
633      DataArrayView::ShapeType dataPointShape = getDataPointShape();
634    
635      //
636      // create the numeric array to be returned
637      boost::python::numeric::array numArray(0.0);
638    
639      //
640      // the shape of the returned numeric array will be the same
641      // as that of the data point
642      int arrayRank = dataPointRank;
643      DataArrayView::ShapeType arrayShape = dataPointShape;
644    
645      //
646      // resize the numeric array to the shape just calculated
647      if (arrayRank==0) {
648        numArray.resize(1);
649      }
650      if (arrayRank==1) {
651        numArray.resize(arrayShape[0]);
652      }
653      if (arrayRank==2) {
654        numArray.resize(arrayShape[0],arrayShape[1]);
655      }
656      if (arrayRank==3) {
657        numArray.resize(arrayShape[0],arrayShape[1],arrayShape[2]);
658      }
659      if (arrayRank==4) {
660        numArray.resize(arrayShape[0],arrayShape[1],arrayShape[2],arrayShape[3]);
661      }
662    
663      // added for the MPI communication
664      length=1;
665      for( i=0; i<arrayRank; i++ ) length *= arrayShape[i];
666      double *tmpData = new double[length];
667    
668      //
669      // load the values for the data point into the numeric array.
670    
671        // updated for the MPI case
672        if( get_MPIRank()==procNo ){
673                 if (getNumDataPointsPerSample()>0) {
674                    int sampleNo = dataPointNo/getNumDataPointsPerSample();
675                    int dataPointNoInSample = dataPointNo - sampleNo * getNumDataPointsPerSample();
676                    //
677                    // Check a valid sample number has been supplied
678                    if ((sampleNo >= getNumSamples()) || (sampleNo < 0 )) {
679                      throw DataException("Error - Data::convertToNumArray: invalid sampleNo.");
680                    }
681                  
682                    //
683                    // Check a valid data point number has been supplied
684                    if ((dataPointNoInSample >= getNumDataPointsPerSample()) || (dataPointNoInSample < 0)) {
685                      throw DataException("Error - Data::convertToNumArray: invalid dataPointNoInSample.");
686                    }
687                    // TODO: global error handling
688            // create a view of the data if it is stored locally
689            DataArrayView dataPointView = getDataPoint(sampleNo, dataPointNoInSample);
690            
691            // pack the data from the view into tmpData for MPI communication
692            pos=0;
693            switch( dataPointRank ){
694                case 0 :
695                    tmpData[0] = dataPointView();
696                    break;
697                case 1 :        
698                    for( i=0; i<dataPointShape[0]; i++ )
699                        tmpData[i]=dataPointView(i);
700                    break;
701                case 2 :        
702                    for( i=0; i<dataPointShape[0]; i++ )
703                        for( j=0; j<dataPointShape[1]; j++, pos++ )
704                            tmpData[pos]=dataPointView(i,j);
705                    break;
706                case 3 :        
707                    for( i=0; i<dataPointShape[0]; i++ )
708                        for( j=0; j<dataPointShape[1]; j++ )
709                            for( k=0; k<dataPointShape[2]; k++, pos++ )
710                                tmpData[pos]=dataPointView(i,j,k);
711                    break;
712                case 4 :
713                    for( i=0; i<dataPointShape[0]; i++ )
714                        for( j=0; j<dataPointShape[1]; j++ )
715                            for( k=0; k<dataPointShape[2]; k++ )
716                                for( l=0; l<dataPointShape[3]; l++, pos++ )
717                                    tmpData[pos]=dataPointView(i,j,k,l);
718                    break;
719            }
720                }
721        }
722            #ifdef PASO_MPI
723            // broadcast the data to all other processes
724        MPI_Bcast( tmpData, length, MPI_DOUBLE, procNo, get_MPIComm() );
725            #endif
726    
727        // unpack the data
728        switch( dataPointRank ){
729            case 0 :
730                numArray[0]=tmpData[0];
731                break;
732            case 1 :        
733                for( i=0; i<dataPointShape[0]; i++ )
734                    numArray[i]=tmpData[i];
735                break;
736            case 2 :        
737                for( i=0; i<dataPointShape[0]; i++ )
738                    for( j=0; j<dataPointShape[1]; j++ )
739                       numArray[make_tuple(i,j)]=tmpData[i+j*dataPointShape[0]];
740                break;
741            case 3 :        
742                for( i=0; i<dataPointShape[0]; i++ )
743                    for( j=0; j<dataPointShape[1]; j++ )
744                        for( k=0; k<dataPointShape[2]; k++ )
745                            numArray[make_tuple(i,j,k)]=tmpData[i+dataPointShape[0]*(j*+k*dataPointShape[1])];
746                break;
747            case 4 :
748                for( i=0; i<dataPointShape[0]; i++ )
749                    for( j=0; j<dataPointShape[1]; j++ )
750                        for( k=0; k<dataPointShape[2]; k++ )
751                            for( l=0; l<dataPointShape[3]; l++ )
752                                    numArray[make_tuple(i,j,k,l)]=tmpData[i+dataPointShape[0]*(j*+dataPointShape[1]*(k+l*dataPointShape[2]))];
753                break;
754        }
755    
756        delete [] tmpData;  
757      //
758      // return the loaded array
759      return numArray;
760    }
761    
762    
763    
764  boost::python::numeric::array  boost::python::numeric::array
765  Data::integrate() const  Data::integrate() const
766  {  {
# Line 450  Data::integrate() const Line 768  Data::integrate() const
768    int rank = getDataPointRank();    int rank = getDataPointRank();
769    DataArrayView::ShapeType shape = getDataPointShape();    DataArrayView::ShapeType shape = getDataPointShape();
770    
771    
772    //    //
773    // calculate the integral values    // calculate the integral values
774    vector<double> integrals(getDataPointSize());    vector<double> integrals(getDataPointSize());
# Line 460  Data::integrate() const Line 779  Data::integrate() const
779    // and load the array with the integral values    // and load the array with the integral values
780    boost::python::numeric::array bp_array(1.0);    boost::python::numeric::array bp_array(1.0);
781    if (rank==0) {    if (rank==0) {
782        bp_array.resize(1);
783      index = 0;      index = 0;
784      bp_array[0] = integrals[index];      bp_array[0] = integrals[index];
785    }    }
# Line 471  Data::integrate() const Line 791  Data::integrate() const
791      }      }
792    }    }
793    if (rank==2) {    if (rank==2) {
794      bp_array.resize(shape[0],shape[1]);         bp_array.resize(shape[0],shape[1]);
795      for (int i=0; i<shape[0]; i++) {         for (int i=0; i<shape[0]; i++) {
796        for (int j=0; j<shape[1]; j++) {           for (int j=0; j<shape[1]; j++) {
797          index = i + shape[0] * j;             index = i + shape[0] * j;
798          bp_array[i,j] = integrals[index];             bp_array[make_tuple(i,j)] = integrals[index];
799        }           }
800      }         }
801    }    }
802    if (rank==3) {    if (rank==3) {
803      bp_array.resize(shape[0],shape[1],shape[2]);      bp_array.resize(shape[0],shape[1],shape[2]);
# Line 485  Data::integrate() const Line 805  Data::integrate() const
805        for (int j=0; j<shape[1]; j++) {        for (int j=0; j<shape[1]; j++) {
806          for (int k=0; k<shape[2]; k++) {          for (int k=0; k<shape[2]; k++) {
807            index = i + shape[0] * ( j + shape[1] * k );            index = i + shape[0] * ( j + shape[1] * k );
808            bp_array[i,j,k] = integrals[index];            bp_array[make_tuple(i,j,k)] = integrals[index];
809          }          }
810        }        }
811      }      }
# Line 497  Data::integrate() const Line 817  Data::integrate() const
817          for (int k=0; k<shape[2]; k++) {          for (int k=0; k<shape[2]; k++) {
818            for (int l=0; l<shape[3]; l++) {            for (int l=0; l<shape[3]; l++) {
819              index = i + shape[0] * ( j + shape[1] * ( k + shape[2] * l ) );              index = i + shape[0] * ( j + shape[1] * ( k + shape[2] * l ) );
820              bp_array[i,j,k,l] = integrals[index];              bp_array[make_tuple(i,j,k,l)] = integrals[index];
821            }            }
822          }          }
823        }        }
# Line 528  Data::tan() const Line 848  Data::tan() const
848  }  }
849    
850  Data  Data
851  Data::log() const  Data::asin() const
852  {  {
853    return escript::unaryOp(*this,(Data::UnaryDFunPtr)::log10);    return escript::unaryOp(*this,(Data::UnaryDFunPtr)::asin);
854  }  }
855    
856  Data  Data
857  Data::ln() const  Data::acos() const
858  {  {
859    return escript::unaryOp(*this,(Data::UnaryDFunPtr)::log);    return escript::unaryOp(*this,(Data::UnaryDFunPtr)::acos);
860  }  }
861    
862  double  
863  Data::Lsup() const  Data
864    Data::atan() const
865  {  {
866    //    return escript::unaryOp(*this,(Data::UnaryDFunPtr)::atan);
   // set the initial absolute maximum value to zero  
   return algorithm(DataAlgorithmAdapter<AbsMax>(0));  
867  }  }
868    
869  double  Data
870  Data::sup() const  Data::sinh() const
871  {  {
872    //    return escript::unaryOp(*this,(Data::UnaryDFunPtr)::sinh);
   // set the initial maximum value to min possible double  
   return algorithm(DataAlgorithmAdapter<FMax>(numeric_limits<double>::min()));  
873  }  }
874    
875  double  Data
876  Data::inf() const  Data::cosh() const
877  {  {
878    //    return escript::unaryOp(*this,(Data::UnaryDFunPtr)::cosh);
   // set the initial minimum value to max possible double  
   return algorithm(DataAlgorithmAdapter<FMin>(numeric_limits<double>::max()));  
879  }  }
880    
881  Data  Data
882  Data::maxval() const  Data::tanh() const
883  {  {
884    // not implemented - will use dp_algorithm    return escript::unaryOp(*this,(Data::UnaryDFunPtr)::tanh);
   return (*this);  
885  }  }
886    
887    
888  Data  Data
889  Data::minval() const  Data::erf() const
890  {  {
891    // not implemented - will use dp_algorithm  #ifdef _WIN32
892    return (*this);    throw DataException("Error - Data:: erf function is not supported on _WIN32 platforms.");
893    #else
894      return escript::unaryOp(*this,(Data::UnaryDFunPtr)::erf);
895    #endif
896  }  }
897    
898  Data  Data
899  Data::length() const  Data::asinh() const
900  {  {
901    // not implemented - will use dp_algorithm  #ifdef _WIN32
902    return (*this);    return escript::unaryOp(*this,escript::asinh_substitute);
903    #else
904      return escript::unaryOp(*this,(Data::UnaryDFunPtr)::asinh);
905    #endif
906  }  }
907    
908  Data  Data
909  Data::trace() const  Data::acosh() const
910  {  {
911    // not implemented - will use dp_algorithm  #ifdef _WIN32
912    return (*this);    return escript::unaryOp(*this,escript::acosh_substitute);
913    #else
914      return escript::unaryOp(*this,(Data::UnaryDFunPtr)::acosh);
915    #endif
916  }  }
917    
918  Data  Data
919  Data::transpose(int axis) const  Data::atanh() const
920  {  {
921    // not implemented  #ifdef _WIN32
922    return (*this);    return escript::unaryOp(*this,escript::atanh_substitute);
923    #else
924      return escript::unaryOp(*this,(Data::UnaryDFunPtr)::atanh);
925    #endif
926    }
927    
928    Data
929    Data::log10() const
930    {
931      return escript::unaryOp(*this,(Data::UnaryDFunPtr)::log10);
932    }
933    
934    Data
935    Data::log() const
936    {
937      return escript::unaryOp(*this,(Data::UnaryDFunPtr)::log);
938  }  }
939    
940  Data  Data
# Line 619  Data::neg() const Line 958  Data::neg() const
958  Data  Data
959  Data::pos() const  Data::pos() const
960  {  {
961    return (*this);    Data result;
962      // perform a deep copy
963      result.copy(*this);
964      return result;
965  }  }
966    
967  Data  Data
# Line 634  Data::sqrt() const Line 976  Data::sqrt() const
976    return escript::unaryOp(*this,(Data::UnaryDFunPtr)::sqrt);    return escript::unaryOp(*this,(Data::UnaryDFunPtr)::sqrt);
977  }  }
978    
979    double
980    Data::Lsup() const
981    {
982      double localValue, globalValue;
983      //
984      // set the initial absolute maximum value to zero
985    
986      AbsMax abs_max_func;
987      localValue = algorithm(abs_max_func,0);
988    #ifdef PASO_MPI
989      MPI_Allreduce( &localValue, &globalValue, 1, MPI_DOUBLE, MPI_MAX, MPI_COMM_WORLD );
990      return globalValue;
991    #else
992      return localValue;
993    #endif
994    }
995    
996    double
997    Data::Linf() const
998    {
999      double localValue, globalValue;
1000      //
1001      // set the initial absolute minimum value to max double
1002      AbsMin abs_min_func;
1003      localValue = algorithm(abs_min_func,numeric_limits<double>::max());
1004    
1005    #ifdef PASO_MPI
1006      MPI_Allreduce( &localValue, &globalValue, 1, MPI_DOUBLE, MPI_MIN, MPI_COMM_WORLD );
1007      return globalValue;
1008    #else
1009      return localValue;
1010    #endif
1011    }
1012    
1013    double
1014    Data::sup() const
1015    {
1016      double localValue, globalValue;
1017      //
1018      // set the initial maximum value to min possible double
1019      FMax fmax_func;
1020      localValue = algorithm(fmax_func,numeric_limits<double>::max()*-1);
1021    #ifdef PASO_MPI
1022      MPI_Allreduce( &localValue, &globalValue, 1, MPI_DOUBLE, MPI_MAX, MPI_COMM_WORLD );
1023      return globalValue;
1024    #else
1025      return localValue;
1026    #endif
1027    }
1028    
1029    double
1030    Data::inf() const
1031    {
1032      double localValue, globalValue;
1033      //
1034      // set the initial minimum value to max possible double
1035      FMin fmin_func;
1036      localValue = algorithm(fmin_func,numeric_limits<double>::max());
1037    #ifdef PASO_MPI
1038      MPI_Allreduce( &localValue, &globalValue, 1, MPI_DOUBLE, MPI_MIN, MPI_COMM_WORLD );
1039      return globalValue;
1040    #else
1041      return localValue;
1042    #endif
1043    }
1044    
1045    /* TODO */
1046    /* global reduction */
1047    Data
1048    Data::maxval() const
1049    {
1050      //
1051      // set the initial maximum value to min possible double
1052      FMax fmax_func;
1053      return dp_algorithm(fmax_func,numeric_limits<double>::max()*-1);
1054    }
1055    
1056    Data
1057    Data::minval() const
1058    {
1059      //
1060      // set the initial minimum value to max possible double
1061      FMin fmin_func;
1062      return dp_algorithm(fmin_func,numeric_limits<double>::max());
1063    }
1064    
1065    Data
1066    Data::swapaxes(const int axis0, const int axis1) const
1067    {
1068         int axis0_tmp,axis1_tmp;
1069         DataArrayView::ShapeType s=getDataPointShape();
1070         DataArrayView::ShapeType ev_shape;
1071         // Here's the equivalent of python s_out=s[axis_offset:]+s[:axis_offset]
1072         // which goes thru all shape vector elements starting with axis_offset (at index=rank wrap around to 0)
1073         int rank=getDataPointRank();
1074         if (rank<2) {
1075            throw DataException("Error - Data::swapaxes argument must have at least rank 2.");
1076         }
1077         if (axis0<0 || axis0>rank-1) {
1078            throw DataException("Error - Data::swapaxes: axis0 must be between 0 and rank-1=" + rank-1);
1079         }
1080         if (axis1<0 || axis1>rank-1) {
1081             throw DataException("Error - Data::swapaxes: axis1 must be between 0 and rank-1=" + rank-1);
1082         }
1083         if (axis0 == axis1) {
1084             throw DataException("Error - Data::swapaxes: axis indices must be different.");
1085         }
1086         if (axis0 > axis1) {
1087             axis0_tmp=axis1;
1088             axis1_tmp=axis0;
1089         } else {
1090             axis0_tmp=axis0;
1091             axis1_tmp=axis1;
1092         }
1093         for (int i=0; i<rank; i++) {
1094           if (i == axis0_tmp) {
1095              ev_shape.push_back(s[axis1_tmp]);
1096           } else if (i == axis1_tmp) {
1097              ev_shape.push_back(s[axis0_tmp]);
1098           } else {
1099              ev_shape.push_back(s[i]);
1100           }
1101         }
1102         Data ev(0.,ev_shape,getFunctionSpace());
1103         ev.typeMatchRight(*this);
1104         m_data->swapaxes(ev.m_data.get(), axis0_tmp, axis1_tmp);
1105         return ev;
1106    
1107    }
1108    
1109    Data
1110    Data::symmetric() const
1111    {
1112         // check input
1113         DataArrayView::ShapeType s=getDataPointShape();
1114         if (getDataPointRank()==2) {
1115            if(s[0] != s[1])
1116               throw DataException("Error - Data::symmetric can only be calculated for rank 2 object with equal first and second dimension.");
1117         }
1118         else if (getDataPointRank()==4) {
1119            if(!(s[0] == s[2] && s[1] == s[3]))
1120               throw DataException("Error - Data::symmetric can only be calculated for rank 4 object with dim0==dim2 and dim1==dim3.");
1121         }
1122         else {
1123            throw DataException("Error - Data::symmetric can only be calculated for rank 2 or 4 object.");
1124         }
1125         Data ev(0.,getDataPointShape(),getFunctionSpace());
1126         ev.typeMatchRight(*this);
1127         m_data->symmetric(ev.m_data.get());
1128         return ev;
1129    }
1130    
1131    Data
1132    Data::nonsymmetric() const
1133    {
1134         // check input
1135         DataArrayView::ShapeType s=getDataPointShape();
1136         if (getDataPointRank()==2) {
1137            if(s[0] != s[1])
1138               throw DataException("Error - Data::nonsymmetric can only be calculated for rank 2 object with equal first and second dimension.");
1139            DataArrayView::ShapeType ev_shape;
1140            ev_shape.push_back(s[0]);
1141            ev_shape.push_back(s[1]);
1142            Data ev(0.,ev_shape,getFunctionSpace());
1143            ev.typeMatchRight(*this);
1144            m_data->nonsymmetric(ev.m_data.get());
1145            return ev;
1146         }
1147         else if (getDataPointRank()==4) {
1148            if(!(s[0] == s[2] && s[1] == s[3]))
1149               throw DataException("Error - Data::nonsymmetric can only be calculated for rank 4 object with dim0==dim2 and dim1==dim3.");
1150            DataArrayView::ShapeType ev_shape;
1151            ev_shape.push_back(s[0]);
1152            ev_shape.push_back(s[1]);
1153            ev_shape.push_back(s[2]);
1154            ev_shape.push_back(s[3]);
1155            Data ev(0.,ev_shape,getFunctionSpace());
1156            ev.typeMatchRight(*this);
1157            m_data->nonsymmetric(ev.m_data.get());
1158            return ev;
1159         }
1160         else {
1161            throw DataException("Error - Data::nonsymmetric can only be calculated for rank 2 or 4 object.");
1162         }
1163    }
1164    
1165    Data
1166    Data::trace(int axis_offset) const
1167    {
1168         DataArrayView::ShapeType s=getDataPointShape();
1169         if (getDataPointRank()==2) {
1170            DataArrayView::ShapeType ev_shape;
1171            Data ev(0.,ev_shape,getFunctionSpace());
1172            ev.typeMatchRight(*this);
1173            m_data->trace(ev.m_data.get(), axis_offset);
1174            return ev;
1175         }
1176         if (getDataPointRank()==3) {
1177            DataArrayView::ShapeType ev_shape;
1178            if (axis_offset==0) {
1179              int s2=s[2];
1180              ev_shape.push_back(s2);
1181            }
1182            else if (axis_offset==1) {
1183              int s0=s[0];
1184              ev_shape.push_back(s0);
1185            }
1186            Data ev(0.,ev_shape,getFunctionSpace());
1187            ev.typeMatchRight(*this);
1188            m_data->trace(ev.m_data.get(), axis_offset);
1189            return ev;
1190         }
1191         if (getDataPointRank()==4) {
1192            DataArrayView::ShapeType ev_shape;
1193            if (axis_offset==0) {
1194              ev_shape.push_back(s[2]);
1195              ev_shape.push_back(s[3]);
1196            }
1197            else if (axis_offset==1) {
1198              ev_shape.push_back(s[0]);
1199              ev_shape.push_back(s[3]);
1200            }
1201        else if (axis_offset==2) {
1202          ev_shape.push_back(s[0]);
1203          ev_shape.push_back(s[1]);
1204        }
1205            Data ev(0.,ev_shape,getFunctionSpace());
1206            ev.typeMatchRight(*this);
1207        m_data->trace(ev.m_data.get(), axis_offset);
1208            return ev;
1209         }
1210         else {
1211            throw DataException("Error - Data::trace can only be calculated for rank 2, 3 or 4 object.");
1212         }
1213    }
1214    
1215    Data
1216    Data::transpose(int axis_offset) const
1217    {
1218         DataArrayView::ShapeType s=getDataPointShape();
1219         DataArrayView::ShapeType ev_shape;
1220         // Here's the equivalent of python s_out=s[axis_offset:]+s[:axis_offset]
1221         // which goes thru all shape vector elements starting with axis_offset (at index=rank wrap around to 0)
1222         int rank=getDataPointRank();
1223         if (axis_offset<0 || axis_offset>rank) {
1224            throw DataException("Error - Data::transpose must have 0 <= axis_offset <= rank=" + rank);
1225         }
1226         for (int i=0; i<rank; i++) {
1227           int index = (axis_offset+i)%rank;
1228           ev_shape.push_back(s[index]); // Append to new shape
1229         }
1230         Data ev(0.,ev_shape,getFunctionSpace());
1231         ev.typeMatchRight(*this);
1232         m_data->transpose(ev.m_data.get(), axis_offset);
1233         return ev;
1234    }
1235    
1236    Data
1237    Data::eigenvalues() const
1238    {
1239         // check input
1240         DataArrayView::ShapeType s=getDataPointShape();
1241         if (getDataPointRank()!=2)
1242            throw DataException("Error - Data::eigenvalues can only be calculated for rank 2 object.");
1243         if(s[0] != s[1])
1244            throw DataException("Error - Data::eigenvalues can only be calculated for object with equal first and second dimension.");
1245         // create return
1246         DataArrayView::ShapeType ev_shape(1,s[0]);
1247         Data ev(0.,ev_shape,getFunctionSpace());
1248         ev.typeMatchRight(*this);
1249         m_data->eigenvalues(ev.m_data.get());
1250         return ev;
1251    }
1252    
1253    const boost::python::tuple
1254    Data::eigenvalues_and_eigenvectors(const double tol) const
1255    {
1256         DataArrayView::ShapeType s=getDataPointShape();
1257         if (getDataPointRank()!=2)
1258            throw DataException("Error - Data::eigenvalues and eigenvectors can only be calculated for rank 2 object.");
1259         if(s[0] != s[1])
1260            throw DataException("Error - Data::eigenvalues and eigenvectors can only be calculated for object with equal first and second dimension.");
1261         // create return
1262         DataArrayView::ShapeType ev_shape(1,s[0]);
1263         Data ev(0.,ev_shape,getFunctionSpace());
1264         ev.typeMatchRight(*this);
1265         DataArrayView::ShapeType V_shape(2,s[0]);
1266         Data V(0.,V_shape,getFunctionSpace());
1267         V.typeMatchRight(*this);
1268         m_data->eigenvalues_and_eigenvectors(ev.m_data.get(),V.m_data.get(),tol);
1269         return make_tuple(boost::python::object(ev),boost::python::object(V));
1270    }
1271    
1272    const boost::python::tuple
1273    Data::minGlobalDataPoint() const
1274    {
1275      // NB: calc_minGlobalDataPoint( had to be split off from minGlobalDataPoint( as boost::make_tuple causes an
1276      // abort (for unknown reasons) if there are openmp directives with it in the
1277      // surrounding function
1278    
1279      int DataPointNo;
1280      int ProcNo;
1281      calc_minGlobalDataPoint(ProcNo,DataPointNo);
1282      return make_tuple(ProcNo,DataPointNo);
1283    }
1284    
1285    void
1286    Data::calc_minGlobalDataPoint(int& ProcNo,
1287                            int& DataPointNo) const
1288    {
1289      int i,j;
1290      int lowi=0,lowj=0;
1291      double min=numeric_limits<double>::max();
1292    
1293      Data temp=minval();
1294    
1295      int numSamples=temp.getNumSamples();
1296      int numDPPSample=temp.getNumDataPointsPerSample();
1297    
1298      double next,local_min;
1299      int local_lowi,local_lowj;
1300    
1301      #pragma omp parallel private(next,local_min,local_lowi,local_lowj)
1302      {
1303        local_min=min;
1304        #pragma omp for private(i,j) schedule(static)
1305        for (i=0; i<numSamples; i++) {
1306          for (j=0; j<numDPPSample; j++) {
1307            next=temp.getDataPoint(i,j)();
1308            if (next<local_min) {
1309              local_min=next;
1310              local_lowi=i;
1311              local_lowj=j;
1312            }
1313          }
1314        }
1315        #pragma omp critical
1316        if (local_min<min) {
1317          min=local_min;
1318          lowi=local_lowi;
1319          lowj=local_lowj;
1320        }
1321      }
1322    
1323    #ifdef PASO_MPI
1324        // determine the processor on which the minimum occurs
1325        next = temp.getDataPoint(lowi,lowj)();
1326        int lowProc = 0;
1327        double *globalMins = new double[get_MPISize()+1];
1328        int error = MPI_Gather ( &next, 1, MPI_DOUBLE, globalMins, 1, MPI_DOUBLE, 0, get_MPIComm() );
1329        
1330        if( get_MPIRank()==0 ){
1331            next = globalMins[lowProc];
1332            for( i=1; i<get_MPISize(); i++ )
1333                if( next>globalMins[i] ){
1334                    lowProc = i;
1335                    next = globalMins[i];
1336                }
1337        }
1338        MPI_Bcast( &lowProc, 1, MPI_DOUBLE, 0, get_MPIComm() );
1339    
1340        delete [] globalMins;
1341        ProcNo = lowProc;
1342    #else
1343        ProcNo = 0;
1344    #endif
1345      DataPointNo = lowj + lowi * numDPPSample;
1346    }
1347    
1348    void
1349    Data::saveDX(std::string fileName) const
1350    {
1351      boost::python::dict args;
1352      args["data"]=boost::python::object(this);
1353      getDomain().saveDX(fileName,args);
1354      return;
1355    }
1356    
1357    void
1358    Data::saveVTK(std::string fileName) const
1359    {
1360      boost::python::dict args;
1361      args["data"]=boost::python::object(this);
1362      getDomain().saveVTK(fileName,args);
1363      return;
1364    }
1365    
1366  Data&  Data&
1367  Data::operator+=(const Data& right)  Data::operator+=(const Data& right)
1368  {  {
1369      if (isProtected()) {
1370            throw DataException("Error - attempt to update protected Data object.");
1371      }
1372    binaryOp(right,plus<double>());    binaryOp(right,plus<double>());
1373    return (*this);    return (*this);
1374  }  }
# Line 644  Data::operator+=(const Data& right) Line 1376  Data::operator+=(const Data& right)
1376  Data&  Data&
1377  Data::operator+=(const boost::python::object& right)  Data::operator+=(const boost::python::object& right)
1378  {  {
1379    binaryOp(right,plus<double>());    Data tmp(right,getFunctionSpace(),false);
1380      binaryOp(tmp,plus<double>());
1381    return (*this);    return (*this);
1382  }  }
1383    
1384  Data&  Data&
1385  Data::operator-=(const Data& right)  Data::operator-=(const Data& right)
1386  {  {
1387      if (isProtected()) {
1388            throw DataException("Error - attempt to update protected Data object.");
1389      }
1390    binaryOp(right,minus<double>());    binaryOp(right,minus<double>());
1391    return (*this);    return (*this);
1392  }  }
# Line 658  Data::operator-=(const Data& right) Line 1394  Data::operator-=(const Data& right)
1394  Data&  Data&
1395  Data::operator-=(const boost::python::object& right)  Data::operator-=(const boost::python::object& right)
1396  {  {
1397    binaryOp(right,minus<double>());    Data tmp(right,getFunctionSpace(),false);
1398      binaryOp(tmp,minus<double>());
1399    return (*this);    return (*this);
1400  }  }
1401    
1402  Data&  Data&
1403  Data::operator*=(const Data& right)  Data::operator*=(const Data& right)
1404  {  {
1405      if (isProtected()) {
1406            throw DataException("Error - attempt to update protected Data object.");
1407      }
1408    binaryOp(right,multiplies<double>());    binaryOp(right,multiplies<double>());
1409    return (*this);    return (*this);
1410  }  }
# Line 672  Data::operator*=(const Data& right) Line 1412  Data::operator*=(const Data& right)
1412  Data&  Data&
1413  Data::operator*=(const boost::python::object& right)  Data::operator*=(const boost::python::object& right)
1414  {  {
1415    binaryOp(right,multiplies<double>());    Data tmp(right,getFunctionSpace(),false);
1416      binaryOp(tmp,multiplies<double>());
1417    return (*this);    return (*this);
1418  }  }
1419    
1420  Data&  Data&
1421  Data::operator/=(const Data& right)  Data::operator/=(const Data& right)
1422  {  {
1423      if (isProtected()) {
1424            throw DataException("Error - attempt to update protected Data object.");
1425      }
1426    binaryOp(right,divides<double>());    binaryOp(right,divides<double>());
1427    return (*this);    return (*this);
1428  }  }
# Line 686  Data::operator/=(const Data& right) Line 1430  Data::operator/=(const Data& right)
1430  Data&  Data&
1431  Data::operator/=(const boost::python::object& right)  Data::operator/=(const boost::python::object& right)
1432  {  {
1433    binaryOp(right,divides<double>());    Data tmp(right,getFunctionSpace(),false);
1434      binaryOp(tmp,divides<double>());
1435    return (*this);    return (*this);
1436  }  }
1437    
1438  Data  Data
1439    Data::rpowO(const boost::python::object& left) const
1440    {
1441      Data left_d(left,*this);
1442      return left_d.powD(*this);
1443    }
1444    
1445    Data
1446  Data::powO(const boost::python::object& right) const  Data::powO(const boost::python::object& right) const
1447  {  {
1448    Data result;    Data tmp(right,getFunctionSpace(),false);
1449    result.copy(*this);    return powD(tmp);
   result.binaryOp(right,(Data::BinaryDFunPtr)::pow);  
   return result;  
1450  }  }
1451    
1452  Data  Data
1453  Data::powD(const Data& right) const  Data::powD(const Data& right) const
1454  {  {
1455    Data result;    Data result;
1456    result.copy(*this);    if (getDataPointRank()<right.getDataPointRank()) {
1457    result.binaryOp(right,(Data::BinaryDFunPtr)::pow);       result.copy(right);
1458         result.binaryOp(*this,escript::rpow);
1459      } else {
1460         result.copy(*this);
1461         result.binaryOp(right,(Data::BinaryDFunPtr)::pow);
1462      }
1463    return result;    return result;
1464  }  }
1465    
1466    
1467  //  //
1468  // NOTE: It is essential to specify the namepsace this operator belongs to  // NOTE: It is essential to specify the namespace this operator belongs to
1469  Data  Data
1470  escript::operator+(const Data& left, const Data& right)  escript::operator+(const Data& left, const Data& right)
1471  {  {
1472    Data result;    Data result;
1473    //    //
1474    // perform a deep copy    // perform a deep copy
1475    result.copy(left);    if (left.getDataPointRank()<right.getDataPointRank()) {
1476    result+=right;       result.copy(right);
1477         result+=left;
1478      } else {
1479         result.copy(left);
1480         result+=right;
1481      }
1482    return result;    return result;
1483  }  }
1484    
1485  //  //
1486  // NOTE: It is essential to specify the namepsace this operator belongs to  // NOTE: It is essential to specify the namespace this operator belongs to
1487  Data  Data
1488  escript::operator-(const Data& left, const Data& right)  escript::operator-(const Data& left, const Data& right)
1489  {  {
1490    Data result;    Data result;
1491    //    //
1492    // perform a deep copy    // perform a deep copy
1493    result.copy(left);    if (left.getDataPointRank()<right.getDataPointRank()) {
1494    result-=right;       result=right.neg();
1495         result+=left;
1496      } else {
1497         result.copy(left);
1498         result-=right;
1499      }
1500    return result;    return result;
1501  }  }
1502    
1503  //  //
1504  // NOTE: It is essential to specify the namepsace this operator belongs to  // NOTE: It is essential to specify the namespace this operator belongs to
1505  Data  Data
1506  escript::operator*(const Data& left, const Data& right)  escript::operator*(const Data& left, const Data& right)
1507  {  {
1508    Data result;    Data result;
1509    //    //
1510    // perform a deep copy    // perform a deep copy
1511    result.copy(left);    if (left.getDataPointRank()<right.getDataPointRank()) {
1512    result*=right;       result.copy(right);
1513         result*=left;
1514      } else {
1515         result.copy(left);
1516         result*=right;
1517      }
1518    return result;    return result;
1519  }  }
1520    
1521  //  //
1522  // NOTE: It is essential to specify the namepsace this operator belongs to  // NOTE: It is essential to specify the namespace this operator belongs to
1523  Data  Data
1524  escript::operator/(const Data& left, const Data& right)  escript::operator/(const Data& left, const Data& right)
1525  {  {
1526    Data result;    Data result;
1527    //    //
1528    // perform a deep copy    // perform a deep copy
1529    result.copy(left);    if (left.getDataPointRank()<right.getDataPointRank()) {
1530    result/=right;       result=right.oneOver();
1531         result*=left;
1532      } else {
1533         result.copy(left);
1534         result/=right;
1535      }
1536    return result;    return result;
1537  }  }
1538    
1539  //  //
1540  // NOTE: It is essential to specify the namepsace this operator belongs to  // NOTE: It is essential to specify the namespace this operator belongs to
1541  Data  Data
1542  escript::operator+(const Data& left, const boost::python::object& right)  escript::operator+(const Data& left, const boost::python::object& right)
1543  {  {
1544    //    return left+Data(right,left.getFunctionSpace(),false);
   // Convert to DataArray format if possible  
   DataArray temp(right);  
   Data result;  
   //  
   // perform a deep copy  
   result.copy(left);  
   result+=right;  
   return result;  
1545  }  }
1546    
1547  //  //
1548  // NOTE: It is essential to specify the namepsace this operator belongs to  // NOTE: It is essential to specify the namespace this operator belongs to
1549  Data  Data
1550  escript::operator-(const Data& left, const boost::python::object& right)  escript::operator-(const Data& left, const boost::python::object& right)
1551  {  {
1552    //    return left-Data(right,left.getFunctionSpace(),false);
   // Convert to DataArray format if possible  
   DataArray temp(right);  
   Data result;  
   //  
   // perform a deep copy  
   result.copy(left);  
   result-=right;  
   return result;  
1553  }  }
1554    
1555  //  //
1556  // NOTE: It is essential to specify the namepsace this operator belongs to  // NOTE: It is essential to specify the namespace this operator belongs to
1557  Data  Data
1558  escript::operator*(const Data& left, const boost::python::object& right)  escript::operator*(const Data& left, const boost::python::object& right)
1559  {  {
1560    //    return left*Data(right,left.getFunctionSpace(),false);
   // Convert to DataArray format if possible  
   DataArray temp(right);  
   Data result;  
   //  
   // perform a deep copy  
   result.copy(left);  
   result*=right;  
   return result;  
1561  }  }
1562    
1563  //  //
1564  // NOTE: It is essential to specify the namepsace this operator belongs to  // NOTE: It is essential to specify the namespace this operator belongs to
1565  Data  Data
1566  escript::operator/(const Data& left, const boost::python::object& right)  escript::operator/(const Data& left, const boost::python::object& right)
1567  {  {
1568    //    return left/Data(right,left.getFunctionSpace(),false);
   // Convert to DataArray format if possible  
   DataArray temp(right);  
   Data result;  
   //  
   // perform a deep copy  
   result.copy(left);  
   result/=right;  
   return result;  
1569  }  }
1570    
1571  //  //
1572  // NOTE: It is essential to specify the namepsace this operator belongs to  // NOTE: It is essential to specify the namespace this operator belongs to
1573  Data  Data
1574  escript::operator+(const boost::python::object& left, const Data& right)  escript::operator+(const boost::python::object& left, const Data& right)
1575  {  {
1576    //    return Data(left,right.getFunctionSpace(),false)+right;
   // Construct the result using the given value and the other parameters  
   // from right  
   Data result(left,right);  
   result+=right;  
   return result;  
1577  }  }
1578    
1579  //  //
1580  // NOTE: It is essential to specify the namepsace this operator belongs to  // NOTE: It is essential to specify the namespace this operator belongs to
1581  Data  Data
1582  escript::operator-(const boost::python::object& left, const Data& right)  escript::operator-(const boost::python::object& left, const Data& right)
1583  {  {
1584    //    return Data(left,right.getFunctionSpace(),false)-right;
   // Construct the result using the given value and the other parameters  
   // from right  
   Data result(left,right);  
   result-=right;  
   return result;  
1585  }  }
1586    
1587  //  //
1588  // NOTE: It is essential to specify the namepsace this operator belongs to  // NOTE: It is essential to specify the namespace this operator belongs to
1589  Data  Data
1590  escript::operator*(const boost::python::object& left, const Data& right)  escript::operator*(const boost::python::object& left, const Data& right)
1591  {  {
1592    //    return Data(left,right.getFunctionSpace(),false)*right;
   // Construct the result using the given value and the other parameters  
   // from right  
   Data result(left,right);  
   result*=right;  
   return result;  
1593  }  }
1594    
1595  //  //
1596  // NOTE: It is essential to specify the namepsace this operator belongs to  // NOTE: It is essential to specify the namespace this operator belongs to
1597  Data  Data
1598  escript::operator/(const boost::python::object& left, const Data& right)  escript::operator/(const boost::python::object& left, const Data& right)
1599  {  {
1600    //    return Data(left,right.getFunctionSpace(),false)/right;
   // Construct the result using the given value and the other parameters  
   // from right  
   Data result(left,right);  
   result/=right;  
   return result;  
1601  }  }
1602    
1603  //  //
 // NOTE: It is essential to specify the namepsace this operator belongs to  
1604  //bool escript::operator==(const Data& left, const Data& right)  //bool escript::operator==(const Data& left, const Data& right)
1605  //{  //{
1606  //  /*  //  /*
# Line 922  escript::operator/(const boost::python:: Line 1645  escript::operator/(const boost::python::
1645  //  return ret;  //  return ret;
1646  //}  //}
1647    
1648    /* TODO */
1649    /* global reduction */
1650  Data  Data
1651  Data::getItem(const boost::python::object& key) const  Data::getItem(const boost::python::object& key) const
1652  {  {
# Line 936  Data::getItem(const boost::python::objec Line 1661  Data::getItem(const boost::python::objec
1661    return getSlice(slice_region);    return getSlice(slice_region);
1662  }  }
1663    
1664    /* TODO */
1665    /* global reduction */
1666  Data  Data
1667  Data::getSlice(const DataArrayView::RegionType& region) const  Data::getSlice(const DataArrayView::RegionType& region) const
1668  {  {
1669    return Data(*this,region);    return Data(*this,region);
1670  }  }
1671    
1672    /* TODO */
1673    /* global reduction */
1674  void  void
1675  Data::setItemO(const boost::python::object& key,  Data::setItemO(const boost::python::object& key,
1676                 const boost::python::object& value)                 const boost::python::object& value)
# Line 950  Data::setItemO(const boost::python::obje Line 1679  Data::setItemO(const boost::python::obje
1679    setItemD(key,tempData);    setItemD(key,tempData);
1680  }  }
1681    
1682    /* TODO */
1683    /* global reduction */
1684  void  void
1685  Data::setItemD(const boost::python::object& key,  Data::setItemD(const boost::python::object& key,
1686                 const Data& value)                 const Data& value)
1687  {  {
1688    const DataArrayView& view=getPointDataView();    const DataArrayView& view=getPointDataView();
1689    
1690    DataArrayView::RegionType slice_region=view.getSliceRegion(key);    DataArrayView::RegionType slice_region=view.getSliceRegion(key);
1691    if (slice_region.size()!=view.getRank()) {    if (slice_region.size()!=view.getRank()) {
1692      throw DataException("Error - slice size does not match Data rank.");      throw DataException("Error - slice size does not match Data rank.");
1693    }    }
1694    setSlice(value,slice_region);    if (getFunctionSpace()!=value.getFunctionSpace()) {
1695         setSlice(Data(value,getFunctionSpace()),slice_region);
1696      } else {
1697         setSlice(value,slice_region);
1698      }
1699  }  }
1700    
1701    /* TODO */
1702    /* global reduction */
1703  void  void
1704  Data::setSlice(const Data& value,  Data::setSlice(const Data& value,
1705                 const DataArrayView::RegionType& region)                 const DataArrayView::RegionType& region)
1706  {  {
1707      if (isProtected()) {
1708            throw DataException("Error - attempt to update protected Data object.");
1709      }
1710    Data tempValue(value);    Data tempValue(value);
1711    typeMatchLeft(tempValue);    typeMatchLeft(tempValue);
1712    typeMatchRight(tempValue);    typeMatchRight(tempValue);
# Line 1001  Data::typeMatchRight(const Data& right) Line 1742  Data::typeMatchRight(const Data& right)
1742  }  }
1743    
1744  void  void
1745    Data::setTaggedValueByName(std::string name,
1746                               const boost::python::object& value)
1747    {
1748         if (getFunctionSpace().getDomain().isValidTagName(name)) {
1749            int tagKey=getFunctionSpace().getDomain().getTag(name);
1750            setTaggedValue(tagKey,value);
1751         }
1752    }
1753    void
1754  Data::setTaggedValue(int tagKey,  Data::setTaggedValue(int tagKey,
1755                       const boost::python::object& value)                       const boost::python::object& value)
1756  {  {
1757      if (isProtected()) {
1758            throw DataException("Error - attempt to update protected Data object.");
1759      }
1760    //    //
1761    // Ensure underlying data object is of type DataTagged    // Ensure underlying data object is of type DataTagged
1762    tag();    tag();
# Line 1021  Data::setTaggedValue(int tagKey, Line 1774  Data::setTaggedValue(int tagKey,
1774    m_data->setTaggedValue(tagKey,valueDataArray.getView());    m_data->setTaggedValue(tagKey,valueDataArray.getView());
1775  }  }
1776    
 /*  
 Note: this version removed for now. Not needed, and breaks escript.cpp  
1777  void  void
1778  Data::setTaggedValue(int tagKey,  Data::setTaggedValueFromCPP(int tagKey,
1779                       const DataArrayView& value)                              const DataArrayView& value)
1780  {  {
1781      if (isProtected()) {
1782            throw DataException("Error - attempt to update protected Data object.");
1783      }
1784    //    //
1785    // Ensure underlying data object is of type DataTagged    // Ensure underlying data object is of type DataTagged
1786    tag();    tag();
# Line 1039  Data::setTaggedValue(int tagKey, Line 1793  Data::setTaggedValue(int tagKey,
1793    // Call DataAbstract::setTaggedValue    // Call DataAbstract::setTaggedValue
1794    m_data->setTaggedValue(tagKey,value);    m_data->setTaggedValue(tagKey,value);
1795  }  }
1796  */  
1797    int
1798    Data::getTagNumber(int dpno)
1799    {
1800      return m_data->getTagNumber(dpno);
1801    }
1802    
1803    void
1804    Data::archiveData(const std::string fileName)
1805    {
1806      cout << "Archiving Data object to: " << fileName << endl;
1807    
1808      //
1809      // Determine type of this Data object
1810      int dataType = -1;
1811    
1812      if (isEmpty()) {
1813        dataType = 0;
1814        cout << "\tdataType: DataEmpty" << endl;
1815      }
1816      if (isConstant()) {
1817        dataType = 1;
1818        cout << "\tdataType: DataConstant" << endl;
1819      }
1820      if (isTagged()) {
1821        dataType = 2;
1822        cout << "\tdataType: DataTagged" << endl;
1823      }
1824      if (isExpanded()) {
1825        dataType = 3;
1826        cout << "\tdataType: DataExpanded" << endl;
1827      }
1828    
1829      if (dataType == -1) {
1830        throw DataException("archiveData Error: undefined dataType");
1831      }
1832    
1833      //
1834      // Collect data items common to all Data types
1835      int noSamples = getNumSamples();
1836      int noDPPSample = getNumDataPointsPerSample();
1837      int functionSpaceType = getFunctionSpace().getTypeCode();
1838      int dataPointRank = getDataPointRank();
1839      int dataPointSize = getDataPointSize();
1840      int dataLength = getLength();
1841      DataArrayView::ShapeType dataPointShape = getDataPointShape();
1842      vector<int> referenceNumbers(noSamples);
1843      for (int sampleNo=0; sampleNo<noSamples; sampleNo++) {
1844        referenceNumbers[sampleNo] = getFunctionSpace().getReferenceIDOfSample(sampleNo);
1845      }
1846      vector<int> tagNumbers(noSamples);
1847      if (isTagged()) {
1848        for (int sampleNo=0; sampleNo<noSamples; sampleNo++) {
1849          tagNumbers[sampleNo] = getFunctionSpace().getTagFromSampleNo(sampleNo);
1850        }
1851      }
1852    
1853      cout << "\tnoSamples: " << noSamples << " noDPPSample: " << noDPPSample << endl;
1854      cout << "\tfunctionSpaceType: " << functionSpaceType << endl;
1855      cout << "\trank: " << dataPointRank << " size: " << dataPointSize << " length: " << dataLength << endl;
1856    
1857      //
1858      // Flatten Shape to an array of integers suitable for writing to file
1859      int flatShape[4] = {0,0,0,0};
1860      cout << "\tshape: < ";
1861      for (int dim=0; dim<dataPointRank; dim++) {
1862        flatShape[dim] = dataPointShape[dim];
1863        cout << dataPointShape[dim] << " ";
1864      }
1865      cout << ">" << endl;
1866    
1867      //
1868      // Open archive file
1869      ofstream archiveFile;
1870      archiveFile.open(fileName.data(), ios::out);
1871    
1872      if (!archiveFile.good()) {
1873        throw DataException("archiveData Error: problem opening archive file");
1874      }
1875    
1876      //
1877      // Write common data items to archive file
1878      archiveFile.write(reinterpret_cast<char *>(&dataType),sizeof(int));
1879      archiveFile.write(reinterpret_cast<char *>(&noSamples),sizeof(int));
1880      archiveFile.write(reinterpret_cast<char *>(&noDPPSample),sizeof(int));
1881      archiveFile.write(reinterpret_cast<char *>(&functionSpaceType),sizeof(int));
1882      archiveFile.write(reinterpret_cast<char *>(&dataPointRank),sizeof(int));
1883      archiveFile.write(reinterpret_cast<char *>(&dataPointSize),sizeof(int));
1884      archiveFile.write(reinterpret_cast<char *>(&dataLength),sizeof(int));
1885      for (int dim = 0; dim < 4; dim++) {
1886        archiveFile.write(reinterpret_cast<char *>(&flatShape[dim]),sizeof(int));
1887      }
1888      for (int sampleNo=0; sampleNo<noSamples; sampleNo++) {
1889        archiveFile.write(reinterpret_cast<char *>(&referenceNumbers[sampleNo]),sizeof(int));
1890      }
1891      if (isTagged()) {
1892        for (int sampleNo=0; sampleNo<noSamples; sampleNo++) {
1893          archiveFile.write(reinterpret_cast<char *>(&tagNumbers[sampleNo]),sizeof(int));
1894        }
1895      }
1896    
1897      if (!archiveFile.good()) {
1898        throw DataException("archiveData Error: problem writing to archive file");
1899      }
1900    
1901      //
1902      // Archive underlying data values for each Data type
1903      int noValues;
1904      switch (dataType) {
1905        case 0:
1906          // DataEmpty
1907          noValues = 0;
1908          archiveFile.write(reinterpret_cast<char *>(&noValues),sizeof(int));
1909          cout << "\tnoValues: " << noValues << endl;
1910          break;
1911        case 1:
1912          // DataConstant
1913          noValues = m_data->getLength();
1914          archiveFile.write(reinterpret_cast<char *>(&noValues),sizeof(int));
1915          cout << "\tnoValues: " << noValues << endl;
1916          if (m_data->archiveData(archiveFile,noValues)) {
1917            throw DataException("archiveData Error: problem writing data to archive file");
1918          }
1919          break;
1920        case 2:
1921          // DataTagged
1922          noValues = m_data->getLength();
1923          archiveFile.write(reinterpret_cast<char *>(&noValues),sizeof(int));
1924          cout << "\tnoValues: " << noValues << endl;
1925          if (m_data->archiveData(archiveFile,noValues)) {
1926            throw DataException("archiveData Error: problem writing data to archive file");
1927          }
1928          break;
1929        case 3:
1930          // DataExpanded
1931          noValues = m_data->getLength();
1932          archiveFile.write(reinterpret_cast<char *>(&noValues),sizeof(int));
1933          cout << "\tnoValues: " << noValues << endl;
1934          if (m_data->archiveData(archiveFile,noValues)) {
1935            throw DataException("archiveData Error: problem writing data to archive file");
1936          }
1937          break;
1938      }
1939    
1940      if (!archiveFile.good()) {
1941        throw DataException("archiveData Error: problem writing data to archive file");
1942      }
1943    
1944      //
1945      // Close archive file
1946      archiveFile.close();
1947    
1948      if (!archiveFile.good()) {
1949        throw DataException("archiveData Error: problem closing archive file");
1950      }
1951    
1952    }
1953    
1954    void
1955    Data::extractData(const std::string fileName,
1956                      const FunctionSpace& fspace)
1957    {
1958      //
1959      // Can only extract Data to an object which is initially DataEmpty
1960      if (!isEmpty()) {
1961        throw DataException("extractData Error: can only extract to DataEmpty object");
1962      }
1963    
1964      cout << "Extracting Data object from: " << fileName << endl;
1965    
1966      int dataType;
1967      int noSamples;
1968      int noDPPSample;
1969      int functionSpaceType;
1970      int dataPointRank;
1971      int dataPointSize;
1972      int dataLength;
1973      DataArrayView::ShapeType dataPointShape;
1974      int flatShape[4];
1975    
1976      //
1977      // Open the archive file
1978      ifstream archiveFile;
1979      archiveFile.open(fileName.data(), ios::in);
1980    
1981      if (!archiveFile.good()) {
1982        throw DataException("extractData Error: problem opening archive file");
1983      }
1984    
1985      //
1986      // Read common data items from archive file
1987      archiveFile.read(reinterpret_cast<char *>(&dataType),sizeof(int));
1988      archiveFile.read(reinterpret_cast<char *>(&noSamples),sizeof(int));
1989      archiveFile.read(reinterpret_cast<char *>(&noDPPSample),sizeof(int));
1990      archiveFile.read(reinterpret_cast<char *>(&functionSpaceType),sizeof(int));
1991      archiveFile.read(reinterpret_cast<char *>(&dataPointRank),sizeof(int));
1992      archiveFile.read(reinterpret_cast<char *>(&dataPointSize),sizeof(int));
1993      archiveFile.read(reinterpret_cast<char *>(&dataLength),sizeof(int));
1994      for (int dim = 0; dim < 4; dim++) {
1995        archiveFile.read(reinterpret_cast<char *>(&flatShape[dim]),sizeof(int));
1996        if (flatShape[dim]>0) {
1997          dataPointShape.push_back(flatShape[dim]);
1998        }
1999      }
2000      vector<int> referenceNumbers(noSamples);
2001      for (int sampleNo=0; sampleNo<noSamples; sampleNo++) {
2002        archiveFile.read(reinterpret_cast<char *>(&referenceNumbers[sampleNo]),sizeof(int));
2003      }
2004      vector<int> tagNumbers(noSamples);
2005      if (dataType==2) {
2006        for (int sampleNo=0; sampleNo<noSamples; sampleNo++) {
2007          archiveFile.read(reinterpret_cast<char *>(&tagNumbers[sampleNo]),sizeof(int));
2008        }
2009      }
2010    
2011      if (!archiveFile.good()) {
2012        throw DataException("extractData Error: problem reading from archive file");
2013      }
2014    
2015      //
2016      // Verify the values just read from the archive file
2017      switch (dataType) {
2018        case 0:
2019          cout << "\tdataType: DataEmpty" << endl;
2020          break;
2021        case 1:
2022          cout << "\tdataType: DataConstant" << endl;
2023          break;
2024        case 2:
2025          cout << "\tdataType: DataTagged" << endl;
2026          break;
2027        case 3:
2028          cout << "\tdataType: DataExpanded" << endl;
2029          break;
2030        default:
2031          throw DataException("extractData Error: undefined dataType read from archive file");
2032          break;
2033      }
2034    
2035      cout << "\tnoSamples: " << noSamples << " noDPPSample: " << noDPPSample << endl;
2036      cout << "\tfunctionSpaceType: " << functionSpaceType << endl;
2037      cout << "\trank: " << dataPointRank << " size: " << dataPointSize << " length: " << dataLength << endl;
2038      cout << "\tshape: < ";
2039      for (int dim = 0; dim < dataPointRank; dim++) {
2040        cout << dataPointShape[dim] << " ";
2041      }
2042      cout << ">" << endl;
2043    
2044      //
2045      // Verify that supplied FunctionSpace object is compatible with this Data object.
2046      if ( (fspace.getTypeCode()!=functionSpaceType) ||
2047           (fspace.getNumSamples()!=noSamples) ||
2048           (fspace.getNumDPPSample()!=noDPPSample)
2049         ) {
2050        throw DataException("extractData Error: incompatible FunctionSpace");
2051      }
2052      for (int sampleNo=0; sampleNo<noSamples; sampleNo++) {
2053        if (referenceNumbers[sampleNo] != fspace.getReferenceIDOfSample(sampleNo)) {
2054          throw DataException("extractData Error: incompatible FunctionSpace");
2055        }
2056      }
2057      if (dataType==2) {
2058        for (int sampleNo=0; sampleNo<noSamples; sampleNo++) {
2059          if (tagNumbers[sampleNo] != fspace.getTagFromSampleNo(sampleNo)) {
2060            throw DataException("extractData Error: incompatible FunctionSpace");
2061          }
2062        }
2063      }
2064    
2065      //
2066      // Construct a DataVector to hold underlying data values
2067      DataVector dataVec(dataLength);
2068    
2069      //
2070      // Load this DataVector with the appropriate values
2071      int noValues;
2072      archiveFile.read(reinterpret_cast<char *>(&noValues),sizeof(int));
2073      cout << "\tnoValues: " << noValues << endl;
2074      switch (dataType) {
2075        case 0:
2076          // DataEmpty
2077          if (noValues != 0) {
2078            throw DataException("extractData Error: problem reading data from archive file");
2079          }
2080          break;
2081        case 1:
2082          // DataConstant
2083          if (dataVec.extractData(archiveFile,noValues)) {
2084            throw DataException("extractData Error: problem reading data from archive file");
2085          }
2086          break;
2087        case 2:
2088          // DataTagged
2089          if (dataVec.extractData(archiveFile,noValues)) {
2090            throw DataException("extractData Error: problem reading data from archive file");
2091          }
2092          break;
2093        case 3:
2094          // DataExpanded
2095          if (dataVec.extractData(archiveFile,noValues)) {
2096            throw DataException("extractData Error: problem reading data from archive file");
2097          }
2098          break;
2099      }
2100    
2101      if (!archiveFile.good()) {
2102        throw DataException("extractData Error: problem reading from archive file");
2103      }
2104    
2105      //
2106      // Close archive file
2107      archiveFile.close();
2108    
2109      if (!archiveFile.good()) {
2110        throw DataException("extractData Error: problem closing archive file");
2111      }
2112    
2113      //
2114      // Construct an appropriate Data object
2115      DataAbstract* tempData;
2116      switch (dataType) {
2117        case 0:
2118          // DataEmpty
2119          tempData=new DataEmpty();
2120          break;
2121        case 1:
2122          // DataConstant
2123          tempData=new DataConstant(fspace,dataPointShape,dataVec);
2124          break;
2125        case 2:
2126          // DataTagged
2127          tempData=new DataTagged(fspace,dataPointShape,tagNumbers,dataVec);
2128          break;
2129        case 3:
2130          // DataExpanded
2131          tempData=new DataExpanded(fspace,dataPointShape,dataVec);
2132          break;
2133      }
2134      shared_ptr<DataAbstract> temp_data(tempData);
2135      m_data=temp_data;
2136    }
2137    
2138  ostream& escript::operator<<(ostream& o, const Data& data)  ostream& escript::operator<<(ostream& o, const Data& data)
2139  {  {
2140    o << data.toString();    o << data.toString();
2141    return o;    return o;
2142  }  }
2143    
2144    Data
2145    escript::C_GeneralTensorProduct(Data& arg_0,
2146                         Data& arg_1,
2147                         int axis_offset,
2148                         int transpose)
2149    {
2150      // General tensor product: res(SL x SR) = arg_0(SL x SM) * arg_1(SM x SR)
2151      // SM is the product of the last axis_offset entries in arg_0.getShape().
2152    
2153    
2154      // Interpolate if necessary and find an appropriate function space
2155      Data arg_0_Z, arg_1_Z;
2156      if (arg_0.getFunctionSpace()!=arg_1.getFunctionSpace()) {
2157        if (arg_0.probeInterpolation(arg_1.getFunctionSpace())) {
2158          arg_0_Z = arg_0.interpolate(arg_1.getFunctionSpace());
2159          arg_1_Z = Data(arg_1);
2160        }
2161        else if (arg_1.probeInterpolation(arg_0.getFunctionSpace())) {
2162          arg_1_Z=arg_1.interpolate(arg_0.getFunctionSpace());
2163          arg_0_Z =Data(arg_0);
2164        }
2165        else {
2166          throw DataException("Error - C_GeneralTensorProduct: arguments have incompatible function spaces.");
2167        }
2168      } else {
2169          arg_0_Z = Data(arg_0);
2170          arg_1_Z = Data(arg_1);
2171      }
2172      // Get rank and shape of inputs
2173      int rank0 = arg_0_Z.getDataPointRank();
2174      int rank1 = arg_1_Z.getDataPointRank();
2175      DataArrayView::ShapeType shape0 = arg_0_Z.getDataPointShape();
2176      DataArrayView::ShapeType shape1 = arg_1_Z.getDataPointShape();
2177    
2178      // Prepare for the loops of the product and verify compatibility of shapes
2179      int start0=0, start1=0;
2180      if (transpose == 0)       {}
2181      else if (transpose == 1)  { start0 = axis_offset; }
2182      else if (transpose == 2)  { start1 = rank1-axis_offset; }
2183      else              { throw DataException("C_GeneralTensorProduct: Error - transpose should be 0, 1 or 2"); }
2184    
2185      // Adjust the shapes for transpose
2186      DataArrayView::ShapeType tmpShape0;
2187      DataArrayView::ShapeType tmpShape1;
2188      for (int i=0; i<rank0; i++)   { tmpShape0.push_back( shape0[(i+start0)%rank0] ); }
2189      for (int i=0; i<rank1; i++)   { tmpShape1.push_back( shape1[(i+start1)%rank1] ); }
2190    
2191    #if 0
2192      // For debugging: show shape after transpose
2193      char tmp[100];
2194      std::string shapeStr;
2195      shapeStr = "(";
2196      for (int i=0; i<rank0; i++)   { sprintf(tmp, "%d,", tmpShape0[i]); shapeStr += tmp; }
2197      shapeStr += ")";
2198      cout << "C_GeneralTensorProduct: Shape of arg0 is " << shapeStr << endl;
2199      shapeStr = "(";
2200      for (int i=0; i<rank1; i++)   { sprintf(tmp, "%d,", tmpShape1[i]); shapeStr += tmp; }
2201      shapeStr += ")";
2202      cout << "C_GeneralTensorProduct: Shape of arg1 is " << shapeStr << endl;
2203    #endif
2204    
2205      // Prepare for the loops of the product
2206      int SL=1, SM=1, SR=1;
2207      for (int i=0; i<rank0-axis_offset; i++)   {
2208        SL *= tmpShape0[i];
2209      }
2210      for (int i=rank0-axis_offset; i<rank0; i++)   {
2211        if (tmpShape0[i] != tmpShape1[i-(rank0-axis_offset)]) {
2212          throw DataException("C_GeneralTensorProduct: Error - incompatible shapes");
2213        }
2214        SM *= tmpShape0[i];
2215      }
2216      for (int i=axis_offset; i<rank1; i++)     {
2217        SR *= tmpShape1[i];
2218      }
2219    
2220      // Define the shape of the output
2221      DataArrayView::ShapeType shape2;
2222      for (int i=0; i<rank0-axis_offset; i++) { shape2.push_back(tmpShape0[i]); } // First part of arg_0_Z
2223      for (int i=axis_offset; i<rank1; i++)   { shape2.push_back(tmpShape1[i]); } // Last part of arg_1_Z
2224    
2225      // Declare output Data object
2226      Data res;
2227    
2228      if      (arg_0_Z.isConstant()   && arg_1_Z.isConstant()) {
2229        res = Data(0.0, shape2, arg_1_Z.getFunctionSpace());    // DataConstant output
2230        double *ptr_0 = &((arg_0_Z.getPointDataView().getData())[0]);
2231        double *ptr_1 = &((arg_1_Z.getPointDataView().getData())[0]);
2232        double *ptr_2 = &((res.getPointDataView().getData())[0]);
2233        matrix_matrix_product(SL, SM, SR, ptr_0, ptr_1, ptr_2, transpose);
2234      }
2235      else if (arg_0_Z.isConstant()   && arg_1_Z.isTagged()) {
2236    
2237        // Prepare the DataConstant input
2238        DataConstant* tmp_0=dynamic_cast<DataConstant*>(arg_0_Z.borrowData());
2239        if (tmp_0==0) { throw DataException("GTP Programming error - casting to DataConstant."); }
2240    
2241        // Borrow DataTagged input from Data object
2242        DataTagged* tmp_1=dynamic_cast<DataTagged*>(arg_1_Z.borrowData());
2243        if (tmp_1==0) { throw DataException("GTP_1 Programming error - casting to DataTagged."); }
2244    
2245        // Prepare a DataTagged output 2
2246        res = Data(0.0, shape2, arg_1_Z.getFunctionSpace());    // DataTagged output
2247        res.tag();
2248        DataTagged* tmp_2=dynamic_cast<DataTagged*>(res.borrowData());
2249        if (tmp_2==0) { throw DataException("GTP Programming error - casting to DataTagged."); }
2250    
2251        // Prepare offset into DataConstant
2252        int offset_0 = tmp_0->getPointOffset(0,0);
2253        double *ptr_0 = &((arg_0_Z.getPointDataView().getData())[offset_0]);
2254        // Get the views
2255        DataArrayView view_1 = tmp_1->getDefaultValue();
2256        DataArrayView view_2 = tmp_2->getDefaultValue();
2257        // Get the pointers to the actual data
2258        double *ptr_1 = &((view_1.getData())[0]);
2259        double *ptr_2 = &((view_2.getData())[0]);
2260        // Compute an MVP for the default
2261        matrix_matrix_product(SL, SM, SR, ptr_0, ptr_1, ptr_2, transpose);
2262        // Compute an MVP for each tag
2263        const DataTagged::DataMapType& lookup_1=tmp_1->getTagLookup();
2264        DataTagged::DataMapType::const_iterator i; // i->first is a tag, i->second is an offset into memory
2265        for (i=lookup_1.begin();i!=lookup_1.end();i++) {
2266          tmp_2->addTaggedValue(i->first,tmp_2->getDefaultValue());
2267          DataArrayView view_1 = tmp_1->getDataPointByTag(i->first);
2268          DataArrayView view_2 = tmp_2->getDataPointByTag(i->first);
2269          double *ptr_1 = &view_1.getData(0);
2270          double *ptr_2 = &view_2.getData(0);
2271          matrix_matrix_product(SL, SM, SR, ptr_0, ptr_1, ptr_2, transpose);
2272        }
2273    
2274      }
2275      else if (arg_0_Z.isConstant()   && arg_1_Z.isExpanded()) {
2276    
2277        res = Data(0.0, shape2, arg_1_Z.getFunctionSpace(),true); // DataExpanded output
2278        DataConstant* tmp_0=dynamic_cast<DataConstant*>(arg_0_Z.borrowData());
2279        DataExpanded* tmp_1=dynamic_cast<DataExpanded*>(arg_1_Z.borrowData());
2280        DataExpanded* tmp_2=dynamic_cast<DataExpanded*>(res.borrowData());
2281        if (tmp_0==0) { throw DataException("GTP Programming error - casting to DataConstant."); }
2282        if (tmp_1==0) { throw DataException("GTP Programming error - casting to DataExpanded."); }
2283        if (tmp_2==0) { throw DataException("GTP Programming error - casting to DataExpanded."); }
2284        int sampleNo_1,dataPointNo_1;
2285        int numSamples_1 = arg_1_Z.getNumSamples();
2286        int numDataPointsPerSample_1 = arg_1_Z.getNumDataPointsPerSample();
2287        int offset_0 = tmp_0->getPointOffset(0,0);
2288        #pragma omp parallel for private(sampleNo_1,dataPointNo_1) schedule(static)
2289        for (sampleNo_1 = 0; sampleNo_1 < numSamples_1; sampleNo_1++) {
2290          for (dataPointNo_1 = 0; dataPointNo_1 < numDataPointsPerSample_1; dataPointNo_1++) {
2291            int offset_1 = tmp_1->getPointOffset(sampleNo_1,dataPointNo_1);
2292            int offset_2 = tmp_2->getPointOffset(sampleNo_1,dataPointNo_1);
2293            double *ptr_0 = &((arg_0_Z.getPointDataView().getData())[offset_0]);
2294            double *ptr_1 = &((arg_1_Z.getPointDataView().getData())[offset_1]);
2295            double *ptr_2 = &((res.getPointDataView().getData())[offset_2]);
2296            matrix_matrix_product(SL, SM, SR, ptr_0, ptr_1, ptr_2, transpose);
2297          }
2298        }
2299    
2300      }
2301      else if (arg_0_Z.isTagged()     && arg_1_Z.isConstant()) {
2302    
2303        // Borrow DataTagged input from Data object
2304        DataTagged* tmp_0=dynamic_cast<DataTagged*>(arg_0_Z.borrowData());
2305        if (tmp_0==0) { throw DataException("GTP_0 Programming error - casting to DataTagged."); }
2306    
2307        // Prepare the DataConstant input
2308        DataConstant* tmp_1=dynamic_cast<DataConstant*>(arg_1_Z.borrowData());
2309        if (tmp_1==0) { throw DataException("GTP Programming error - casting to DataConstant."); }
2310    
2311        // Prepare a DataTagged output 2
2312        res = Data(0.0, shape2, arg_0_Z.getFunctionSpace());    // DataTagged output
2313        res.tag();
2314        DataTagged* tmp_2=dynamic_cast<DataTagged*>(res.borrowData());
2315        if (tmp_2==0) { throw DataException("GTP Programming error - casting to DataTagged."); }
2316    
2317        // Prepare offset into DataConstant
2318        int offset_1 = tmp_1->getPointOffset(0,0);
2319        double *ptr_1 = &((arg_1_Z.getPointDataView().getData())[offset_1]);
2320        // Get the views
2321        DataArrayView view_0 = tmp_0->getDefaultValue();
2322        DataArrayView view_2 = tmp_2->getDefaultValue();
2323        // Get the pointers to the actual data
2324        double *ptr_0 = &((view_0.getData())[0]);
2325        double *ptr_2 = &((view_2.getData())[0]);
2326        // Compute an MVP for the default
2327        matrix_matrix_product(SL, SM, SR, ptr_0, ptr_1, ptr_2, transpose);
2328        // Compute an MVP for each tag
2329        const DataTagged::DataMapType& lookup_0=tmp_0->getTagLookup();
2330        DataTagged::DataMapType::const_iterator i; // i->first is a tag, i->second is an offset into memory
2331        for (i=lookup_0.begin();i!=lookup_0.end();i++) {
2332          tmp_2->addTaggedValue(i->first,tmp_2->getDefaultValue());
2333          DataArrayView view_0 = tmp_0->getDataPointByTag(i->first);
2334          DataArrayView view_2 = tmp_2->getDataPointByTag(i->first);
2335          double *ptr_0 = &view_0.getData(0);
2336          double *ptr_2 = &view_2.getData(0);
2337          matrix_matrix_product(SL, SM, SR, ptr_0, ptr_1, ptr_2, transpose);
2338        }
2339    
2340      }
2341      else if (arg_0_Z.isTagged()     && arg_1_Z.isTagged()) {
2342    
2343        // Borrow DataTagged input from Data object
2344        DataTagged* tmp_0=dynamic_cast<DataTagged*>(arg_0_Z.borrowData());
2345        if (tmp_0==0) { throw DataException("GTP Programming error - casting to DataTagged."); }
2346    
2347        // Borrow DataTagged input from Data object
2348        DataTagged* tmp_1=dynamic_cast<DataTagged*>(arg_1_Z.borrowData());
2349        if (tmp_1==0) { throw DataException("GTP Programming error - casting to DataTagged."); }
2350    
2351        // Prepare a DataTagged output 2
2352        res = Data(0.0, shape2, arg_1_Z.getFunctionSpace());
2353        res.tag();  // DataTagged output
2354        DataTagged* tmp_2=dynamic_cast<DataTagged*>(res.borrowData());
2355        if (tmp_2==0) { throw DataException("GTP Programming error - casting to DataTagged."); }
2356    
2357        // Get the views
2358        DataArrayView view_0 = tmp_0->getDefaultValue();
2359        DataArrayView view_1 = tmp_1->getDefaultValue();
2360        DataArrayView view_2 = tmp_2->getDefaultValue();
2361        // Get the pointers to the actual data
2362        double *ptr_0 = &((view_0.getData())[0]);
2363        double *ptr_1 = &((view_1.getData())[0]);
2364        double *ptr_2 = &((view_2.getData())[0]);
2365        // Compute an MVP for the default
2366        matrix_matrix_product(SL, SM, SR, ptr_0, ptr_1, ptr_2, transpose);
2367        // Merge the tags
2368        DataTagged::DataMapType::const_iterator i; // i->first is a tag, i->second is an offset into memory
2369        const DataTagged::DataMapType& lookup_0=tmp_0->getTagLookup();
2370        const DataTagged::DataMapType& lookup_1=tmp_1->getTagLookup();
2371        for (i=lookup_0.begin();i!=lookup_0.end();i++) {
2372          tmp_2->addTaggedValue(i->first,tmp_2->getDefaultValue()); // use tmp_2 to get correct shape
2373        }
2374        for (i=lookup_1.begin();i!=lookup_1.end();i++) {
2375          tmp_2->addTaggedValue(i->first,tmp_2->getDefaultValue());
2376        }
2377        // Compute an MVP for each tag
2378        const DataTagged::DataMapType& lookup_2=tmp_2->getTagLookup();
2379        for (i=lookup_2.begin();i!=lookup_2.end();i++) {
2380          DataArrayView view_0 = tmp_0->getDataPointByTag(i->first);
2381          DataArrayView view_1 = tmp_1->getDataPointByTag(i->first);
2382          DataArrayView view_2 = tmp_2->getDataPointByTag(i->first);
2383          double *ptr_0 = &view_0.getData(0);
2384          double *ptr_1 = &view_1.getData(0);
2385          double *ptr_2 = &view_2.getData(0);
2386          matrix_matrix_product(SL, SM, SR, ptr_0, ptr_1, ptr_2, transpose);
2387        }
2388    
2389      }
2390      else if (arg_0_Z.isTagged()     && arg_1_Z.isExpanded()) {
2391    
2392        // After finding a common function space above the two inputs have the same numSamples and num DPPS
2393        res = Data(0.0, shape2, arg_1_Z.getFunctionSpace(),true); // DataExpanded output
2394        DataTagged*   tmp_0=dynamic_cast<DataTagged*>(arg_0_Z.borrowData());
2395        DataExpanded* tmp_1=dynamic_cast<DataExpanded*>(arg_1_Z.borrowData());
2396        DataExpanded* tmp_2=dynamic_cast<DataExpanded*>(res.borrowData());
2397        if (tmp_0==0) { throw DataException("GTP Programming error - casting to DataTagged."); }
2398        if (tmp_1==0) { throw DataException("GTP Programming error - casting to DataExpanded."); }
2399        if (tmp_2==0) { throw DataException("GTP Programming error - casting to DataExpanded."); }
2400        int sampleNo_0,dataPointNo_0;
2401        int numSamples_0 = arg_0_Z.getNumSamples();
2402        int numDataPointsPerSample_0 = arg_0_Z.getNumDataPointsPerSample();
2403        #pragma omp parallel for private(sampleNo_0,dataPointNo_0) schedule(static)
2404        for (sampleNo_0 = 0; sampleNo_0 < numSamples_0; sampleNo_0++) {
2405          int offset_0 = tmp_0->getPointOffset(sampleNo_0,0); // They're all the same, so just use #0
2406          double *ptr_0 = &((arg_0_Z.getPointDataView().getData())[offset_0]);
2407          for (dataPointNo_0 = 0; dataPointNo_0 < numDataPointsPerSample_0; dataPointNo_0++) {
2408            int offset_1 = tmp_1->getPointOffset(sampleNo_0,dataPointNo_0);
2409            int offset_2 = tmp_2->getPointOffset(sampleNo_0,dataPointNo_0);
2410            double *ptr_1 = &((arg_1_Z.getPointDataView().getData())[offset_1]);
2411            double *ptr_2 = &((res.getPointDataView().getData())[offset_2]);
2412            matrix_matrix_product(SL, SM, SR, ptr_0, ptr_1, ptr_2, transpose);
2413          }
2414        }
2415    
2416      }
2417      else if (arg_0_Z.isExpanded()   && arg_1_Z.isConstant()) {
2418    
2419        res = Data(0.0, shape2, arg_1_Z.getFunctionSpace(),true); // DataExpanded output
2420        DataExpanded* tmp_0=dynamic_cast<DataExpanded*>(arg_0_Z.borrowData());
2421        DataConstant* tmp_1=dynamic_cast<DataConstant*>(arg_1_Z.borrowData());
2422        DataExpanded* tmp_2=dynamic_cast<DataExpanded*>(res.borrowData());
2423        if (tmp_0==0) { throw DataException("GTP Programming error - casting to DataExpanded."); }
2424        if (tmp_1==0) { throw DataException("GTP Programming error - casting to DataConstant."); }
2425        if (tmp_2==0) { throw DataException("GTP Programming error - casting to DataExpanded."); }
2426        int sampleNo_0,dataPointNo_0;
2427        int numSamples_0 = arg_0_Z.getNumSamples();
2428        int numDataPointsPerSample_0 = arg_0_Z.getNumDataPointsPerSample();
2429        int offset_1 = tmp_1->getPointOffset(0,0);
2430        #pragma omp parallel for private(sampleNo_0,dataPointNo_0) schedule(static)
2431        for (sampleNo_0 = 0; sampleNo_0 < numSamples_0; sampleNo_0++) {
2432          for (dataPointNo_0 = 0; dataPointNo_0 < numDataPointsPerSample_0; dataPointNo_0++) {
2433            int offset_0 = tmp_0->getPointOffset(sampleNo_0,dataPointNo_0);
2434            int offset_2 = tmp_2->getPointOffset(sampleNo_0,dataPointNo_0);
2435            double *ptr_0 = &((arg_0_Z.getPointDataView().getData())[offset_0]);
2436            double *ptr_1 = &((arg_1_Z.getPointDataView().getData())[offset_1]);
2437            double *ptr_2 = &((res.getPointDataView().getData())[offset_2]);
2438            matrix_matrix_product(SL, SM, SR, ptr_0, ptr_1, ptr_2, transpose);
2439          }
2440        }
2441    
2442    
2443      }
2444      else if (arg_0_Z.isExpanded()   && arg_1_Z.isTagged()) {
2445    
2446        // After finding a common function space above the two inputs have the same numSamples and num DPPS
2447        res = Data(0.0, shape2, arg_1_Z.getFunctionSpace(),true); // DataExpanded output
2448        DataExpanded* tmp_0=dynamic_cast<DataExpanded*>(arg_0_Z.borrowData());
2449        DataTagged*   tmp_1=dynamic_cast<DataTagged*>(arg_1_Z.borrowData());
2450        DataExpanded* tmp_2=dynamic_cast<DataExpanded*>(res.borrowData());
2451        if (tmp_0==0) { throw DataException("GTP Programming error - casting to DataExpanded."); }
2452        if (tmp_1==0) { throw DataException("GTP Programming error - casting to DataTagged."); }
2453        if (tmp_2==0) { throw DataException("GTP Programming error - casting to DataExpanded."); }
2454        int sampleNo_0,dataPointNo_0;
2455        int numSamples_0 = arg_0_Z.getNumSamples();
2456        int numDataPointsPerSample_0 = arg_0_Z.getNumDataPointsPerSample();
2457        #pragma omp parallel for private(sampleNo_0,dataPointNo_0) schedule(static)
2458        for (sampleNo_0 = 0; sampleNo_0 < numSamples_0; sampleNo_0++) {
2459          int offset_1 = tmp_1->getPointOffset(sampleNo_0,0);
2460          double *ptr_1 = &((arg_1_Z.getPointDataView().getData())[offset_1]);
2461          for (dataPointNo_0 = 0; dataPointNo_0 < numDataPointsPerSample_0; dataPointNo_0++) {
2462            int offset_0 = tmp_0->getPointOffset(sampleNo_0,dataPointNo_0);
2463            int offset_2 = tmp_2->getPointOffset(sampleNo_0,dataPointNo_0);
2464            double *ptr_0 = &((arg_0_Z.getPointDataView().getData())[offset_0]);
2465            double *ptr_2 = &((res.getPointDataView().getData())[offset_2]);
2466            matrix_matrix_product(SL, SM, SR, ptr_0, ptr_1, ptr_2, transpose);
2467          }
2468        }
2469    
2470      }
2471      else if (arg_0_Z.isExpanded()   && arg_1_Z.isExpanded()) {
2472    
2473        // After finding a common function space above the two inputs have the same numSamples and num DPPS
2474        res = Data(0.0, shape2, arg_1_Z.getFunctionSpace(),true); // DataExpanded output
2475        DataExpanded* tmp_0=dynamic_cast<DataExpanded*>(arg_0_Z.borrowData());
2476        DataExpanded* tmp_1=dynamic_cast<DataExpanded*>(arg_1_Z.borrowData());
2477        DataExpanded* tmp_2=dynamic_cast<DataExpanded*>(res.borrowData());
2478        if (tmp_0==0) { throw DataException("GTP Programming error - casting to DataExpanded."); }
2479        if (tmp_1==0) { throw DataException("GTP Programming error - casting to DataExpanded."); }
2480        if (tmp_2==0) { throw DataException("GTP Programming error - casting to DataExpanded."); }
2481        int sampleNo_0,dataPointNo_0;
2482        int numSamples_0 = arg_0_Z.getNumSamples();
2483        int numDataPointsPerSample_0 = arg_0_Z.getNumDataPointsPerSample();
2484        #pragma omp parallel for private(sampleNo_0,dataPointNo_0) schedule(static)
2485        for (sampleNo_0 = 0; sampleNo_0 < numSamples_0; sampleNo_0++) {
2486          for (dataPointNo_0 = 0; dataPointNo_0 < numDataPointsPerSample_0; dataPointNo_0++) {
2487            int offset_0 = tmp_0->getPointOffset(sampleNo_0,dataPointNo_0);
2488            int offset_1 = tmp_1->getPointOffset(sampleNo_0,dataPointNo_0);
2489            int offset_2 = tmp_2->getPointOffset(sampleNo_0,dataPointNo_0);
2490            double *ptr_0 = &((arg_0_Z.getPointDataView().getData())[offset_0]);
2491            double *ptr_1 = &((arg_1_Z.getPointDataView().getData())[offset_1]);
2492            double *ptr_2 = &((res.getPointDataView().getData())[offset_2]);
2493            matrix_matrix_product(SL, SM, SR, ptr_0, ptr_1, ptr_2, transpose);
2494          }
2495        }
2496    
2497      }
2498      else {
2499        throw DataException("Error - C_GeneralTensorProduct: unknown combination of inputs");
2500      }
2501    
2502      return res;
2503    }
2504    
2505    DataAbstract*
2506    Data::borrowData() const
2507    {
2508      return m_data.get();
2509    }
2510    
2511    /* Member functions specific to the MPI implementation */
2512    
2513    void
2514    Data::print()
2515    {
2516      int i,j;
2517      
2518      printf( "Data is %dX%d\n", getNumSamples(), getNumDataPointsPerSample() );
2519      for( i=0; i<getNumSamples(); i++ )
2520      {
2521        printf( "[%6d]", i );
2522        for( j=0; j<getNumDataPointsPerSample(); j++ )
2523          printf( "\t%10.7g", (getSampleData(i))[j] );
2524        printf( "\n" );
2525      }
2526    }
2527    
2528    int
2529    Data::get_MPISize() const
2530    {
2531        int error, size;
2532    #ifdef PASO_MPI
2533        error = MPI_Comm_size( get_MPIComm(), &size );
2534    #else
2535        size = 1;
2536    #endif
2537        return size;
2538    }
2539    
2540    int
2541    Data::get_MPIRank() const
2542    {
2543        int error, rank;
2544    #ifdef PASO_MPI
2545        error = MPI_Comm_rank( get_MPIComm(), &rank );
2546    #else
2547        rank = 0;
2548    #endif
2549        return rank;
2550    }
2551    
2552    MPI_Comm
2553    Data::get_MPIComm() const
2554    {
2555    #ifdef PASO_MPI
2556        return MPI_COMM_WORLD;
2557    #else
2558        return -1;
2559    #endif
2560    }
2561    

Legend:
Removed from v.97  
changed lines
  Added in v.1092

  ViewVC Help
Powered by ViewVC 1.1.26