/[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

revision 790 by bcumming, Wed Jul 26 23:12:34 2006 UTC revision 826 by gross, Tue Aug 29 09:12:33 2006 UTC
# Line 1316  Data::minval() const Line 1316  Data::minval() const
1316  }  }
1317    
1318  Data  Data
1319  Data::trace() const  Data::swapaxes(const int axis0, const int axis1) const
1320  {  {
1321  #if defined DOPROF       int axis0_tmp,axis1_tmp;
1322    profData->reduction2++;       #if defined DOPROF
1323  #endif       profData->unary++;
1324    Trace trace_func;       #endif
1325    return dp_algorithm(trace_func,0);       DataArrayView::ShapeType s=getDataPointShape();
1326         DataArrayView::ShapeType ev_shape;
1327         // Here's the equivalent of python s_out=s[axis_offset:]+s[:axis_offset]
1328         // which goes thru all shape vector elements starting with axis_offset (at index=rank wrap around to 0)
1329         int rank=getDataPointRank();
1330         if (rank<2) {
1331            throw DataException("Error - Data::swapaxes argument must have at least rank 2.");
1332         }
1333         if (axis0<0 || axis0>rank-1) {
1334            throw DataException("Error - Data::swapaxes: axis0 must be between 0 and rank-1=" + rank-1);
1335         }
1336         if (axis1<0 || axis1>rank-1) {
1337             throw DataException("Error - Data::swapaxes: axis1 must be between 0 and rank-1=" + rank-1);
1338         }
1339         if (axis0 == axis1) {
1340             throw DataException("Error - Data::swapaxes: axis indices must be different.");
1341         }
1342         if (axis0 > axis1) {
1343             axis0_tmp=axis1;
1344             axis1_tmp=axis0;
1345         } else {
1346             axis0_tmp=axis0;
1347             axis1_tmp=axis1;
1348         }
1349         for (int i=0; i<rank; i++) {
1350           if (i == axis0_tmp) {
1351              ev_shape.push_back(s[axis1_tmp]);
1352           } else if (i == axis1_tmp) {
1353              ev_shape.push_back(s[axis0_tmp]);
1354           } else {
1355              ev_shape.push_back(s[i]);
1356           }
1357         }
1358         Data ev(0.,ev_shape,getFunctionSpace());
1359         ev.typeMatchRight(*this);
1360         m_data->swapaxes(ev.m_data.get(), axis0_tmp, axis1_tmp);
1361         return ev;
1362    
1363  }  }
1364    
1365  Data  Data
# Line 1388  Data::nonsymmetric() const Line 1425  Data::nonsymmetric() const
1425  }  }
1426    
1427  Data  Data
1428  Data::matrixtrace(int axis_offset) const  Data::trace(int axis_offset) const
1429  {  {
1430       #if defined DOPROF       #if defined DOPROF
1431          profData->unary++;          profData->unary++;
# Line 1398  Data::matrixtrace(int axis_offset) const Line 1435  Data::matrixtrace(int axis_offset) const
1435          DataArrayView::ShapeType ev_shape;          DataArrayView::ShapeType ev_shape;
1436          Data ev(0.,ev_shape,getFunctionSpace());          Data ev(0.,ev_shape,getFunctionSpace());
1437          ev.typeMatchRight(*this);          ev.typeMatchRight(*this);
1438          m_data->matrixtrace(ev.m_data.get(), axis_offset);          m_data->trace(ev.m_data.get(), axis_offset);
1439          return ev;          return ev;
1440       }       }
1441       if (getDataPointRank()==3) {       if (getDataPointRank()==3) {
# Line 1413  Data::matrixtrace(int axis_offset) const Line 1450  Data::matrixtrace(int axis_offset) const
1450          }          }
1451          Data ev(0.,ev_shape,getFunctionSpace());          Data ev(0.,ev_shape,getFunctionSpace());
1452          ev.typeMatchRight(*this);          ev.typeMatchRight(*this);
1453          m_data->matrixtrace(ev.m_data.get(), axis_offset);          m_data->trace(ev.m_data.get(), axis_offset);
1454          return ev;          return ev;
1455       }       }
1456       if (getDataPointRank()==4) {       if (getDataPointRank()==4) {
# Line 1432  Data::matrixtrace(int axis_offset) const Line 1469  Data::matrixtrace(int axis_offset) const
1469      }      }
1470          Data ev(0.,ev_shape,getFunctionSpace());          Data ev(0.,ev_shape,getFunctionSpace());
1471          ev.typeMatchRight(*this);          ev.typeMatchRight(*this);
1472      m_data->matrixtrace(ev.m_data.get(), axis_offset);      m_data->trace(ev.m_data.get(), axis_offset);
1473          return ev;          return ev;
1474       }       }
1475       else {       else {
1476          throw DataException("Error - Data::matrixtrace can only be calculated for rank 2, 3 or 4 object.");          throw DataException("Error - Data::trace can only be calculated for rank 2, 3 or 4 object.");
1477       }       }
1478  }  }
1479    
1480  Data  Data
1481  Data::transpose(int axis_offset) const  Data::transpose(int axis_offset) const
1482  {  {
1483  #if defined DOPROF       #if defined DOPROF
1484       profData->reduction2++;       profData->unary++;
1485  #endif       #endif
1486       DataArrayView::ShapeType s=getDataPointShape();       DataArrayView::ShapeType s=getDataPointShape();
1487       DataArrayView::ShapeType ev_shape;       DataArrayView::ShapeType ev_shape;
1488       // Here's the equivalent of python s_out=s[axis_offset:]+s[:axis_offset]       // Here's the equivalent of python s_out=s[axis_offset:]+s[:axis_offset]
# Line 2541  ostream& escript::operator<<(ostream& o, Line 2578  ostream& escript::operator<<(ostream& o,
2578    return o;    return o;
2579  }  }
2580    
2581    Data
2582    escript::C_GeneralTensorProduct(Data& arg_0,
2583                         Data& arg_1,
2584                         int axis_offset,
2585                         int transpose)
2586    {
2587      // General tensor product: res(SL x SR) = arg_0(SL x SM) * arg_1(SM x SR)
2588      // SM is the product of the last axis_offset entries in arg_0.getShape().
2589    
2590      #if defined DOPROF
2591        // profData->binary++;
2592      #endif
2593    
2594      // Interpolate if necessary and find an appropriate function space
2595      Data arg_0_Z, arg_1_Z;
2596      if (arg_0.getFunctionSpace()!=arg_1.getFunctionSpace()) {
2597        if (arg_0.probeInterpolation(arg_1.getFunctionSpace())) {
2598          arg_0_Z = arg_0.interpolate(arg_1.getFunctionSpace());
2599          arg_1_Z = Data(arg_1);
2600        }
2601        else if (arg_1.probeInterpolation(arg_0.getFunctionSpace())) {
2602          arg_1_Z=arg_1.interpolate(arg_0.getFunctionSpace());
2603          arg_0_Z =Data(arg_0);
2604        }
2605        else {
2606          throw DataException("Error - C_GeneralTensorProduct: arguments have incompatible function spaces.");
2607        }
2608      } else {
2609          arg_0_Z = Data(arg_0);
2610          arg_1_Z = Data(arg_1);
2611      }
2612      // Get rank and shape of inputs
2613      int rank0 = arg_0_Z.getDataPointRank();
2614      int rank1 = arg_1_Z.getDataPointRank();
2615      DataArrayView::ShapeType shape0 = arg_0_Z.getDataPointShape();
2616      DataArrayView::ShapeType shape1 = arg_1_Z.getDataPointShape();
2617    
2618      // Prepare for the loops of the product and verify compatibility of shapes
2619      int start0=0, start1=0;
2620      if (transpose == 0)       {}
2621      else if (transpose == 1)  { start0 = axis_offset; }
2622      else if (transpose == 2)  { start1 = rank1-axis_offset; }
2623      else              { throw DataException("C_GeneralTensorProduct: Error - transpose should be 0, 1 or 2"); }
2624    
2625      // Adjust the shapes for transpose
2626      DataArrayView::ShapeType tmpShape0;
2627      DataArrayView::ShapeType tmpShape1;
2628      for (int i=0; i<rank0; i++)   { tmpShape0.push_back( shape0[(i+start0)%rank0] ); }
2629      for (int i=0; i<rank1; i++)   { tmpShape1.push_back( shape1[(i+start1)%rank1] ); }
2630    
2631    #if 0
2632      // For debugging: show shape after transpose
2633      char tmp[100];
2634      std::string shapeStr;
2635      shapeStr = "(";
2636      for (int i=0; i<rank0; i++)   { sprintf(tmp, "%d,", tmpShape0[i]); shapeStr += tmp; }
2637      shapeStr += ")";
2638      cout << "C_GeneralTensorProduct: Shape of arg0 is " << shapeStr << endl;
2639      shapeStr = "(";
2640      for (int i=0; i<rank1; i++)   { sprintf(tmp, "%d,", tmpShape1[i]); shapeStr += tmp; }
2641      shapeStr += ")";
2642      cout << "C_GeneralTensorProduct: Shape of arg1 is " << shapeStr << endl;
2643    #endif
2644    
2645      // Prepare for the loops of the product
2646      int SL=1, SM=1, SR=1;
2647      for (int i=0; i<rank0-axis_offset; i++)   {
2648        SL *= tmpShape0[i];
2649      }
2650      for (int i=rank0-axis_offset; i<rank0; i++)   {
2651        if (tmpShape0[i] != tmpShape1[i-(rank0-axis_offset)]) {
2652          throw DataException("C_GeneralTensorProduct: Error - incompatible shapes");
2653        }
2654        SM *= tmpShape0[i];
2655      }
2656      for (int i=axis_offset; i<rank1; i++)     {
2657        SR *= tmpShape1[i];
2658      }
2659    
2660      // Define the shape of the output
2661      DataArrayView::ShapeType shape2;
2662      for (int i=0; i<rank0-axis_offset; i++) { shape2.push_back(tmpShape0[i]); } // First part of arg_0_Z
2663      for (int i=axis_offset; i<rank1; i++)   { shape2.push_back(tmpShape1[i]); } // Last part of arg_1_Z
2664    
2665      // Declare output Data object
2666      Data res;
2667    
2668      if      (arg_0_Z.isConstant()   && arg_1_Z.isConstant()) {
2669        res = Data(0.0, shape2, arg_1_Z.getFunctionSpace());    // DataConstant output
2670        double *ptr_0 = &((arg_0_Z.getPointDataView().getData())[0]);
2671        double *ptr_1 = &((arg_1_Z.getPointDataView().getData())[0]);
2672        double *ptr_2 = &((res.getPointDataView().getData())[0]);
2673        matrix_matrix_product(SL, SM, SR, ptr_0, ptr_1, ptr_2, transpose);
2674      }
2675      else if (arg_0_Z.isConstant()   && arg_1_Z.isTagged()) {
2676    
2677        // Prepare the DataConstant input
2678        DataConstant* tmp_0=dynamic_cast<DataConstant*>(arg_0_Z.borrowData());
2679        if (tmp_0==0) { throw DataException("GTP Programming error - casting to DataConstant."); }
2680    
2681        // Borrow DataTagged input from Data object
2682        DataTagged* tmp_1=dynamic_cast<DataTagged*>(arg_1_Z.borrowData());
2683        if (tmp_1==0) { throw DataException("GTP_1 Programming error - casting to DataTagged."); }
2684    
2685        // Prepare a DataTagged output 2
2686        res = Data(0.0, shape2, arg_1_Z.getFunctionSpace());    // DataTagged output
2687        res.tag();
2688        DataTagged* tmp_2=dynamic_cast<DataTagged*>(res.borrowData());
2689        if (tmp_2==0) { throw DataException("GTP Programming error - casting to DataTagged."); }
2690    
2691        // Prepare offset into DataConstant
2692        int offset_0 = tmp_0->getPointOffset(0,0);
2693        double *ptr_0 = &((arg_0_Z.getPointDataView().getData())[offset_0]);
2694        // Get the views
2695        DataArrayView view_1 = tmp_1->getDefaultValue();
2696        DataArrayView view_2 = tmp_2->getDefaultValue();
2697        // Get the pointers to the actual data
2698        double *ptr_1 = &((view_1.getData())[0]);
2699        double *ptr_2 = &((view_2.getData())[0]);
2700        // Compute an MVP for the default
2701        matrix_matrix_product(SL, SM, SR, ptr_0, ptr_1, ptr_2, transpose);
2702        // Compute an MVP for each tag
2703        const DataTagged::DataMapType& lookup_1=tmp_1->getTagLookup();
2704        DataTagged::DataMapType::const_iterator i; // i->first is a tag, i->second is an offset into memory
2705        for (i=lookup_1.begin();i!=lookup_1.end();i++) {
2706          tmp_2->addTaggedValue(i->first,tmp_2->getDefaultValue());
2707          DataArrayView view_1 = tmp_1->getDataPointByTag(i->first);
2708          DataArrayView view_2 = tmp_2->getDataPointByTag(i->first);
2709          double *ptr_1 = &view_1.getData(0);
2710          double *ptr_2 = &view_2.getData(0);
2711          matrix_matrix_product(SL, SM, SR, ptr_0, ptr_1, ptr_2, transpose);
2712        }
2713    
2714      }
2715      else if (arg_0_Z.isConstant()   && arg_1_Z.isExpanded()) {
2716    
2717        res = Data(0.0, shape2, arg_1_Z.getFunctionSpace(),true); // DataExpanded output
2718        DataConstant* tmp_0=dynamic_cast<DataConstant*>(arg_0_Z.borrowData());
2719        DataExpanded* tmp_1=dynamic_cast<DataExpanded*>(arg_1_Z.borrowData());
2720        DataExpanded* tmp_2=dynamic_cast<DataExpanded*>(res.borrowData());
2721        if (tmp_0==0) { throw DataException("GTP Programming error - casting to DataConstant."); }
2722        if (tmp_1==0) { throw DataException("GTP Programming error - casting to DataExpanded."); }
2723        if (tmp_2==0) { throw DataException("GTP Programming error - casting to DataExpanded."); }
2724        int sampleNo_1,dataPointNo_1;
2725        int numSamples_1 = arg_1_Z.getNumSamples();
2726        int numDataPointsPerSample_1 = arg_1_Z.getNumDataPointsPerSample();
2727        int offset_0 = tmp_0->getPointOffset(0,0);
2728        #pragma omp parallel for private(sampleNo_1,dataPointNo_1) schedule(static)
2729        for (sampleNo_1 = 0; sampleNo_1 < numSamples_1; sampleNo_1++) {
2730          for (dataPointNo_1 = 0; dataPointNo_1 < numDataPointsPerSample_1; dataPointNo_1++) {
2731            int offset_1 = tmp_1->getPointOffset(sampleNo_1,dataPointNo_1);
2732            int offset_2 = tmp_2->getPointOffset(sampleNo_1,dataPointNo_1);
2733            double *ptr_0 = &((arg_0_Z.getPointDataView().getData())[offset_0]);
2734            double *ptr_1 = &((arg_1_Z.getPointDataView().getData())[offset_1]);
2735            double *ptr_2 = &((res.getPointDataView().getData())[offset_2]);
2736            matrix_matrix_product(SL, SM, SR, ptr_0, ptr_1, ptr_2, transpose);
2737          }
2738        }
2739    
2740      }
2741      else if (arg_0_Z.isTagged()     && arg_1_Z.isConstant()) {
2742    
2743        // Borrow DataTagged input from Data object
2744        DataTagged* tmp_0=dynamic_cast<DataTagged*>(arg_0_Z.borrowData());
2745        if (tmp_0==0) { throw DataException("GTP_0 Programming error - casting to DataTagged."); }
2746    
2747        // Prepare the DataConstant input
2748        DataConstant* tmp_1=dynamic_cast<DataConstant*>(arg_1_Z.borrowData());
2749        if (tmp_1==0) { throw DataException("GTP Programming error - casting to DataConstant."); }
2750    
2751        // Prepare a DataTagged output 2
2752        res = Data(0.0, shape2, arg_0_Z.getFunctionSpace());    // DataTagged output
2753        res.tag();
2754        DataTagged* tmp_2=dynamic_cast<DataTagged*>(res.borrowData());
2755        if (tmp_2==0) { throw DataException("GTP Programming error - casting to DataTagged."); }
2756    
2757        // Prepare offset into DataConstant
2758        int offset_1 = tmp_1->getPointOffset(0,0);
2759        double *ptr_1 = &((arg_1_Z.getPointDataView().getData())[offset_1]);
2760        // Get the views
2761        DataArrayView view_0 = tmp_0->getDefaultValue();
2762        DataArrayView view_2 = tmp_2->getDefaultValue();
2763        // Get the pointers to the actual data
2764        double *ptr_0 = &((view_0.getData())[0]);
2765        double *ptr_2 = &((view_2.getData())[0]);
2766        // Compute an MVP for the default
2767        matrix_matrix_product(SL, SM, SR, ptr_0, ptr_1, ptr_2, transpose);
2768        // Compute an MVP for each tag
2769        const DataTagged::DataMapType& lookup_0=tmp_0->getTagLookup();
2770        DataTagged::DataMapType::const_iterator i; // i->first is a tag, i->second is an offset into memory
2771        for (i=lookup_0.begin();i!=lookup_0.end();i++) {
2772          tmp_2->addTaggedValue(i->first,tmp_2->getDefaultValue());
2773          DataArrayView view_0 = tmp_0->getDataPointByTag(i->first);
2774          DataArrayView view_2 = tmp_2->getDataPointByTag(i->first);
2775          double *ptr_0 = &view_0.getData(0);
2776          double *ptr_2 = &view_2.getData(0);
2777          matrix_matrix_product(SL, SM, SR, ptr_0, ptr_1, ptr_2, transpose);
2778        }
2779    
2780      }
2781      else if (arg_0_Z.isTagged()     && arg_1_Z.isTagged()) {
2782    
2783        // Borrow DataTagged input from Data object
2784        DataTagged* tmp_0=dynamic_cast<DataTagged*>(arg_0_Z.borrowData());
2785        if (tmp_0==0) { throw DataException("GTP Programming error - casting to DataTagged."); }
2786    
2787        // Borrow DataTagged input from Data object
2788        DataTagged* tmp_1=dynamic_cast<DataTagged*>(arg_1_Z.borrowData());
2789        if (tmp_1==0) { throw DataException("GTP Programming error - casting to DataTagged."); }
2790    
2791        // Prepare a DataTagged output 2
2792        res = Data(0.0, shape2, arg_1_Z.getFunctionSpace());
2793        res.tag();  // DataTagged output
2794        DataTagged* tmp_2=dynamic_cast<DataTagged*>(res.borrowData());
2795        if (tmp_2==0) { throw DataException("GTP Programming error - casting to DataTagged."); }
2796    
2797        // Get the views
2798        DataArrayView view_0 = tmp_0->getDefaultValue();
2799        DataArrayView view_1 = tmp_1->getDefaultValue();
2800        DataArrayView view_2 = tmp_2->getDefaultValue();
2801        // Get the pointers to the actual data
2802        double *ptr_0 = &((view_0.getData())[0]);
2803        double *ptr_1 = &((view_1.getData())[0]);
2804        double *ptr_2 = &((view_2.getData())[0]);
2805        // Compute an MVP for the default
2806        matrix_matrix_product(SL, SM, SR, ptr_0, ptr_1, ptr_2, transpose);
2807        // Merge the tags
2808        DataTagged::DataMapType::const_iterator i; // i->first is a tag, i->second is an offset into memory
2809        const DataTagged::DataMapType& lookup_0=tmp_0->getTagLookup();
2810        const DataTagged::DataMapType& lookup_1=tmp_1->getTagLookup();
2811        for (i=lookup_0.begin();i!=lookup_0.end();i++) {
2812          tmp_2->addTaggedValue(i->first,tmp_2->getDefaultValue()); // use tmp_2 to get correct shape
2813        }
2814        for (i=lookup_1.begin();i!=lookup_1.end();i++) {
2815          tmp_2->addTaggedValue(i->first,tmp_2->getDefaultValue());
2816        }
2817        // Compute an MVP for each tag
2818        const DataTagged::DataMapType& lookup_2=tmp_2->getTagLookup();
2819        for (i=lookup_2.begin();i!=lookup_2.end();i++) {
2820          DataArrayView view_0 = tmp_0->getDataPointByTag(i->first);
2821          DataArrayView view_1 = tmp_1->getDataPointByTag(i->first);
2822          DataArrayView view_2 = tmp_2->getDataPointByTag(i->first);
2823          double *ptr_0 = &view_0.getData(0);
2824          double *ptr_1 = &view_1.getData(0);
2825          double *ptr_2 = &view_2.getData(0);
2826          matrix_matrix_product(SL, SM, SR, ptr_0, ptr_1, ptr_2, transpose);
2827        }
2828    
2829      }
2830      else if (arg_0_Z.isTagged()     && arg_1_Z.isExpanded()) {
2831    
2832        // After finding a common function space above the two inputs have the same numSamples and num DPPS
2833        res = Data(0.0, shape2, arg_1_Z.getFunctionSpace(),true); // DataExpanded output
2834        DataTagged*   tmp_0=dynamic_cast<DataTagged*>(arg_0_Z.borrowData());
2835        DataExpanded* tmp_1=dynamic_cast<DataExpanded*>(arg_1_Z.borrowData());
2836        DataExpanded* tmp_2=dynamic_cast<DataExpanded*>(res.borrowData());
2837        if (tmp_0==0) { throw DataException("GTP Programming error - casting to DataTagged."); }
2838        if (tmp_1==0) { throw DataException("GTP Programming error - casting to DataExpanded."); }
2839        if (tmp_2==0) { throw DataException("GTP Programming error - casting to DataExpanded."); }
2840        int sampleNo_0,dataPointNo_0;
2841        int numSamples_0 = arg_0_Z.getNumSamples();
2842        int numDataPointsPerSample_0 = arg_0_Z.getNumDataPointsPerSample();
2843        #pragma omp parallel for private(sampleNo_0,dataPointNo_0) schedule(static)
2844        for (sampleNo_0 = 0; sampleNo_0 < numSamples_0; sampleNo_0++) {
2845          int offset_0 = tmp_0->getPointOffset(sampleNo_0,0); // They're all the same, so just use #0
2846          double *ptr_0 = &((arg_0_Z.getPointDataView().getData())[offset_0]);
2847          for (dataPointNo_0 = 0; dataPointNo_0 < numDataPointsPerSample_0; dataPointNo_0++) {
2848            int offset_1 = tmp_1->getPointOffset(sampleNo_0,dataPointNo_0);
2849            int offset_2 = tmp_2->getPointOffset(sampleNo_0,dataPointNo_0);
2850            double *ptr_1 = &((arg_1_Z.getPointDataView().getData())[offset_1]);
2851            double *ptr_2 = &((res.getPointDataView().getData())[offset_2]);
2852            matrix_matrix_product(SL, SM, SR, ptr_0, ptr_1, ptr_2, transpose);
2853          }
2854        }
2855    
2856      }
2857      else if (arg_0_Z.isExpanded()   && arg_1_Z.isConstant()) {
2858    
2859        res = Data(0.0, shape2, arg_1_Z.getFunctionSpace(),true); // DataExpanded output
2860        DataExpanded* tmp_0=dynamic_cast<DataExpanded*>(arg_0_Z.borrowData());
2861        DataConstant* tmp_1=dynamic_cast<DataConstant*>(arg_1_Z.borrowData());
2862        DataExpanded* tmp_2=dynamic_cast<DataExpanded*>(res.borrowData());
2863        if (tmp_0==0) { throw DataException("GTP Programming error - casting to DataExpanded."); }
2864        if (tmp_1==0) { throw DataException("GTP Programming error - casting to DataConstant."); }
2865        if (tmp_2==0) { throw DataException("GTP Programming error - casting to DataExpanded."); }
2866        int sampleNo_0,dataPointNo_0;
2867        int numSamples_0 = arg_0_Z.getNumSamples();
2868        int numDataPointsPerSample_0 = arg_0_Z.getNumDataPointsPerSample();
2869        int offset_1 = tmp_1->getPointOffset(0,0);
2870        #pragma omp parallel for private(sampleNo_0,dataPointNo_0) schedule(static)
2871        for (sampleNo_0 = 0; sampleNo_0 < numSamples_0; sampleNo_0++) {
2872          for (dataPointNo_0 = 0; dataPointNo_0 < numDataPointsPerSample_0; dataPointNo_0++) {
2873            int offset_0 = tmp_0->getPointOffset(sampleNo_0,dataPointNo_0);
2874            int offset_2 = tmp_2->getPointOffset(sampleNo_0,dataPointNo_0);
2875            double *ptr_0 = &((arg_0_Z.getPointDataView().getData())[offset_0]);
2876            double *ptr_1 = &((arg_1_Z.getPointDataView().getData())[offset_1]);
2877            double *ptr_2 = &((res.getPointDataView().getData())[offset_2]);
2878            matrix_matrix_product(SL, SM, SR, ptr_0, ptr_1, ptr_2, transpose);
2879          }
2880        }
2881    
2882    
2883      }
2884      else if (arg_0_Z.isExpanded()   && arg_1_Z.isTagged()) {
2885    
2886        // After finding a common function space above the two inputs have the same numSamples and num DPPS
2887        res = Data(0.0, shape2, arg_1_Z.getFunctionSpace(),true); // DataExpanded output
2888        DataExpanded* tmp_0=dynamic_cast<DataExpanded*>(arg_0_Z.borrowData());
2889        DataTagged*   tmp_1=dynamic_cast<DataTagged*>(arg_1_Z.borrowData());
2890        DataExpanded* tmp_2=dynamic_cast<DataExpanded*>(res.borrowData());
2891        if (tmp_0==0) { throw DataException("GTP Programming error - casting to DataExpanded."); }
2892        if (tmp_1==0) { throw DataException("GTP Programming error - casting to DataTagged."); }
2893        if (tmp_2==0) { throw DataException("GTP Programming error - casting to DataExpanded."); }
2894        int sampleNo_0,dataPointNo_0;
2895        int numSamples_0 = arg_0_Z.getNumSamples();
2896        int numDataPointsPerSample_0 = arg_0_Z.getNumDataPointsPerSample();
2897        #pragma omp parallel for private(sampleNo_0,dataPointNo_0) schedule(static)
2898        for (sampleNo_0 = 0; sampleNo_0 < numSamples_0; sampleNo_0++) {
2899          int offset_1 = tmp_1->getPointOffset(sampleNo_0,0);
2900          double *ptr_1 = &((arg_1_Z.getPointDataView().getData())[offset_1]);
2901          for (dataPointNo_0 = 0; dataPointNo_0 < numDataPointsPerSample_0; dataPointNo_0++) {
2902            int offset_0 = tmp_0->getPointOffset(sampleNo_0,dataPointNo_0);
2903            int offset_2 = tmp_2->getPointOffset(sampleNo_0,dataPointNo_0);
2904            double *ptr_0 = &((arg_0_Z.getPointDataView().getData())[offset_0]);
2905            double *ptr_2 = &((res.getPointDataView().getData())[offset_2]);
2906            matrix_matrix_product(SL, SM, SR, ptr_0, ptr_1, ptr_2, transpose);
2907          }
2908        }
2909    
2910      }
2911      else if (arg_0_Z.isExpanded()   && arg_1_Z.isExpanded()) {
2912    
2913        // After finding a common function space above the two inputs have the same numSamples and num DPPS
2914        res = Data(0.0, shape2, arg_1_Z.getFunctionSpace(),true); // DataExpanded output
2915        DataExpanded* tmp_0=dynamic_cast<DataExpanded*>(arg_0_Z.borrowData());
2916        DataExpanded* tmp_1=dynamic_cast<DataExpanded*>(arg_1_Z.borrowData());
2917        DataExpanded* tmp_2=dynamic_cast<DataExpanded*>(res.borrowData());
2918        if (tmp_0==0) { throw DataException("GTP Programming error - casting to DataExpanded."); }
2919        if (tmp_1==0) { throw DataException("GTP Programming error - casting to DataExpanded."); }
2920        if (tmp_2==0) { throw DataException("GTP Programming error - casting to DataExpanded."); }
2921        int sampleNo_0,dataPointNo_0;
2922        int numSamples_0 = arg_0_Z.getNumSamples();
2923        int numDataPointsPerSample_0 = arg_0_Z.getNumDataPointsPerSample();
2924        #pragma omp parallel for private(sampleNo_0,dataPointNo_0) schedule(static)
2925        for (sampleNo_0 = 0; sampleNo_0 < numSamples_0; sampleNo_0++) {
2926          for (dataPointNo_0 = 0; dataPointNo_0 < numDataPointsPerSample_0; dataPointNo_0++) {
2927            int offset_0 = tmp_0->getPointOffset(sampleNo_0,dataPointNo_0);
2928            int offset_1 = tmp_1->getPointOffset(sampleNo_0,dataPointNo_0);
2929            int offset_2 = tmp_2->getPointOffset(sampleNo_0,dataPointNo_0);
2930            double *ptr_0 = &((arg_0_Z.getPointDataView().getData())[offset_0]);
2931            double *ptr_1 = &((arg_1_Z.getPointDataView().getData())[offset_1]);
2932            double *ptr_2 = &((res.getPointDataView().getData())[offset_2]);
2933            matrix_matrix_product(SL, SM, SR, ptr_0, ptr_1, ptr_2, transpose);
2934          }
2935        }
2936    
2937      }
2938      else {
2939        throw DataException("Error - C_GeneralTensorProduct: unknown combination of inputs");
2940      }
2941    
2942      return res;
2943    }
2944    
2945    DataAbstract*
2946    Data::borrowData() const
2947    {
2948      return m_data.get();
2949    }
2950    
2951  /* Member functions specific to the MPI implementation */  /* Member functions specific to the MPI implementation */
2952    
2953  void  void

Legend:
Removed from v.790  
changed lines
  Added in v.826

  ViewVC Help
Powered by ViewVC 1.1.26