/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. Tapenade 3.4 (r3375) - 10 Feb 2010 15:08
  3. */
  4. #include "GlobalDeclarations_b.c"
  5. /*
  6. Differentiation of inv in reverse (adjoint) mode:
  7. gradient of useful results: *qt
  8. with respect to varying inputs: *qt *a
  9. RW status of diff variables: *qt:in-out *a:out
  10. */
  11. void inv_b(double *a, double *ab, double *qt, double *qtb, double *r, double *
  12. rb, int na) {
  13. int n, m, k;
  14. double rnn, rnm;
  15. double rnnb, rnmb;
  16. int res;
  17. int res0;
  18. int res1;
  19. double tmp;
  20. int res2;
  21. int res3;
  22. double tmp0;
  23. int res4;
  24. int res5;
  25. int res6;
  26. int res7;
  27. int res8;
  28. double tmp1;
  29. int adFrom;
  30. double tmpb;
  31. double tmp0b;
  32. double tmp1b;
  33. pushreal8_(*r);
  34. pushreal8_(*qt);
  35. qr(a, qt, r, na);
  36. for (n = na-1; n >= 0; --n) {
  37. res = myindex(n, n, na);
  38. pushreal8_(rnn);
  39. rnn = r[res];
  40. for (m = 0; m <= na-1; ++m) {
  41. res0 = myindex(n, m, na);
  42. res1 = myindex(n, m, na);
  43. tmp = r[res1]/rnn;
  44. pushreal8_(r[res0]);
  45. r[res0] = tmp;
  46. res2 = myindex(n, m, na);
  47. res3 = myindex(n, m, na);
  48. tmp0 = qt[res3]/rnn;
  49. pushreal8_(qt[res2]);
  50. qt[res2] = tmp0;
  51. }
  52. adFrom = n + 1;
  53. for (m = adFrom; m <= na-1; ++m) {
  54. res4 = myindex(n, m, na);
  55. pushreal8_(rnm);
  56. rnm = r[res4];
  57. res5 = myindex(n, m, na);
  58. pushreal8_(r[res5]);
  59. r[res5] = 0;
  60. for (k = 0; k <= na-1; ++k) {
  61. res6 = myindex(n, k, na);
  62. res7 = myindex(n, k, na);
  63. res8 = myindex(m, k, na);
  64. tmp1 = qt[res7] - qt[res8]*rnm;
  65. pushreal8_(qt[res6]);
  66. qt[res6] = tmp1;
  67. }
  68. }
  69. pushinteger4_(adFrom);
  70. }
  71. *rb = 0.0;
  72. for (n = 0; n <= na-1; ++n) {
  73. popinteger4_(&adFrom);
  74. for (m = na-1; m >= adFrom; --m) {
  75. rnmb = 0.0;
  76. for (k = na-1; k >= 0; --k) {
  77. popreal8_(&qt[res6]);
  78. tmp1b = qtb[res6];
  79. qtb[res6] = 0.0;
  80. qtb[res7] = qtb[res7] + tmp1b;
  81. qtb[res8] = qtb[res8] - rnm*tmp1b;
  82. rnmb = rnmb - qt[res8]*tmp1b;
  83. }
  84. popreal8_(&r[res5]);
  85. rb[res5] = 0.0;
  86. popreal8_(&rnm);
  87. rb[res4] = rb[res4] + rnmb;
  88. }
  89. rnnb = 0.0;
  90. for (m = na-1; m >= 0; --m) {
  91. popreal8_(&qt[res2]);
  92. tmp0b = qtb[res2];
  93. qtb[res2] = 0.0;
  94. qtb[res3] = qtb[res3] + tmp0b/rnn;
  95. popreal8_(&r[res0]);
  96. tmpb = rb[res0];
  97. rnnb = rnnb - r[res1]*tmpb/(rnn*rnn) - qt[res3]*tmp0b/(rnn*rnn);
  98. rb[res0] = 0.0;
  99. rb[res1] = rb[res1] + tmpb/rnn;
  100. }
  101. popreal8_(&rnn);
  102. rb[res] = rb[res] + rnnb;
  103. }
  104. popreal8_(qt);
  105. popreal8_(r);
  106. qr_b(a, ab, qt, qtb, r, rb, na);
  107. }