fork download
copy 
  1. #include <cstdlib>
  2. #include <cstring>
  3. #include <iomanip>
  4. #include <iostream>
  5. #include <iterator>
  6.  
  7.  
  8. void errPropag6_orig(double G[21], const double F[6][6], int nF);
  9. void errPropag6_trasat(double G[21], const double F[6][6], int nF);
  10.  
  11. void print(const double (&G)[21], const double (&F)[6][6]);
  12. void print(const double (&G)[21], const double (&G1)[21], const double (&G2)[21]);
  13. double my_rand(const double min, const double max);
  14.  
  15.  
  16.  
  17. int main(int argc, char **argv)
  18. {
  19.    double myG[21] = {};
  20.    double myF[6][6] = {};
  21.  
  22.    srand(1);
  23.  
  24.    // Generate input for the functions
  25.    for (int i = 0; i < 21; i++)
  26.       myG[i] = my_rand(-0.5, 0.5);
  27.  
  28.    for (int i=0; i<6; i++)
  29.       for (int j=0; j<6; j++)
  30.          myF[i][j] = my_rand(-0.5, 0.5);
  31.  
  32.    std::cout << "\n" << "input:\n";
  33.    print(myG, myF);
  34.  
  35.    double myG_orig[21];
  36.    std::copy( std::begin(myG), std::end(myG), std::begin(myG_orig) );
  37.  
  38.    errPropag6_orig(myG_orig, myF, 6);
  39.  
  40.    double myG_trasat[21];
  41.    std::copy( std::begin(myG), std::end(myG), std::begin(myG_trasat) );
  42.  
  43.    errPropag6_trasat(myG_trasat, myF, 6);
  44.  
  45.    std::cout << "\n" << "output:\n";
  46.  
  47.    print(myG, myG_orig, myG_trasat);
  48.  
  49.    return EXIT_SUCCESS;
  50. }
  51.  
  52.  
  53.  
  54. static const int idx66[6][6] =
  55.   {{ 0, 1, 3, 6,10,15},
  56.    { 1, 2, 4, 7,11,16},
  57.    { 3, 4, 5, 8,12,17},
  58.    { 6, 7, 8, 9,13,18},
  59.    {10,11,12,13,14,19},
  60.    {15,16,17,18,19,20}};
  61.  
  62.  
  63.  
  64. void errPropag6_orig(double G[21], const double F[6][6], int nF)
  65. {
  66.   enum {NP=6,NE=21};
  67.  
  68.   double g[NE];      memcpy(g,    G,sizeof( g));
  69.   double fg[NP][NP]; memset(fg[0],0,sizeof(fg));
  70.  
  71.   //double myF[6][6];  memcpy(myF[0],F[0],sizeof(fg));
  72.   //for (int i=0;i<NP;i++) {myF[i][i]+=1.;}
  73.  
  74.   for (int i=0;i<nF;i++) {
  75.   for (int j=0;j<nF;j++) {
  76.     if (!F[i][j])    continue;
  77.     for (int k=0;k<NP;k++) {
  78.       int jk = idx66[j][k];
  79.       if (!g[jk]) continue;
  80.       fg[i][k] += F[i][j]*g[jk];
  81.   }}}
  82.  
  83.   for (int i=0;i<NP;i++) {
  84.   for (int k=i;k<NP;k++) {
  85.     int ik = idx66[i][k];
  86.     double s = 0;
  87.     for (int j=0;j<NP;j++) {
  88.       if (!F[k][j])  continue;
  89.       s += fg[i][j]*F[k][j];
  90.     }
  91.     G[ik] += (s + fg[i][k] + fg[k][i]);
  92.   }}
  93.  
  94. }
  95.  
  96.  
  97.  
  98. inline double* vzero(double *a, int n1)
  99. {
  100.    //to be documented
  101.    if (n1 <= 0) return 0;
  102.    return (double *)memset(a,0,n1*sizeof(double));
  103. }
  104.  
  105.  
  106.  
  107. #define TCL_TRASAT(a, s, r__, m, n) \
  108.    int imax,  k;                    \
  109.    int ia, mn, ir, is, iaa;         \
  110.    double sum;                      \
  111.    --r__;    --s;    --a;           \
  112.    imax = (m * m + m) / 2;          \
  113.    vzero(&r__[1], imax);            \
  114.    mn = m * n;                      \
  115.    int ind = 0;                     \
  116.    int i__ = 0;                     \
  117.    do {                             \
  118.       ind += i__;                   \
  119.       ia = 0; ir = 0;               \
  120.       do {                          \
  121.          is = ind;                  \
  122.          sum = 0.;   k = 0;         \
  123.          do {                       \
  124.             if (k > i__) is += k;   \
  125.             else         ++is;      \
  126.             ++ia;                   \
  127.             sum += s[is] * a[ia];   \
  128.             ++k;                    \
  129.          } while (k < n);           \
  130.          iaa = i__ + 1;             \
  131.          do {                       \
  132.             ++ir;                   \
  133.             r__[ir] += sum * a[iaa];\
  134.             iaa += n;               \
  135.          } while (iaa <= ia);       \
  136.       } while (ia < mn);            \
  137.       ++i__;                        \
  138.    } while (i__ < n);               \
  139.    ++r__;
  140.  
  141.  
  142.  
  143. double* trasat(const double *a, const double *s, double *r__, int m, int n)
  144. {
  145.    TCL_TRASAT(a, s, r__, m, n)
  146.    return r__;
  147. } /* trasat_ */
  148.  
  149.  
  150.  
  151. void errPropag6_trasat(double G[21], const double F[6][6], int nF)
  152. {
  153.   double r[21] = {};
  154.  
  155.   trasat(&F[0][0], G, r, nF, nF);
  156.  
  157.   std::memcpy(G, r, sizeof(double)*21);
  158. }
  159.  
  160.  
  161.  
  162. void print(const double (&G)[21], const double (&F)[6][6])
  163. {
  164.    std::cout << "G: ";
  165.    std::copy(G, G + 21, std::ostream_iterator<double>(std::cout, ", "));
  166.    std::cout << "\n";
  167.  
  168.    std::cout << "F: ";
  169.    std::copy(&F[0][0], &F[0][0] + 36, std::ostream_iterator<double>(std::cout, ", "));
  170.    std::cout << "\n";
  171. }
  172.  
  173.  
  174.  
  175. void print(const double (&G)[21], const double (&G1)[21], const double (&G2)[21])
  176. {
  177.    std::cout << std::fixed << std::setprecision(6) << std::right;
  178.  
  179.    std::cout << std::setw(17) << "index"
  180.              << std::setw(17) << "G in"
  181.              << std::setw(17) << "G out: orig"
  182.              << std::setw(17) << "G out: trasat" << "\n";
  183.  
  184.    for (int i = 0; i < 21; i++)
  185.    {
  186.       std::cout << std::setw(17) << i
  187.                 << std::setw(17) << G[i]
  188.                 << std::setw(17) << G1[i]
  189.                 << std::setw(17) << G2[i] << "\n";
  190.    }
  191. }
  192.  
  193.  
  194.  
  195. double my_rand(const double min, const double max)
  196. {
  197.    return ( std::rand()/(static_cast<double>(RAND_MAX)+1) ) * (max-min) + min;
  198. }
