/[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 812 by gross, Thu Aug 10 01:12:16 2006 UTC revision 813 by ksteube, Mon Aug 21 02:08:47 2006 UTC
# Line 2578  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: arg_2(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      if (arg_0.getFunctionSpace()!=arg_1.getFunctionSpace()) {
2596        if (arg_0.probeInterpolation(arg_1.getFunctionSpace())) {
2597          Data arg_0_tmp = arg_0.interpolate(arg_1.getFunctionSpace());
2598          arg_1=arg_0_tmp;
2599        }
2600        else if (arg_1.probeInterpolation(arg_0.getFunctionSpace())) {
2601          arg_1=arg_1.interpolate(arg_0.getFunctionSpace());
2602        }
2603        else {
2604          throw DataException("Error - C_GeneralTensorProduct: arguments have incompatible function spaces.");
2605        }
2606      }
2607    
2608      // Get rank and shape of inputs
2609      int rank0 = arg_0.getDataPointRank();
2610      int rank1 = arg_1.getDataPointRank();
2611      DataArrayView::ShapeType shape0 = arg_0.getDataPointShape();
2612      DataArrayView::ShapeType shape1 = arg_1.getDataPointShape();
2613    
2614      // Prepare for the loops of the product and verify compatibility of shapes
2615      int start0=0, start1=0;
2616      if (transpose == 0)       {}
2617      else if (transpose == 1)  { start0 = axis_offset; }
2618      else if (transpose == 2)  { start1 = rank1-axis_offset; }
2619      else              { throw DataException("C_GeneralTensorProduct: Error - transpose should be 0, 1 or 2"); }
2620    
2621      // Adjust the shapes for transpose
2622      DataArrayView::ShapeType tmpShape0;
2623      DataArrayView::ShapeType tmpShape1;
2624      for (int i=0; i<rank0; i++)   { tmpShape0.push_back( shape0[(i+start0)%rank0] ); }
2625      for (int i=0; i<rank1; i++)   { tmpShape1.push_back( shape1[(i+start1)%rank1] ); }
2626    
2627    #if 0
2628      // For debugging: show shape after transpose
2629      char tmp[100];
2630      std::string shapeStr;
2631      shapeStr = "(";
2632      for (int i=0; i<rank0; i++)   { sprintf(tmp, "%d,", tmpShape0[i]); shapeStr += tmp; }
2633      shapeStr += ")";
2634      cout << "C_GeneralTensorProduct: Shape of arg0 is " << shapeStr << endl;
2635      shapeStr = "(";
2636      for (int i=0; i<rank1; i++)   { sprintf(tmp, "%d,", tmpShape1[i]); shapeStr += tmp; }
2637      shapeStr += ")";
2638      cout << "C_GeneralTensorProduct: Shape of arg1 is " << shapeStr << endl;
2639    #endif
2640    
2641      // Prepare for the loops of the product
2642      int SL=1, SM=1, SR=1;
2643      for (int i=0; i<rank0-axis_offset; i++)   {
2644        SL *= tmpShape0[i];
2645      }
2646      for (int i=rank0-axis_offset; i<rank0; i++)   {
2647        if (tmpShape0[i] != tmpShape1[i-(rank0-axis_offset)]) {
2648          throw DataException("C_GeneralTensorProduct: Error - incompatible shapes");
2649        }
2650        SM *= tmpShape0[i];
2651      }
2652      for (int i=axis_offset; i<rank1; i++)     {
2653        SR *= tmpShape1[i];
2654      }
2655    
2656      // Define the shape of the output
2657      DataArrayView::ShapeType shape2;
2658      for (int i=0; i<rank0-axis_offset; i++) { shape2.push_back(tmpShape0[i]); } // First part of arg_0
2659      for (int i=axis_offset; i<rank1; i++)   { shape2.push_back(tmpShape1[i]); } // Last part of arg_1
2660    
2661      // Declare output Data object
2662      Data arg_2;
2663    
2664      if      (arg_0.isConstant()   && arg_1.isConstant()) {
2665        arg_2 = Data(0.0, shape2, arg_1.getFunctionSpace());    // DataConstant output
2666        double *ptr_0 = &((arg_0.getPointDataView().getData())[0]);
2667        double *ptr_1 = &((arg_1.getPointDataView().getData())[0]);
2668        double *ptr_2 = &((arg_2.getPointDataView().getData())[0]);
2669        matrix_matrix_product(SL, SM, SR, ptr_0, ptr_1, ptr_2, transpose);
2670      }
2671      else if (arg_0.isConstant()   && arg_1.isTagged()) {
2672    
2673        // Prepare the DataConstant input
2674        DataConstant* tmp_0=dynamic_cast<DataConstant*>(arg_0.borrowData());
2675        if (tmp_0==0) { throw DataException("GTP Programming error - casting to DataConstant."); }
2676    
2677        // Borrow DataTagged input from Data object
2678        DataTagged* tmp_1=dynamic_cast<DataTagged*>(arg_1.borrowData());
2679        if (tmp_1==0) { throw DataException("GTP_1 Programming error - casting to DataTagged."); }
2680    
2681        // Prepare a DataTagged output 2
2682        arg_2 = Data(0.0, shape2, arg_1.getFunctionSpace());    // DataTagged output
2683        arg_2.tag();
2684        DataTagged* tmp_2=dynamic_cast<DataTagged*>(arg_2.borrowData());
2685        if (tmp_2==0) { throw DataException("GTP Programming error - casting to DataTagged."); }
2686    
2687        // Prepare offset into DataConstant
2688        int offset_0 = tmp_0->getPointOffset(0,0);
2689        double *ptr_0 = &((arg_0.getPointDataView().getData())[offset_0]);
2690        // Get the views
2691        DataArrayView view_1 = tmp_1->getDefaultValue();
2692        DataArrayView view_2 = tmp_2->getDefaultValue();
2693        // Get the pointers to the actual data
2694        double *ptr_1 = &((view_1.getData())[0]);
2695        double *ptr_2 = &((view_2.getData())[0]);
2696        // Compute an MVP for the default
2697        matrix_matrix_product(SL, SM, SR, ptr_0, ptr_1, ptr_2, transpose);
2698        // Compute an MVP for each tag
2699        const DataTagged::DataMapType& lookup_1=tmp_1->getTagLookup();
2700        DataTagged::DataMapType::const_iterator i; // i->first is a tag, i->second is an offset into memory
2701        for (i=lookup_1.begin();i!=lookup_1.end();i++) {
2702          tmp_2->addTaggedValue(i->first,tmp_2->getDefaultValue());
2703          DataArrayView view_1 = tmp_1->getDataPointByTag(i->first);
2704          DataArrayView view_2 = tmp_2->getDataPointByTag(i->first);
2705          double *ptr_1 = &view_1.getData(0);
2706          double *ptr_2 = &view_2.getData(0);
2707          matrix_matrix_product(SL, SM, SR, ptr_0, ptr_1, ptr_2, transpose);
2708        }
2709    
2710      }
2711      else if (arg_0.isConstant()   && arg_1.isExpanded()) {
2712    
2713        arg_2 = Data(0.0, shape2, arg_1.getFunctionSpace(),true); // DataExpanded output
2714        DataConstant* tmp_0=dynamic_cast<DataConstant*>(arg_0.borrowData());
2715        DataExpanded* tmp_1=dynamic_cast<DataExpanded*>(arg_1.borrowData());
2716        DataExpanded* tmp_2=dynamic_cast<DataExpanded*>(arg_2.borrowData());
2717        if (tmp_0==0) { throw DataException("GTP Programming error - casting to DataConstant."); }
2718        if (tmp_1==0) { throw DataException("GTP Programming error - casting to DataExpanded."); }
2719        if (tmp_2==0) { throw DataException("GTP Programming error - casting to DataExpanded."); }
2720        int sampleNo_1,dataPointNo_1;
2721        int numSamples_1 = arg_1.getNumSamples();
2722        int numDataPointsPerSample_1 = arg_1.getNumDataPointsPerSample();
2723        int offset_0 = tmp_0->getPointOffset(0,0);
2724        #pragma omp parallel for private(sampleNo_1,dataPointNo_1) schedule(static)
2725        for (sampleNo_1 = 0; sampleNo_1 < numSamples_1; sampleNo_1++) {
2726          for (dataPointNo_1 = 0; dataPointNo_1 < numDataPointsPerSample_1; dataPointNo_1++) {
2727            int offset_1 = tmp_1->getPointOffset(sampleNo_1,dataPointNo_1);
2728            int offset_2 = tmp_2->getPointOffset(sampleNo_1,dataPointNo_1);
2729            double *ptr_0 = &((arg_0.getPointDataView().getData())[offset_0]);
2730            double *ptr_1 = &((arg_1.getPointDataView().getData())[offset_1]);
2731            double *ptr_2 = &((arg_2.getPointDataView().getData())[offset_2]);
2732            matrix_matrix_product(SL, SM, SR, ptr_0, ptr_1, ptr_2, transpose);
2733          }
2734        }
2735    
2736      }
2737      else if (arg_0.isTagged()     && arg_1.isConstant()) {
2738    
2739        // Borrow DataTagged input from Data object
2740        DataTagged* tmp_0=dynamic_cast<DataTagged*>(arg_0.borrowData());
2741        if (tmp_0==0) { throw DataException("GTP_0 Programming error - casting to DataTagged."); }
2742    
2743        // Prepare the DataConstant input
2744        DataConstant* tmp_1=dynamic_cast<DataConstant*>(arg_1.borrowData());
2745        if (tmp_1==0) { throw DataException("GTP Programming error - casting to DataConstant."); }
2746    
2747        // Prepare a DataTagged output 2
2748        arg_2 = Data(0.0, shape2, arg_0.getFunctionSpace());    // DataTagged output
2749        arg_2.tag();
2750        DataTagged* tmp_2=dynamic_cast<DataTagged*>(arg_2.borrowData());
2751        if (tmp_2==0) { throw DataException("GTP Programming error - casting to DataTagged."); }
2752    
2753        // Prepare offset into DataConstant
2754        int offset_1 = tmp_1->getPointOffset(0,0);
2755        double *ptr_1 = &((arg_1.getPointDataView().getData())[offset_1]);
2756        // Get the views
2757        DataArrayView view_0 = tmp_0->getDefaultValue();
2758        DataArrayView view_2 = tmp_2->getDefaultValue();
2759        // Get the pointers to the actual data
2760        double *ptr_0 = &((view_0.getData())[0]);
2761        double *ptr_2 = &((view_2.getData())[0]);
2762        // Compute an MVP for the default
2763        matrix_matrix_product(SL, SM, SR, ptr_0, ptr_1, ptr_2, transpose);
2764        // Compute an MVP for each tag
2765        const DataTagged::DataMapType& lookup_0=tmp_0->getTagLookup();
2766        DataTagged::DataMapType::const_iterator i; // i->first is a tag, i->second is an offset into memory
2767        for (i=lookup_0.begin();i!=lookup_0.end();i++) {
2768          tmp_2->addTaggedValue(i->first,tmp_2->getDefaultValue());
2769          DataArrayView view_0 = tmp_0->getDataPointByTag(i->first);
2770          DataArrayView view_2 = tmp_2->getDataPointByTag(i->first);
2771          double *ptr_0 = &view_0.getData(0);
2772          double *ptr_2 = &view_2.getData(0);
2773          matrix_matrix_product(SL, SM, SR, ptr_0, ptr_1, ptr_2, transpose);
2774        }
2775    
2776      }
2777      else if (arg_0.isTagged()     && arg_1.isTagged()) {
2778    
2779        // Borrow DataTagged input from Data object
2780        DataTagged* tmp_0=dynamic_cast<DataTagged*>(arg_0.borrowData());
2781        if (tmp_0==0) { throw DataException("GTP Programming error - casting to DataTagged."); }
2782    
2783        // Borrow DataTagged input from Data object
2784        DataTagged* tmp_1=dynamic_cast<DataTagged*>(arg_1.borrowData());
2785        if (tmp_1==0) { throw DataException("GTP Programming error - casting to DataTagged."); }
2786    
2787        // Prepare a DataTagged output 2
2788        arg_2 = Data(0.0, shape2, arg_1.getFunctionSpace());
2789        arg_2.tag();    // DataTagged output
2790        DataTagged* tmp_2=dynamic_cast<DataTagged*>(arg_2.borrowData());
2791        if (tmp_2==0) { throw DataException("GTP Programming error - casting to DataTagged."); }
2792    
2793        // Get the views
2794        DataArrayView view_0 = tmp_0->getDefaultValue();
2795        DataArrayView view_1 = tmp_1->getDefaultValue();
2796        DataArrayView view_2 = tmp_2->getDefaultValue();
2797        // Get the pointers to the actual data
2798        double *ptr_0 = &((view_0.getData())[0]);
2799        double *ptr_1 = &((view_1.getData())[0]);
2800        double *ptr_2 = &((view_2.getData())[0]);
2801        // Compute an MVP for the default
2802        matrix_matrix_product(SL, SM, SR, ptr_0, ptr_1, ptr_2, transpose);
2803        // Merge the tags
2804        DataTagged::DataMapType::const_iterator i; // i->first is a tag, i->second is an offset into memory
2805        const DataTagged::DataMapType& lookup_0=tmp_0->getTagLookup();
2806        const DataTagged::DataMapType& lookup_1=tmp_1->getTagLookup();
2807        for (i=lookup_0.begin();i!=lookup_0.end();i++) {
2808          tmp_2->addTaggedValue(i->first,tmp_2->getDefaultValue()); // use tmp_2 to get correct shape
2809        }
2810        for (i=lookup_1.begin();i!=lookup_1.end();i++) {
2811          tmp_2->addTaggedValue(i->first,tmp_2->getDefaultValue());
2812        }
2813        // Compute an MVP for each tag
2814        const DataTagged::DataMapType& lookup_2=tmp_2->getTagLookup();
2815        for (i=lookup_2.begin();i!=lookup_2.end();i++) {
2816          DataArrayView view_0 = tmp_0->getDataPointByTag(i->first);
2817          DataArrayView view_1 = tmp_1->getDataPointByTag(i->first);
2818          DataArrayView view_2 = tmp_2->getDataPointByTag(i->first);
2819          double *ptr_0 = &view_0.getData(0);
2820          double *ptr_1 = &view_1.getData(0);
2821          double *ptr_2 = &view_2.getData(0);
2822          matrix_matrix_product(SL, SM, SR, ptr_0, ptr_1, ptr_2, transpose);
2823        }
2824    
2825      }
2826      else if (arg_0.isTagged()     && arg_1.isExpanded()) {
2827    
2828        // After finding a common function space above the two inputs have the same numSamples and num DPPS
2829        arg_2 = Data(0.0, shape2, arg_1.getFunctionSpace(),true); // DataExpanded output
2830        DataTagged*   tmp_0=dynamic_cast<DataTagged*>(arg_0.borrowData());
2831        DataExpanded* tmp_1=dynamic_cast<DataExpanded*>(arg_1.borrowData());
2832        DataExpanded* tmp_2=dynamic_cast<DataExpanded*>(arg_2.borrowData());
2833        if (tmp_0==0) { throw DataException("GTP Programming error - casting to DataTagged."); }
2834        if (tmp_1==0) { throw DataException("GTP Programming error - casting to DataExpanded."); }
2835        if (tmp_2==0) { throw DataException("GTP Programming error - casting to DataExpanded."); }
2836        int sampleNo_0,dataPointNo_0;
2837        int numSamples_0 = arg_0.getNumSamples();
2838        int numDataPointsPerSample_0 = arg_0.getNumDataPointsPerSample();
2839        #pragma omp parallel for private(sampleNo_0,dataPointNo_0) schedule(static)
2840        for (sampleNo_0 = 0; sampleNo_0 < numSamples_0; sampleNo_0++) {
2841          int offset_0 = tmp_0->getPointOffset(sampleNo_0,0); // They're all the same, so just use #0
2842          double *ptr_0 = &((arg_0.getPointDataView().getData())[offset_0]);
2843          for (dataPointNo_0 = 0; dataPointNo_0 < numDataPointsPerSample_0; dataPointNo_0++) {
2844            int offset_1 = tmp_1->getPointOffset(sampleNo_0,dataPointNo_0);
2845            int offset_2 = tmp_2->getPointOffset(sampleNo_0,dataPointNo_0);
2846            double *ptr_1 = &((arg_1.getPointDataView().getData())[offset_1]);
2847            double *ptr_2 = &((arg_2.getPointDataView().getData())[offset_2]);
2848            matrix_matrix_product(SL, SM, SR, ptr_0, ptr_1, ptr_2, transpose);
2849          }
2850        }
2851    
2852      }
2853      else if (arg_0.isExpanded()   && arg_1.isConstant()) {
2854    
2855        arg_2 = Data(0.0, shape2, arg_1.getFunctionSpace(),true); // DataExpanded output
2856        DataExpanded* tmp_0=dynamic_cast<DataExpanded*>(arg_0.borrowData());
2857        DataConstant* tmp_1=dynamic_cast<DataConstant*>(arg_1.borrowData());
2858        DataExpanded* tmp_2=dynamic_cast<DataExpanded*>(arg_2.borrowData());
2859        if (tmp_0==0) { throw DataException("GTP Programming error - casting to DataExpanded."); }
2860        if (tmp_1==0) { throw DataException("GTP Programming error - casting to DataConstant."); }
2861        if (tmp_2==0) { throw DataException("GTP Programming error - casting to DataExpanded."); }
2862        int sampleNo_0,dataPointNo_0;
2863        int numSamples_0 = arg_0.getNumSamples();
2864        int numDataPointsPerSample_0 = arg_0.getNumDataPointsPerSample();
2865        int offset_1 = tmp_1->getPointOffset(0,0);
2866        #pragma omp parallel for private(sampleNo_0,dataPointNo_0) schedule(static)
2867        for (sampleNo_0 = 0; sampleNo_0 < numSamples_0; sampleNo_0++) {
2868          for (dataPointNo_0 = 0; dataPointNo_0 < numDataPointsPerSample_0; dataPointNo_0++) {
2869            int offset_0 = tmp_0->getPointOffset(sampleNo_0,dataPointNo_0);
2870            int offset_2 = tmp_2->getPointOffset(sampleNo_0,dataPointNo_0);
2871            double *ptr_0 = &((arg_0.getPointDataView().getData())[offset_0]);
2872            double *ptr_1 = &((arg_1.getPointDataView().getData())[offset_1]);
2873            double *ptr_2 = &((arg_2.getPointDataView().getData())[offset_2]);
2874            matrix_matrix_product(SL, SM, SR, ptr_0, ptr_1, ptr_2, transpose);
2875          }
2876        }
2877    
2878    
2879      }
2880      else if (arg_0.isExpanded()   && arg_1.isTagged()) {
2881    
2882        // After finding a common function space above the two inputs have the same numSamples and num DPPS
2883        arg_2 = Data(0.0, shape2, arg_1.getFunctionSpace(),true); // DataExpanded output
2884        DataExpanded* tmp_0=dynamic_cast<DataExpanded*>(arg_0.borrowData());
2885        DataTagged*   tmp_1=dynamic_cast<DataTagged*>(arg_1.borrowData());
2886        DataExpanded* tmp_2=dynamic_cast<DataExpanded*>(arg_2.borrowData());
2887        if (tmp_0==0) { throw DataException("GTP Programming error - casting to DataExpanded."); }
2888        if (tmp_1==0) { throw DataException("GTP Programming error - casting to DataTagged."); }
2889        if (tmp_2==0) { throw DataException("GTP Programming error - casting to DataExpanded."); }
2890        int sampleNo_0,dataPointNo_0;
2891        int numSamples_0 = arg_0.getNumSamples();
2892        int numDataPointsPerSample_0 = arg_0.getNumDataPointsPerSample();
2893        #pragma omp parallel for private(sampleNo_0,dataPointNo_0) schedule(static)
2894        for (sampleNo_0 = 0; sampleNo_0 < numSamples_0; sampleNo_0++) {
2895          int offset_1 = tmp_1->getPointOffset(sampleNo_0,0);
2896          double *ptr_1 = &((arg_1.getPointDataView().getData())[offset_1]);
2897          for (dataPointNo_0 = 0; dataPointNo_0 < numDataPointsPerSample_0; dataPointNo_0++) {
2898            int offset_0 = tmp_0->getPointOffset(sampleNo_0,dataPointNo_0);
2899            int offset_2 = tmp_2->getPointOffset(sampleNo_0,dataPointNo_0);
2900            double *ptr_0 = &((arg_0.getPointDataView().getData())[offset_0]);
2901            double *ptr_2 = &((arg_2.getPointDataView().getData())[offset_2]);
2902            matrix_matrix_product(SL, SM, SR, ptr_0, ptr_1, ptr_2, transpose);
2903          }
2904        }
2905    
2906      }
2907      else if (arg_0.isExpanded()   && arg_1.isExpanded()) {
2908    
2909        // After finding a common function space above the two inputs have the same numSamples and num DPPS
2910        arg_2 = Data(0.0, shape2, arg_1.getFunctionSpace(),true); // DataExpanded output
2911        DataExpanded* tmp_0=dynamic_cast<DataExpanded*>(arg_0.borrowData());
2912        DataExpanded* tmp_1=dynamic_cast<DataExpanded*>(arg_1.borrowData());
2913        DataExpanded* tmp_2=dynamic_cast<DataExpanded*>(arg_2.borrowData());
2914        if (tmp_0==0) { throw DataException("GTP Programming error - casting to DataExpanded."); }
2915        if (tmp_1==0) { throw DataException("GTP Programming error - casting to DataExpanded."); }
2916        if (tmp_2==0) { throw DataException("GTP Programming error - casting to DataExpanded."); }
2917        int sampleNo_0,dataPointNo_0;
2918        int numSamples_0 = arg_0.getNumSamples();
2919        int numDataPointsPerSample_0 = arg_0.getNumDataPointsPerSample();
2920        #pragma omp parallel for private(sampleNo_0,dataPointNo_0) schedule(static)
2921        for (sampleNo_0 = 0; sampleNo_0 < numSamples_0; sampleNo_0++) {
2922          for (dataPointNo_0 = 0; dataPointNo_0 < numDataPointsPerSample_0; dataPointNo_0++) {
2923            int offset_0 = tmp_0->getPointOffset(sampleNo_0,dataPointNo_0);
2924            int offset_1 = tmp_1->getPointOffset(sampleNo_0,dataPointNo_0);
2925            int offset_2 = tmp_2->getPointOffset(sampleNo_0,dataPointNo_0);
2926            double *ptr_0 = &((arg_0.getPointDataView().getData())[offset_0]);
2927            double *ptr_1 = &((arg_1.getPointDataView().getData())[offset_1]);
2928            double *ptr_2 = &((arg_2.getPointDataView().getData())[offset_2]);
2929            matrix_matrix_product(SL, SM, SR, ptr_0, ptr_1, ptr_2, transpose);
2930          }
2931        }
2932    
2933      }
2934      else {
2935        throw DataException("Error - C_GeneralTensorProduct: unknown combination of inputs");
2936      }
2937    
2938      return arg_2;
2939    }
2940    
2941    DataAbstract*
2942    Data::borrowData() const
2943    {
2944      return m_data.get();
2945    }
2946    
2947  /* Member functions specific to the MPI implementation */  /* Member functions specific to the MPI implementation */
2948    
2949  void  void
Legend:
Removed from v.812  
changed lines
  Added in v.813
  ViewVC Help
Powered by ViewVC 1.1.26