/algorithms_for_tapenade/c_code/inv_b.c

  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}