Success #stdin #stdout 0s 16064KB
comments (?)
stdin
copy 
Standard input is empty
stdout
copy 
input:
G: 0.340188, -0.105617, 0.283099, 0.29844, 0.411647, -0.302449, -0.164777, 0.26823, -0.222225, 0.05397, -0.0226029, 0.128871, -0.135216, 0.0134009, 0.45223, 0.416195, 0.135712, 0.217297, -0.358397, 0.106969, -0.483699, 
F: -0.257113, -0.362768, 0.304177, -0.343321, -0.0990556, -0.37021, -0.391191, 0.498925, -0.281743, 0.0129324, 0.339112, 0.11264, -0.203968, 0.137552, 0.0242872, -0.00641701, 0.472775, -0.207483, 0.271358, 0.026745, 0.269914, -0.0997714, 0.391529, -0.216685, -0.147542, 0.307725, 0.419026, -0.430245, 0.449327, 0.0259953, -0.413944, -0.307786, 0.163227, 0.390233, -0.151107, -0.435829, 

output:
            index             G in      G out: orig    G out: trasat
                0         0.340188         0.136731        -0.096171
                1        -0.105617        -0.423876         0.019940
                2         0.283099         0.731080         0.189901
                3         0.298440        -0.316686        -0.133656
                4         0.411647         0.898230         0.235458
                5        -0.302449        -0.415051         0.125760
                6        -0.164777        -0.297413        -0.210375
                7         0.268230         0.683651         0.105343
                8        -0.222225        -0.172629        -0.028366
                9         0.053970        -0.009385        -0.023355
               10        -0.022603        -0.301134        -0.221384
               11         0.128871         0.786228         0.159324
               12        -0.135216         0.112981         0.041210
               13         0.013401         0.167580         0.053599
               14         0.452230         0.971194         0.145870
               15         0.416195         0.254238        -0.061458
               16         0.135712         0.211789         0.205417
               17         0.217297        -0.145527         0.005487
               18        -0.358397        -0.130097        -0.253982
               19         0.106969        -0.056839        -0.262800
               20        -0.483699        -0.569772         0.161517
https://ideone.com/tI1MgY
Test to compare output from errPropag6 and trasat functions
language:
C++14 (gcc 8.3)
created:
7 years ago
visibility:
secret

Share or Embed source code

Discover > Sphere Engine API

The brand new service which powers Ideone!

Discover > IDE Widget

Widget for compiling and running the source code in a web browser!