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