/[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 825 by ksteube, Tue Aug 29 05:22:50 2006 UTC revision 826 by gross, Tue Aug 29 09:12:33 2006 UTC
# Line 2584  escript::C_GeneralTensorProduct(Data& ar Line 2584  escript::C_GeneralTensorProduct(Data& ar
2584                       int axis_offset,                       int axis_offset,
2585                       int transpose)                       int transpose)
2586  {  {
2587    // General tensor product: arg_2(SL x SR) = arg_0(SL x SM) * arg_1(SM x SR)    // 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().    // SM is the product of the last axis_offset entries in arg_0.getShape().
2589    
2590    #if defined DOPROF    #if defined DOPROF
# Line 2592  escript::C_GeneralTensorProduct(Data& ar Line 2592  escript::C_GeneralTensorProduct(Data& ar
2592    #endif    #endif
2593    
2594    // Interpolate if necessary and find an appropriate function space    // Interpolate if necessary and find an appropriate function space
2595    Data arg_Z;    Data arg_0_Z, arg_1_Z;
2596    if (arg_0.getFunctionSpace()!=arg_1.getFunctionSpace()) {    if (arg_0.getFunctionSpace()!=arg_1.getFunctionSpace()) {
2597      if (arg_0.probeInterpolation(arg_1.getFunctionSpace())) {      if (arg_0.probeInterpolation(arg_1.getFunctionSpace())) {
2598        arg_Z = arg_0.interpolate(arg_1.getFunctionSpace());        arg_0_Z = arg_0.interpolate(arg_1.getFunctionSpace());
2599        // arg_0=arg_Z;        arg_1_Z = Data(arg_1);
2600      }      }
2601      else if (arg_1.probeInterpolation(arg_0.getFunctionSpace())) {      else if (arg_1.probeInterpolation(arg_0.getFunctionSpace())) {
2602        arg_1=arg_1.interpolate(arg_0.getFunctionSpace());        arg_1_Z=arg_1.interpolate(arg_0.getFunctionSpace());
2603          arg_0_Z =Data(arg_0);
2604      }      }
2605      else {      else {
2606        throw DataException("Error - C_GeneralTensorProduct: arguments have incompatible function spaces.");        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    // Get rank and shape of inputs
2613    int rank0 = arg_Z.getDataPointRank();    int rank0 = arg_0_Z.getDataPointRank();
2614    int rank1 = arg_1.getDataPointRank();    int rank1 = arg_1_Z.getDataPointRank();
2615    DataArrayView::ShapeType shape0 = arg_Z.getDataPointShape();    DataArrayView::ShapeType shape0 = arg_0_Z.getDataPointShape();
2616    DataArrayView::ShapeType shape1 = arg_1.getDataPointShape();    DataArrayView::ShapeType shape1 = arg_1_Z.getDataPointShape();
2617    
2618    // Prepare for the loops of the product and verify compatibility of shapes    // Prepare for the loops of the product and verify compatibility of shapes
2619    int start0=0, start1=0;    int start0=0, start1=0;
# Line 2655  escript::C_GeneralTensorProduct(Data& ar Line 2659  escript::C_GeneralTensorProduct(Data& ar
2659    
2660    // Define the shape of the output    // Define the shape of the output
2661    DataArrayView::ShapeType shape2;    DataArrayView::ShapeType shape2;
2662    for (int i=0; i<rank0-axis_offset; i++) { shape2.push_back(tmpShape0[i]); } // First part of arg_Z    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    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    // Declare output Data object
2666    Data arg_2;    Data res;
2667    
2668    if      (arg_Z.isConstant()   && arg_1.isConstant()) {    if      (arg_0_Z.isConstant()   && arg_1_Z.isConstant()) {
2669      arg_2 = Data(0.0, shape2, arg_1.getFunctionSpace());    // DataConstant output      res = Data(0.0, shape2, arg_1_Z.getFunctionSpace());    // DataConstant output
2670      double *ptr_0 = &((arg_Z.getPointDataView().getData())[0]);      double *ptr_0 = &((arg_0_Z.getPointDataView().getData())[0]);
2671      double *ptr_1 = &((arg_1.getPointDataView().getData())[0]);      double *ptr_1 = &((arg_1_Z.getPointDataView().getData())[0]);
2672      double *ptr_2 = &((arg_2.getPointDataView().getData())[0]);      double *ptr_2 = &((res.getPointDataView().getData())[0]);
2673      matrix_matrix_product(SL, SM, SR, ptr_0, ptr_1, ptr_2, transpose);      matrix_matrix_product(SL, SM, SR, ptr_0, ptr_1, ptr_2, transpose);
2674    }    }
2675    else if (arg_Z.isConstant()   && arg_1.isTagged()) {    else if (arg_0_Z.isConstant()   && arg_1_Z.isTagged()) {
2676    
2677      // Prepare the DataConstant input      // Prepare the DataConstant input
2678      DataConstant* tmp_0=dynamic_cast<DataConstant*>(arg_Z.borrowData());      DataConstant* tmp_0=dynamic_cast<DataConstant*>(arg_0_Z.borrowData());
2679      if (tmp_0==0) { throw DataException("GTP Programming error - casting to DataConstant."); }      if (tmp_0==0) { throw DataException("GTP Programming error - casting to DataConstant."); }
2680    
2681      // Borrow DataTagged input from Data object      // Borrow DataTagged input from Data object
2682      DataTagged* tmp_1=dynamic_cast<DataTagged*>(arg_1.borrowData());      DataTagged* tmp_1=dynamic_cast<DataTagged*>(arg_1_Z.borrowData());
2683      if (tmp_1==0) { throw DataException("GTP_1 Programming error - casting to DataTagged."); }      if (tmp_1==0) { throw DataException("GTP_1 Programming error - casting to DataTagged."); }
2684    
2685      // Prepare a DataTagged output 2      // Prepare a DataTagged output 2
2686      arg_2 = Data(0.0, shape2, arg_1.getFunctionSpace());    // DataTagged output      res = Data(0.0, shape2, arg_1_Z.getFunctionSpace());    // DataTagged output
2687      arg_2.tag();      res.tag();
2688      DataTagged* tmp_2=dynamic_cast<DataTagged*>(arg_2.borrowData());      DataTagged* tmp_2=dynamic_cast<DataTagged*>(res.borrowData());
2689      if (tmp_2==0) { throw DataException("GTP Programming error - casting to DataTagged."); }      if (tmp_2==0) { throw DataException("GTP Programming error - casting to DataTagged."); }
2690    
2691      // Prepare offset into DataConstant      // Prepare offset into DataConstant
2692      int offset_0 = tmp_0->getPointOffset(0,0);      int offset_0 = tmp_0->getPointOffset(0,0);
2693      double *ptr_0 = &((arg_Z.getPointDataView().getData())[offset_0]);      double *ptr_0 = &((arg_0_Z.getPointDataView().getData())[offset_0]);
2694      // Get the views      // Get the views
2695      DataArrayView view_1 = tmp_1->getDefaultValue();      DataArrayView view_1 = tmp_1->getDefaultValue();
2696      DataArrayView view_2 = tmp_2->getDefaultValue();      DataArrayView view_2 = tmp_2->getDefaultValue();
# Line 2708  escript::C_GeneralTensorProduct(Data& ar Line 2712  escript::C_GeneralTensorProduct(Data& ar
2712      }      }
2713    
2714    }    }
2715    else if (arg_Z.isConstant()   && arg_1.isExpanded()) {    else if (arg_0_Z.isConstant()   && arg_1_Z.isExpanded()) {
2716    
2717      arg_2 = Data(0.0, shape2, arg_1.getFunctionSpace(),true); // DataExpanded output      res = Data(0.0, shape2, arg_1_Z.getFunctionSpace(),true); // DataExpanded output
2718      DataConstant* tmp_0=dynamic_cast<DataConstant*>(arg_Z.borrowData());      DataConstant* tmp_0=dynamic_cast<DataConstant*>(arg_0_Z.borrowData());
2719      DataExpanded* tmp_1=dynamic_cast<DataExpanded*>(arg_1.borrowData());      DataExpanded* tmp_1=dynamic_cast<DataExpanded*>(arg_1_Z.borrowData());
2720      DataExpanded* tmp_2=dynamic_cast<DataExpanded*>(arg_2.borrowData());      DataExpanded* tmp_2=dynamic_cast<DataExpanded*>(res.borrowData());
2721      if (tmp_0==0) { throw DataException("GTP Programming error - casting to DataConstant."); }      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."); }      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."); }      if (tmp_2==0) { throw DataException("GTP Programming error - casting to DataExpanded."); }
2724      int sampleNo_1,dataPointNo_1;      int sampleNo_1,dataPointNo_1;
2725      int numSamples_1 = arg_1.getNumSamples();      int numSamples_1 = arg_1_Z.getNumSamples();
2726      int numDataPointsPerSample_1 = arg_1.getNumDataPointsPerSample();      int numDataPointsPerSample_1 = arg_1_Z.getNumDataPointsPerSample();
2727      int offset_0 = tmp_0->getPointOffset(0,0);      int offset_0 = tmp_0->getPointOffset(0,0);
2728      #pragma omp parallel for private(sampleNo_1,dataPointNo_1) schedule(static)      #pragma omp parallel for private(sampleNo_1,dataPointNo_1) schedule(static)
2729      for (sampleNo_1 = 0; sampleNo_1 < numSamples_1; sampleNo_1++) {      for (sampleNo_1 = 0; sampleNo_1 < numSamples_1; sampleNo_1++) {
2730        for (dataPointNo_1 = 0; dataPointNo_1 < numDataPointsPerSample_1; dataPointNo_1++) {        for (dataPointNo_1 = 0; dataPointNo_1 < numDataPointsPerSample_1; dataPointNo_1++) {
2731          int offset_1 = tmp_1->getPointOffset(sampleNo_1,dataPointNo_1);          int offset_1 = tmp_1->getPointOffset(sampleNo_1,dataPointNo_1);
2732          int offset_2 = tmp_2->getPointOffset(sampleNo_1,dataPointNo_1);          int offset_2 = tmp_2->getPointOffset(sampleNo_1,dataPointNo_1);
2733          double *ptr_0 = &((arg_Z.getPointDataView().getData())[offset_0]);          double *ptr_0 = &((arg_0_Z.getPointDataView().getData())[offset_0]);
2734          double *ptr_1 = &((arg_1.getPointDataView().getData())[offset_1]);          double *ptr_1 = &((arg_1_Z.getPointDataView().getData())[offset_1]);
2735          double *ptr_2 = &((arg_2.getPointDataView().getData())[offset_2]);          double *ptr_2 = &((res.getPointDataView().getData())[offset_2]);
2736          matrix_matrix_product(SL, SM, SR, ptr_0, ptr_1, ptr_2, transpose);          matrix_matrix_product(SL, SM, SR, ptr_0, ptr_1, ptr_2, transpose);
2737        }        }
2738      }      }
2739    
2740    }    }
2741    else if (arg_Z.isTagged()     && arg_1.isConstant()) {    else if (arg_0_Z.isTagged()     && arg_1_Z.isConstant()) {
2742    
2743      // Borrow DataTagged input from Data object      // Borrow DataTagged input from Data object
2744      DataTagged* tmp_0=dynamic_cast<DataTagged*>(arg_Z.borrowData());      DataTagged* tmp_0=dynamic_cast<DataTagged*>(arg_0_Z.borrowData());
2745      if (tmp_0==0) { throw DataException("GTP_0 Programming error - casting to DataTagged."); }      if (tmp_0==0) { throw DataException("GTP_0 Programming error - casting to DataTagged."); }
2746    
2747      // Prepare the DataConstant input      // Prepare the DataConstant input
2748      DataConstant* tmp_1=dynamic_cast<DataConstant*>(arg_1.borrowData());      DataConstant* tmp_1=dynamic_cast<DataConstant*>(arg_1_Z.borrowData());
2749      if (tmp_1==0) { throw DataException("GTP Programming error - casting to DataConstant."); }      if (tmp_1==0) { throw DataException("GTP Programming error - casting to DataConstant."); }
2750    
2751      // Prepare a DataTagged output 2      // Prepare a DataTagged output 2
2752      arg_2 = Data(0.0, shape2, arg_Z.getFunctionSpace());    // DataTagged output      res = Data(0.0, shape2, arg_0_Z.getFunctionSpace());    // DataTagged output
2753      arg_2.tag();      res.tag();
2754      DataTagged* tmp_2=dynamic_cast<DataTagged*>(arg_2.borrowData());      DataTagged* tmp_2=dynamic_cast<DataTagged*>(res.borrowData());
2755      if (tmp_2==0) { throw DataException("GTP Programming error - casting to DataTagged."); }      if (tmp_2==0) { throw DataException("GTP Programming error - casting to DataTagged."); }
2756    
2757      // Prepare offset into DataConstant      // Prepare offset into DataConstant
2758      int offset_1 = tmp_1->getPointOffset(0,0);      int offset_1 = tmp_1->getPointOffset(0,0);
2759      double *ptr_1 = &((arg_1.getPointDataView().getData())[offset_1]);      double *ptr_1 = &((arg_1_Z.getPointDataView().getData())[offset_1]);
2760      // Get the views      // Get the views
2761      DataArrayView view_0 = tmp_0->getDefaultValue();      DataArrayView view_0 = tmp_0->getDefaultValue();
2762      DataArrayView view_2 = tmp_2->getDefaultValue();      DataArrayView view_2 = tmp_2->getDefaultValue();
# Line 2774  escript::C_GeneralTensorProduct(Data& ar Line 2778  escript::C_GeneralTensorProduct(Data& ar
2778      }      }
2779    
2780    }    }
2781    else if (arg_Z.isTagged()     && arg_1.isTagged()) {    else if (arg_0_Z.isTagged()     && arg_1_Z.isTagged()) {
2782    
2783      // Borrow DataTagged input from Data object      // Borrow DataTagged input from Data object
2784      DataTagged* tmp_0=dynamic_cast<DataTagged*>(arg_Z.borrowData());      DataTagged* tmp_0=dynamic_cast<DataTagged*>(arg_0_Z.borrowData());
2785      if (tmp_0==0) { throw DataException("GTP Programming error - casting to DataTagged."); }      if (tmp_0==0) { throw DataException("GTP Programming error - casting to DataTagged."); }
2786    
2787      // Borrow DataTagged input from Data object      // Borrow DataTagged input from Data object
2788      DataTagged* tmp_1=dynamic_cast<DataTagged*>(arg_1.borrowData());      DataTagged* tmp_1=dynamic_cast<DataTagged*>(arg_1_Z.borrowData());
2789      if (tmp_1==0) { throw DataException("GTP Programming error - casting to DataTagged."); }      if (tmp_1==0) { throw DataException("GTP Programming error - casting to DataTagged."); }
2790    
2791      // Prepare a DataTagged output 2      // Prepare a DataTagged output 2
2792      arg_2 = Data(0.0, shape2, arg_1.getFunctionSpace());      res = Data(0.0, shape2, arg_1_Z.getFunctionSpace());
2793      arg_2.tag();    // DataTagged output      res.tag();  // DataTagged output
2794      DataTagged* tmp_2=dynamic_cast<DataTagged*>(arg_2.borrowData());      DataTagged* tmp_2=dynamic_cast<DataTagged*>(res.borrowData());
2795      if (tmp_2==0) { throw DataException("GTP Programming error - casting to DataTagged."); }      if (tmp_2==0) { throw DataException("GTP Programming error - casting to DataTagged."); }
2796    
2797      // Get the views      // Get the views
# Line 2823  escript::C_GeneralTensorProduct(Data& ar Line 2827  escript::C_GeneralTensorProduct(Data& ar
2827      }      }
2828    
2829    }    }
2830    else if (arg_Z.isTagged()     && arg_1.isExpanded()) {    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      // After finding a common function space above the two inputs have the same numSamples and num DPPS
2833      arg_2 = Data(0.0, shape2, arg_1.getFunctionSpace(),true); // DataExpanded output      res = Data(0.0, shape2, arg_1_Z.getFunctionSpace(),true); // DataExpanded output
2834      DataTagged*   tmp_0=dynamic_cast<DataTagged*>(arg_Z.borrowData());      DataTagged*   tmp_0=dynamic_cast<DataTagged*>(arg_0_Z.borrowData());
2835      DataExpanded* tmp_1=dynamic_cast<DataExpanded*>(arg_1.borrowData());      DataExpanded* tmp_1=dynamic_cast<DataExpanded*>(arg_1_Z.borrowData());
2836      DataExpanded* tmp_2=dynamic_cast<DataExpanded*>(arg_2.borrowData());      DataExpanded* tmp_2=dynamic_cast<DataExpanded*>(res.borrowData());
2837      if (tmp_0==0) { throw DataException("GTP Programming error - casting to DataTagged."); }      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."); }      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."); }      if (tmp_2==0) { throw DataException("GTP Programming error - casting to DataExpanded."); }
2840      int sampleNo_0,dataPointNo_0;      int sampleNo_0,dataPointNo_0;
2841      int numSamples_0 = arg_Z.getNumSamples();      int numSamples_0 = arg_0_Z.getNumSamples();
2842      int numDataPointsPerSample_0 = arg_Z.getNumDataPointsPerSample();      int numDataPointsPerSample_0 = arg_0_Z.getNumDataPointsPerSample();
2843      #pragma omp parallel for private(sampleNo_0,dataPointNo_0) schedule(static)      #pragma omp parallel for private(sampleNo_0,dataPointNo_0) schedule(static)
2844      for (sampleNo_0 = 0; sampleNo_0 < numSamples_0; sampleNo_0++) {      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        int offset_0 = tmp_0->getPointOffset(sampleNo_0,0); // They're all the same, so just use #0
2846        double *ptr_0 = &((arg_Z.getPointDataView().getData())[offset_0]);        double *ptr_0 = &((arg_0_Z.getPointDataView().getData())[offset_0]);
2847        for (dataPointNo_0 = 0; dataPointNo_0 < numDataPointsPerSample_0; dataPointNo_0++) {        for (dataPointNo_0 = 0; dataPointNo_0 < numDataPointsPerSample_0; dataPointNo_0++) {
2848          int offset_1 = tmp_1->getPointOffset(sampleNo_0,dataPointNo_0);          int offset_1 = tmp_1->getPointOffset(sampleNo_0,dataPointNo_0);
2849          int offset_2 = tmp_2->getPointOffset(sampleNo_0,dataPointNo_0);          int offset_2 = tmp_2->getPointOffset(sampleNo_0,dataPointNo_0);
2850          double *ptr_1 = &((arg_1.getPointDataView().getData())[offset_1]);          double *ptr_1 = &((arg_1_Z.getPointDataView().getData())[offset_1]);
2851          double *ptr_2 = &((arg_2.getPointDataView().getData())[offset_2]);          double *ptr_2 = &((res.getPointDataView().getData())[offset_2]);
2852          matrix_matrix_product(SL, SM, SR, ptr_0, ptr_1, ptr_2, transpose);          matrix_matrix_product(SL, SM, SR, ptr_0, ptr_1, ptr_2, transpose);
2853        }        }
2854      }      }
2855    
2856    }    }
2857    else if (arg_Z.isExpanded()   && arg_1.isConstant()) {    else if (arg_0_Z.isExpanded()   && arg_1_Z.isConstant()) {
2858    
2859      arg_2 = Data(0.0, shape2, arg_1.getFunctionSpace(),true); // DataExpanded output      res = Data(0.0, shape2, arg_1_Z.getFunctionSpace(),true); // DataExpanded output
2860      DataExpanded* tmp_0=dynamic_cast<DataExpanded*>(arg_Z.borrowData());      DataExpanded* tmp_0=dynamic_cast<DataExpanded*>(arg_0_Z.borrowData());
2861      DataConstant* tmp_1=dynamic_cast<DataConstant*>(arg_1.borrowData());      DataConstant* tmp_1=dynamic_cast<DataConstant*>(arg_1_Z.borrowData());
2862      DataExpanded* tmp_2=dynamic_cast<DataExpanded*>(arg_2.borrowData());      DataExpanded* tmp_2=dynamic_cast<DataExpanded*>(res.borrowData());
2863      if (tmp_0==0) { throw DataException("GTP Programming error - casting to DataExpanded."); }      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."); }      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."); }      if (tmp_2==0) { throw DataException("GTP Programming error - casting to DataExpanded."); }
2866      int sampleNo_0,dataPointNo_0;      int sampleNo_0,dataPointNo_0;
2867      int numSamples_0 = arg_Z.getNumSamples();      int numSamples_0 = arg_0_Z.getNumSamples();
2868      int numDataPointsPerSample_0 = arg_Z.getNumDataPointsPerSample();      int numDataPointsPerSample_0 = arg_0_Z.getNumDataPointsPerSample();
2869      int offset_1 = tmp_1->getPointOffset(0,0);      int offset_1 = tmp_1->getPointOffset(0,0);
2870      #pragma omp parallel for private(sampleNo_0,dataPointNo_0) schedule(static)      #pragma omp parallel for private(sampleNo_0,dataPointNo_0) schedule(static)
2871      for (sampleNo_0 = 0; sampleNo_0 < numSamples_0; sampleNo_0++) {      for (sampleNo_0 = 0; sampleNo_0 < numSamples_0; sampleNo_0++) {
2872        for (dataPointNo_0 = 0; dataPointNo_0 < numDataPointsPerSample_0; dataPointNo_0++) {        for (dataPointNo_0 = 0; dataPointNo_0 < numDataPointsPerSample_0; dataPointNo_0++) {
2873          int offset_0 = tmp_0->getPointOffset(sampleNo_0,dataPointNo_0);          int offset_0 = tmp_0->getPointOffset(sampleNo_0,dataPointNo_0);
2874          int offset_2 = tmp_2->getPointOffset(sampleNo_0,dataPointNo_0);          int offset_2 = tmp_2->getPointOffset(sampleNo_0,dataPointNo_0);
2875          double *ptr_0 = &((arg_Z.getPointDataView().getData())[offset_0]);          double *ptr_0 = &((arg_0_Z.getPointDataView().getData())[offset_0]);
2876          double *ptr_1 = &((arg_1.getPointDataView().getData())[offset_1]);          double *ptr_1 = &((arg_1_Z.getPointDataView().getData())[offset_1]);
2877          double *ptr_2 = &((arg_2.getPointDataView().getData())[offset_2]);          double *ptr_2 = &((res.getPointDataView().getData())[offset_2]);
2878          matrix_matrix_product(SL, SM, SR, ptr_0, ptr_1, ptr_2, transpose);          matrix_matrix_product(SL, SM, SR, ptr_0, ptr_1, ptr_2, transpose);
2879        }        }
2880      }      }
2881    
2882    
2883    }    }
2884    else if (arg_Z.isExpanded()   && arg_1.isTagged()) {    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      // After finding a common function space above the two inputs have the same numSamples and num DPPS
2887      arg_2 = Data(0.0, shape2, arg_1.getFunctionSpace(),true); // DataExpanded output      res = Data(0.0, shape2, arg_1_Z.getFunctionSpace(),true); // DataExpanded output
2888      DataExpanded* tmp_0=dynamic_cast<DataExpanded*>(arg_Z.borrowData());      DataExpanded* tmp_0=dynamic_cast<DataExpanded*>(arg_0_Z.borrowData());
2889      DataTagged*   tmp_1=dynamic_cast<DataTagged*>(arg_1.borrowData());      DataTagged*   tmp_1=dynamic_cast<DataTagged*>(arg_1_Z.borrowData());
2890      DataExpanded* tmp_2=dynamic_cast<DataExpanded*>(arg_2.borrowData());      DataExpanded* tmp_2=dynamic_cast<DataExpanded*>(res.borrowData());
2891      if (tmp_0==0) { throw DataException("GTP Programming error - casting to DataExpanded."); }      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."); }      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."); }      if (tmp_2==0) { throw DataException("GTP Programming error - casting to DataExpanded."); }
2894      int sampleNo_0,dataPointNo_0;      int sampleNo_0,dataPointNo_0;
2895      int numSamples_0 = arg_Z.getNumSamples();      int numSamples_0 = arg_0_Z.getNumSamples();
2896      int numDataPointsPerSample_0 = arg_Z.getNumDataPointsPerSample();      int numDataPointsPerSample_0 = arg_0_Z.getNumDataPointsPerSample();
2897      #pragma omp parallel for private(sampleNo_0,dataPointNo_0) schedule(static)      #pragma omp parallel for private(sampleNo_0,dataPointNo_0) schedule(static)
2898      for (sampleNo_0 = 0; sampleNo_0 < numSamples_0; sampleNo_0++) {      for (sampleNo_0 = 0; sampleNo_0 < numSamples_0; sampleNo_0++) {
2899        int offset_1 = tmp_1->getPointOffset(sampleNo_0,0);        int offset_1 = tmp_1->getPointOffset(sampleNo_0,0);
2900        double *ptr_1 = &((arg_1.getPointDataView().getData())[offset_1]);        double *ptr_1 = &((arg_1_Z.getPointDataView().getData())[offset_1]);
2901        for (dataPointNo_0 = 0; dataPointNo_0 < numDataPointsPerSample_0; dataPointNo_0++) {        for (dataPointNo_0 = 0; dataPointNo_0 < numDataPointsPerSample_0; dataPointNo_0++) {
2902          int offset_0 = tmp_0->getPointOffset(sampleNo_0,dataPointNo_0);          int offset_0 = tmp_0->getPointOffset(sampleNo_0,dataPointNo_0);
2903          int offset_2 = tmp_2->getPointOffset(sampleNo_0,dataPointNo_0);          int offset_2 = tmp_2->getPointOffset(sampleNo_0,dataPointNo_0);
2904          double *ptr_0 = &((arg_Z.getPointDataView().getData())[offset_0]);          double *ptr_0 = &((arg_0_Z.getPointDataView().getData())[offset_0]);
2905          double *ptr_2 = &((arg_2.getPointDataView().getData())[offset_2]);          double *ptr_2 = &((res.getPointDataView().getData())[offset_2]);
2906          matrix_matrix_product(SL, SM, SR, ptr_0, ptr_1, ptr_2, transpose);          matrix_matrix_product(SL, SM, SR, ptr_0, ptr_1, ptr_2, transpose);
2907        }        }
2908      }      }
2909    
2910    }    }
2911    else if (arg_Z.isExpanded()   && arg_1.isExpanded()) {    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      // After finding a common function space above the two inputs have the same numSamples and num DPPS
2914      arg_2 = Data(0.0, shape2, arg_1.getFunctionSpace(),true); // DataExpanded output      res = Data(0.0, shape2, arg_1_Z.getFunctionSpace(),true); // DataExpanded output
2915      DataExpanded* tmp_0=dynamic_cast<DataExpanded*>(arg_Z.borrowData());      DataExpanded* tmp_0=dynamic_cast<DataExpanded*>(arg_0_Z.borrowData());
2916      DataExpanded* tmp_1=dynamic_cast<DataExpanded*>(arg_1.borrowData());      DataExpanded* tmp_1=dynamic_cast<DataExpanded*>(arg_1_Z.borrowData());
2917      DataExpanded* tmp_2=dynamic_cast<DataExpanded*>(arg_2.borrowData());      DataExpanded* tmp_2=dynamic_cast<DataExpanded*>(res.borrowData());
2918      if (tmp_0==0) { throw DataException("GTP Programming error - casting to DataExpanded."); }      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."); }      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."); }      if (tmp_2==0) { throw DataException("GTP Programming error - casting to DataExpanded."); }
2921      int sampleNo_0,dataPointNo_0;      int sampleNo_0,dataPointNo_0;
2922      int numSamples_0 = arg_Z.getNumSamples();      int numSamples_0 = arg_0_Z.getNumSamples();
2923      int numDataPointsPerSample_0 = arg_Z.getNumDataPointsPerSample();      int numDataPointsPerSample_0 = arg_0_Z.getNumDataPointsPerSample();
2924      #pragma omp parallel for private(sampleNo_0,dataPointNo_0) schedule(static)      #pragma omp parallel for private(sampleNo_0,dataPointNo_0) schedule(static)
2925      for (sampleNo_0 = 0; sampleNo_0 < numSamples_0; sampleNo_0++) {      for (sampleNo_0 = 0; sampleNo_0 < numSamples_0; sampleNo_0++) {
2926        for (dataPointNo_0 = 0; dataPointNo_0 < numDataPointsPerSample_0; dataPointNo_0++) {        for (dataPointNo_0 = 0; dataPointNo_0 < numDataPointsPerSample_0; dataPointNo_0++) {
2927          int offset_0 = tmp_0->getPointOffset(sampleNo_0,dataPointNo_0);          int offset_0 = tmp_0->getPointOffset(sampleNo_0,dataPointNo_0);
2928          int offset_1 = tmp_1->getPointOffset(sampleNo_0,dataPointNo_0);          int offset_1 = tmp_1->getPointOffset(sampleNo_0,dataPointNo_0);
2929          int offset_2 = tmp_2->getPointOffset(sampleNo_0,dataPointNo_0);          int offset_2 = tmp_2->getPointOffset(sampleNo_0,dataPointNo_0);
2930          double *ptr_0 = &((arg_Z.getPointDataView().getData())[offset_0]);          double *ptr_0 = &((arg_0_Z.getPointDataView().getData())[offset_0]);
2931          double *ptr_1 = &((arg_1.getPointDataView().getData())[offset_1]);          double *ptr_1 = &((arg_1_Z.getPointDataView().getData())[offset_1]);
2932          double *ptr_2 = &((arg_2.getPointDataView().getData())[offset_2]);          double *ptr_2 = &((res.getPointDataView().getData())[offset_2]);
2933          matrix_matrix_product(SL, SM, SR, ptr_0, ptr_1, ptr_2, transpose);          matrix_matrix_product(SL, SM, SR, ptr_0, ptr_1, ptr_2, transpose);
2934        }        }
2935      }      }
# Line 2935  escript::C_GeneralTensorProduct(Data& ar Line 2939  escript::C_GeneralTensorProduct(Data& ar
2939      throw DataException("Error - C_GeneralTensorProduct: unknown combination of inputs");      throw DataException("Error - C_GeneralTensorProduct: unknown combination of inputs");
2940    }    }
2941    
2942    return arg_2;    return res;
2943  }  }
2944    
2945  DataAbstract*  DataAbstract*

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

  ViewVC Help
Powered by ViewVC 1.1.26