/algorithms_for_tapenade/c_code/inv_b.c

http://github.com/b45ch1/hpsc_hanoi_2009_walter · C · 108 lines · 98 code · 1 blank · 9 comment · 8 complexity · dce60ff27967c900523946897f8f3956 MD5 · raw file 

    

  1. /*        Generated by TAPENADE     (INRIA, Tropics team)
    2. 

  2.     Tapenade 3.4 (r3375) - 10 Feb 2010 15:08
    3. 

  3. */
    4. 

  4. #include "GlobalDeclarations_b.c"
    5. 

  5. 

  6. /*
    7. 

  7.   Differentiation of inv in reverse (adjoint) mode:
    8. 

  8.    gradient     of useful results: *qt
    9. 

  9.    with respect to varying inputs: *qt *a
    10. 

  10.    RW status of diff variables: *qt:in-out *a:out
    11. 

  11. */
    12. 

  12. void inv_b(double *a, double *ab, double *qt, double *qtb, double *r, double *
    13. 

  13.         rb, int na) {
    14. 

  14.     int n, m, k;
    15. 

  15.     double rnn, rnm;
    16. 

  16.     double rnnb, rnmb;
    17. 

  17.     int res;
    18. 

  18.     int res0;
    19. 

  19.     int res1;
    20. 

  20.     double tmp;
    21. 

  21.     int res2;
    22. 

  22.     int res3;
    23. 

  23.     double tmp0;
    24. 

  24.     int res4;
    25. 

  25.     int res5;
    26. 

  26.     int res6;
    27. 

  27.     int res7;
    28. 

  28.     int res8;
    29. 

  29.     double tmp1;
    30. 

  30.     int adFrom;
    31. 

  31.     double tmpb;
    32. 

  32.     double tmp0b;
    33. 

  33.     double tmp1b;
    34. 

  34.     pushreal8_(*r);
    35. 

  35.     pushreal8_(*qt);
    36. 

  36.     qr(a, qt, r, na);
    37. 

  37.     for (n = na-1; n >= 0; --n) {
    38. 

  38.         res = myindex(n, n, na);
    39. 

  39.         pushreal8_(rnn);
    40. 

  40.         rnn = r[res];
    41. 

  41.         for (m = 0; m <= na-1; ++m) {
    42. 

  42.             res0 = myindex(n, m, na);
    43. 

  43.             res1 = myindex(n, m, na);
    44. 

  44.             tmp = r[res1]/rnn;
    45. 

  45.             pushreal8_(r[res0]);
    46. 

  46.             r[res0] = tmp;
    47. 

  47.             res2 = myindex(n, m, na);
    48. 

  48.             res3 = myindex(n, m, na);
    49. 

  49.             tmp0 = qt[res3]/rnn;
    50. 

  50.             pushreal8_(qt[res2]);
    51. 

  51.             qt[res2] = tmp0;
    52. 

  52.         }
    53. 

  53.         adFrom = n + 1;
    54. 

  54.         for (m = adFrom; m <= na-1; ++m) {
    55. 

  55.             res4 = myindex(n, m, na);
    56. 

  56.             pushreal8_(rnm);
    57. 

  57.             rnm = r[res4];
    58. 

  58.             res5 = myindex(n, m, na);
    59. 

  59.             pushreal8_(r[res5]);
    60. 

  60.             r[res5] = 0;
    61. 

  61.             for (k = 0; k <= na-1; ++k) {
    62. 

  62.                 res6 = myindex(n, k, na);
    63. 

  63.                 res7 = myindex(n, k, na);
    64. 

  64.                 res8 = myindex(m, k, na);
    65. 

  65.                 tmp1 = qt[res7] - qt[res8]*rnm;
    66. 

  66.                 pushreal8_(qt[res6]);
    67. 

  67.                 qt[res6] = tmp1;
    68. 

  68.             }
    69. 

  69.         }
    70. 

  70.         pushinteger4_(adFrom);
    71. 

  71.     }
    72. 

  72.     *rb = 0.0;
    73. 

  73.     for (n = 0; n <= na-1; ++n) {
    74. 

  74.         popinteger4_(&adFrom);
    75. 

  75.         for (m = na-1; m >= adFrom; --m) {
    76. 

  76.             rnmb = 0.0;
    77. 

  77.             for (k = na-1; k >= 0; --k) {
    78. 

  78.                 popreal8_(&qt[res6]);
    79. 

  79.                 tmp1b = qtb[res6];
    80. 

  80.                 qtb[res6] = 0.0;
    81. 

  81.                 qtb[res7] = qtb[res7] + tmp1b;
    82. 

  82.                 qtb[res8] = qtb[res8] - rnm*tmp1b;
    83. 

  83.                 rnmb = rnmb - qt[res8]*tmp1b;
    84. 

  84.             }
    85. 

  85.             popreal8_(&r[res5]);
    86. 

  86.             rb[res5] = 0.0;
    87. 

  87.             popreal8_(&rnm);
    88. 

  88.             rb[res4] = rb[res4] + rnmb;
    89. 

  89.         }
    90. 

  90.         rnnb = 0.0;
    91. 

  91.         for (m = na-1; m >= 0; --m) {
    92. 

  92.             popreal8_(&qt[res2]);
    93. 

  93.             tmp0b = qtb[res2];
    94. 

  94.             qtb[res2] = 0.0;
    95. 

  95.             qtb[res3] = qtb[res3] + tmp0b/rnn;
    96. 

  96.             popreal8_(&r[res0]);
    97. 

  97.             tmpb = rb[res0];
    98. 

  98.             rnnb = rnnb - r[res1]*tmpb/(rnn*rnn) - qt[res3]*tmp0b/(rnn*rnn);
    99. 

  99.             rb[res0] = 0.0;
    100. 

  100.             rb[res1] = rb[res1] + tmpb/rnn;
    101. 

  101.         }
    102. 

  102.         popreal8_(&rnn);
    103. 

  103.         rb[res] = rb[res] + rnnb;
    104. 

  104.     }
    105. 

  105.     popreal8_(qt);
    106. 

  106.     popreal8_(r);
    107. 

  107.     qr_b(a, ab, qt, qtb, r, rb, na);
    108. 

  108. }
    109.