/[escript]/branches/clazy/escriptcore/src/DataLazy.cpp
ViewVC logotype

Diff of /branches/clazy/escriptcore/src/DataLazy.cpp

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

revision 2195 by jfenwick, Wed Jan 7 04:13:52 2009 UTC revision 2199 by jfenwick, Thu Jan 8 06:10:52 2009 UTC
# Line 28  Line 28 
28  #include "UnaryFuncs.h"     // for escript::fsign  #include "UnaryFuncs.h"     // for escript::fsign
29  #include "Utils.h"  #include "Utils.h"
30    
31    #include <iomanip>      // for some fancy formatting in debug
32    
33  // #define LAZYDEBUG(X) if (privdebug){X;}  // #define LAZYDEBUG(X) if (privdebug){X;}
34  #define LAZYDEBUG(X)  #define LAZYDEBUG(X)
35  namespace  namespace
# Line 198  resultShape(DataAbstract_ptr left, DataA Line 200  resultShape(DataAbstract_ptr left, DataA
200  // return the shape for "op left"  // return the shape for "op left"
201    
202  DataTypes::ShapeType  DataTypes::ShapeType
203  resultShape(DataAbstract_ptr left, ES_optype op)  resultShape(DataAbstract_ptr left, ES_optype op, int axis_offset)
204  {  {
205      switch(op)      switch(op)
206      {      {
207          case TRANS:          case TRANS:
208          return left->getShape();         {            // for the scoping of variables
209            const DataTypes::ShapeType& s=left->getShape();
210            DataTypes::ShapeType sh;
211            int rank=left->getRank();
212            if (axis_offset<0 || axis_offset>rank)
213            {
214                   throw DataException("Error - Data::transpose must have 0 <= axis_offset <= rank=" + rank);
215                }
216                for (int i=0; i<rank; i++)
217            {
218               int index = (axis_offset+i)%rank;
219                   sh.push_back(s[index]); // Append to new shape
220                }
221            return sh;
222           }
223      break;      break;
224      case TRACE:      case TRACE:
225          return DataTypes::scalarShape;         {
226            int rank=left->getRank();
227            if (rank<2)
228            {
229               throw DataException("Trace can only be computed for objects with rank 2 or greater.");
230            }
231            if ((axis_offset>rank-2) || (axis_offset<0))
232            {
233               throw DataException("Trace: axis offset must lie between 0 and rank-2 inclusive.");
234            }
235            if (rank==2)
236            {
237               return DataTypes::scalarShape;
238            }
239            else if (rank==3)
240            {
241               DataTypes::ShapeType sh;
242                   if (axis_offset==0)
243               {
244                    sh.push_back(left->getShape()[2]);
245                   }
246                   else     // offset==1
247               {
248                sh.push_back(left->getShape()[0]);
249                   }
250               return sh;
251            }
252            else if (rank==4)
253            {
254               DataTypes::ShapeType sh;
255               const DataTypes::ShapeType& s=left->getShape();
256                   if (axis_offset==0)
257               {
258                    sh.push_back(s[2]);
259                    sh.push_back(s[3]);
260                   }
261                   else if (axis_offset==1)
262               {
263                    sh.push_back(s[0]);
264                    sh.push_back(s[3]);
265                   }
266               else     // offset==2
267               {
268                sh.push_back(s[0]);
269                sh.push_back(s[1]);
270               }
271               return sh;
272            }
273            else        // unknown rank
274            {
275               throw DataException("Error - Data::trace can only be calculated for rank 2, 3 or 4 object.");
276            }
277           }
278      break;      break;
279          default:          default:
280      throw DataException("Programmer error - resultShape(left,op) can't compute shapes for operator "+opToString(op)+".");      throw DataException("Programmer error - resultShape(left,op) can't compute shapes for operator "+opToString(op)+".");
# Line 331  DataLazy::DataLazy(DataAbstract_ptr p) Line 399  DataLazy::DataLazy(DataAbstract_ptr p)
399     {     {
400      DataReady_ptr dr=dynamic_pointer_cast<DataReady>(p);      DataReady_ptr dr=dynamic_pointer_cast<DataReady>(p);
401      makeIdentity(dr);      makeIdentity(dr);
402  cout << "Wrapping " << dr.get() << " id=" << m_id.get() << endl;  LAZYDEBUG(cout << "Wrapping " << dr.get() << " id=" << m_id.get() << endl;)
403     }     }
404  LAZYDEBUG(cout << "(1)Lazy created with " << m_samplesize << endl;)  LAZYDEBUG(cout << "(1)Lazy created with " << m_samplesize << endl;)
405  }  }
# Line 377  DataLazy::DataLazy(DataAbstract_ptr left Line 445  DataLazy::DataLazy(DataAbstract_ptr left
445      m_op(op),      m_op(op),
446      m_SL(0), m_SM(0), m_SR(0)      m_SL(0), m_SM(0), m_SR(0)
447  {  {
448  cout << "Forming operator with " << left.get() << " " << right.get() << endl;  LAZYDEBUG(cout << "Forming operator with " << left.get() << " " << right.get() << endl;)
449     if ((getOpgroup(op)!=G_BINARY))     if ((getOpgroup(op)!=G_BINARY))
450     {     {
451      throw DataException("Programmer error - constructor DataLazy(left, right, op) will only process BINARY operations.");      throw DataException("Programmer error - constructor DataLazy(left, right, op) will only process BINARY operations.");
# Line 394  cout << "Forming operator with " << left Line 462  cout << "Forming operator with " << left
462     {     {
463      Data tmp(Data(right),getFunctionSpace());      Data tmp(Data(right),getFunctionSpace());
464      right=tmp.borrowDataPtr();      right=tmp.borrowDataPtr();
465  cout << "Right interpolation required " << right.get() << endl;  LAZYDEBUG(cout << "Right interpolation required " << right.get() << endl;)
466     }     }
467     left->operandCheck(*right);     left->operandCheck(*right);
468    
469     if (left->isLazy())          // the children need to be DataLazy. Wrap them in IDENTITY if required     if (left->isLazy())          // the children need to be DataLazy. Wrap them in IDENTITY if required
470     {     {
471      m_left=dynamic_pointer_cast<DataLazy>(left);      m_left=dynamic_pointer_cast<DataLazy>(left);
472  cout << "Left is " << m_left->toString() << endl;  LAZYDEBUG(cout << "Left is " << m_left->toString() << endl;)
473     }     }
474     else     else
475     {     {
476      m_left=DataLazy_ptr(new DataLazy(left));      m_left=DataLazy_ptr(new DataLazy(left));
477  cout << "Left " << left.get() << " wrapped " << m_left->m_id.get() << endl;  LAZYDEBUG(cout << "Left " << left.get() << " wrapped " << m_left->m_id.get() << endl;)
478     }     }
479     if (right->isLazy())     if (right->isLazy())
480     {     {
481      m_right=dynamic_pointer_cast<DataLazy>(right);      m_right=dynamic_pointer_cast<DataLazy>(right);
482  cout << "Right is " << m_right->toString() << endl;  LAZYDEBUG(cout << "Right is " << m_right->toString() << endl;)
483     }     }
484     else     else
485     {     {
486      m_right=DataLazy_ptr(new DataLazy(right));      m_right=DataLazy_ptr(new DataLazy(right));
487  cout << "Right " << right.get() << " wrapped " << m_right->m_id.get() << endl;  LAZYDEBUG(cout << "Right " << right.get() << " wrapped " << m_right->m_id.get() << endl;)
488     }     }
489     char lt=m_left->m_readytype;     char lt=m_left->m_readytype;
490     char rt=m_right->m_readytype;     char rt=m_right->m_readytype;
# Line 510  LAZYDEBUG(cout << "(4)Lazy created with Line 578  LAZYDEBUG(cout << "(4)Lazy created with
578    
579    
580  DataLazy::DataLazy(DataAbstract_ptr left, ES_optype op, int axis_offset)  DataLazy::DataLazy(DataAbstract_ptr left, ES_optype op, int axis_offset)
581      : parent(left->getFunctionSpace(), resultShape(left,op)),      : parent(left->getFunctionSpace(), resultShape(left,op, axis_offset)),
582      m_op(op),      m_op(op),
583      m_axis_offset(axis_offset),      m_axis_offset(axis_offset),
584      m_transpose(0),      m_transpose(0),
# Line 1184  LAZYDEBUG(cout << " numsteps=" << numste Line 1252  LAZYDEBUG(cout << " numsteps=" << numste
1252  LAZYDEBUG(cout << "onumsteps=" << onumsteps << endl;)  LAZYDEBUG(cout << "onumsteps=" << onumsteps << endl;)
1253  LAZYDEBUG(cout << " DPPS=" << m_left->getNumDPPSample() << "," <<m_right->getNumDPPSample() << endl;)  LAZYDEBUG(cout << " DPPS=" << m_left->getNumDPPSample() << "," <<m_right->getNumDPPSample() << endl;)
1254  LAZYDEBUG(cout << "" << LS << RS << LN << RN << LES << RES <<LEN << REN <<   endl;)  LAZYDEBUG(cout << "" << LS << RS << LN << RN << LES << RES <<LEN << REN <<   endl;)
1255  // LAZYDEBUG(  
 // cout << "Results of bin" << endl;  
 // cout << "Left=";  
 // for (int i=lroffset;i<lroffset+m_left->m_samplesize;++i)  
 // {  
 // cout << (*left)[i] << " ";  
 // }  
 // cout << endl << "Right=";  
 // for (int i=rroffset;i<rroffset+(m_right->m_readytype=='E'?m_right->m_samplesize:m_right->getNoValues());++i)  
 // {  
 // cout << (*right)[i] << " ";  
 // }  
 // cout << endl;  
 // )  
1256    
1257    double* resultp=&(v[offset]);     // results are stored at the vector offset we recieved    double* resultp=&(v[offset]);     // results are stored at the vector offset we recieved
1258    switch(m_op)    switch(m_op)
# Line 1257  LAZYDEBUG(cout << "Resolve TensorProduct Line 1312  LAZYDEBUG(cout << "Resolve TensorProduct
1312    int rightStep=(rightExp?m_right->getNoValues() : 0);    int rightStep=(rightExp?m_right->getNoValues() : 0);
1313    
1314    int resultStep=getNoValues();    int resultStep=getNoValues();
 //   int resultStep=max(leftStep,rightStep);    // only one (at most) should be !=0  
1315      // Get the values of sub-expressions (leave a gap of one sample for the result).      // Get the values of sub-expressions (leave a gap of one sample for the result).
1316    int gap=offset+m_left->getMaxSampleSize();    // actually only needs to be m_left->m_samplesize    int gap=offset+m_samplesize;  
1317    
1318  LAZYDEBUG(cout << "Query left with offset=" << gap << endl;)  LAZYDEBUG(cout << "Query left with offset=" << gap << endl;)
1319    
1320    const ValueType* left=m_left->resolveSample(v,gap,sampleNo,lroffset);    const ValueType* left=m_left->resolveSample(v,gap,sampleNo,lroffset);
1321    gap+=m_right->getMaxSampleSize();    gap+=m_left->getMaxSampleSize();
1322    
1323    
1324  LAZYDEBUG(cout << "Query right with offset=" << gap << endl;)  LAZYDEBUG(cout << "Query right with offset=" << gap << endl;)
# Line 1276  LAZYDEBUG(cerr << "[Left shape]=" << Dat Line 1330  LAZYDEBUG(cerr << "[Left shape]=" << Dat
1330  cout << getNoValues() << endl;)  cout << getNoValues() << endl;)
1331  LAZYDEBUG(cerr << "Result of left=";)  LAZYDEBUG(cerr << "Result of left=";)
1332  LAZYDEBUG(cerr << "[" << lroffset << " .. " << lroffset+m_left->getNoValues() << "]" << endl;  LAZYDEBUG(cerr << "[" << lroffset << " .. " << lroffset+m_left->getNoValues() << "]" << endl;
1333  for (int i=lroffset;i<lroffset+m_left->getNoValues();++i)  
1334    for (int i=lroffset, limit=lroffset+(leftExp?m_left->getNoValues()*m_left->getNumDPPSample():m_left->getNoValues());i<limit;++i)
1335  {  {
1336  cout << (*left)[i] << " ";  cout << "[" << setw(2) << i-lroffset << "] " << setw(10) << (*left)[i] << " ";
1337    if (i%4==0) cout << endl;
1338  })  })
1339  LAZYDEBUG(cerr << "\nResult of right=" << endl;)  LAZYDEBUG(cerr << "\nResult of right=" << endl;)
1340  LAZYDEBUG(cerr << "[" << rroffset << " .. " << rroffset+m_right->m_samplesize << "]" << endl;  LAZYDEBUG(
1341  for (int i=rroffset;i<rroffset+m_right->m_samplesize;++i)  for (int i=rroffset, limit=rroffset+(rightExp?m_right->getNoValues()*m_right->getNumDPPSample():m_right->getNoValues());i<limit;++i)
1342  {  {
1343  cerr << (*right)[i] << " ";  cerr << "[" <<  setw(2)<< i-rroffset << "] " << setw(10) << (*right)[i] << " ";
1344    if (i%4==0) cout << endl;
1345  }  }
1346  cerr << endl;  cerr << endl;
1347  )  )
# Line 1302  LAZYDEBUG(cout << "DPPS=" << m_right->ge Line 1359  LAZYDEBUG(cout << "DPPS=" << m_right->ge
1359      case PROD:      case PROD:
1360      for (int i=0;i<steps;++i,resultp+=resultStep)      for (int i=0;i<steps;++i,resultp+=resultStep)
1361      {      {
1362    
1363  LAZYDEBUG(cout << "lroffset=" << lroffset << "rroffset=" << rroffset << endl;)  LAZYDEBUG(cout << "lroffset=" << lroffset << "rroffset=" << rroffset << endl;)
1364  LAZYDEBUG(cout << "l*=" << left << " r*=" << right << endl;)  LAZYDEBUG(cout << "l*=" << left << " r*=" << right << endl;)
1365  LAZYDEBUG(cout << "m_SL=" << m_SL << " m_SM=" << m_SM << " m_SR=" << m_SR << endl;)  LAZYDEBUG(cout << "m_SL=" << m_SL << " m_SM=" << m_SM << " m_SR=" << m_SR << endl;)
1366    
1367            const double *ptr_0 = &((*left)[lroffset]);            const double *ptr_0 = &((*left)[lroffset]);
1368            const double *ptr_1 = &((*right)[rroffset]);            const double *ptr_1 = &((*right)[rroffset]);
1369    
1370  LAZYDEBUG(cout << DataTypes::pointToString(*left, m_left->getShape(),lroffset,"LEFT") << endl;)  LAZYDEBUG(cout << DataTypes::pointToString(*left, m_left->getShape(),lroffset,"LEFT") << endl;)
1371  LAZYDEBUG(cout << DataTypes::pointToString(*right,m_right->getShape(),rroffset, "RIGHT") << endl;)  LAZYDEBUG(cout << DataTypes::pointToString(*right,m_right->getShape(),rroffset, "RIGHT") << endl;)
1372    
1373            matrix_matrix_product(m_SL, m_SM, m_SR, ptr_0, ptr_1, resultp, m_transpose);            matrix_matrix_product(m_SL, m_SM, m_SR, ptr_0, ptr_1, resultp, m_transpose);
1374    
1375  LAZYDEBUG(cout << "Results--\n";  LAZYDEBUG(cout << "Results--\n";
1376    {
1377      DataVector dv(getNoValues());
1378  for (int z=0;z<getNoValues();++z)  for (int z=0;z<getNoValues();++z)
1379  {  {
1380  cout << resultp[z] << " ";    cout << "[" << setw(2) << z<< "] " << setw(10) << resultp[z] << " ";
1381      if (z%4==0) cout << endl;
1382      dv[z]=resultp[z];
1383  }  }
1384    cout << endl << DataTypes::pointToString(dv,getShape(),0,"RESLT");
1385  cout << "\nWritten to: " << resultp << " resultStep=" << resultStep << endl;  cout << "\nWritten to: " << resultp << " resultStep=" << resultStep << endl;
1386    }
1387  )  )
1388        lroffset+=leftStep;        lroffset+=leftStep;
1389        rroffset+=rightStep;        rroffset+=rightStep;

Legend:
Removed from v.2195  
changed lines
  Added in v.2199

  ViewVC Help
Powered by ViewVC 1.1.26