/algorithms_for_tapenade/fortran_code/inv_b.f

 1C        Generated by TAPENADE     (INRIA, Tropics team)
 2C  Tapenade 3.4 (r3375) - 10 Feb 2010 15:08
 3C
 4C  Differentiation of inv in reverse (adjoint) mode:
 5C   gradient     of useful results: qt
 6C   with respect to varying inputs: qt a
 7C   RW status of diff variables: qt:in-out a:out
 8      SUBROUTINE INV_B(a, ab, qt, qtb, r, na)
 9      IMPLICIT NONE
10C  Hint: SIZE1OFrbINinv should be the value of na*na
11      INTEGER*4 na
12      REAL*8 a(na*na), qt(na*na), r(na*na)
13      REAL*8 ab(na*na), qtb(na*na)
14      INTEGER*4 n, m, k
15      REAL*8 tmp
16      REAL*8 tmpb
17      REAL*8 tmp0
18      INTEGER ad_from
19      REAL*8 tmp0b
20      REAL*8 rb(na*na)
21      INTEGER ii1
22      CALL PUSHREAL8ARRAY(r, na**2)
23      CALL PUSHREAL8ARRAY(qt, na**2)
24      CALL QR(a, qt, r, na)
25      DO n=na-1,0,-1
26        CALL PUSHREAL8(tmp)
27        tmp = r(n*na+n+1)
28        DO m=0,na-1
29          CALL PUSHREAL8(r(n*na+m+1))
30          r(n*na+m+1) = r(n*na+m+1)/tmp
31          CALL PUSHREAL8(qt(n*na+m+1))
32          qt(n*na+m+1) = qt(n*na+m+1)/tmp
33        ENDDO
34        ad_from = n + 1
35        DO m=ad_from,na-1
36          CALL PUSHREAL8(tmp)
37          tmp = r(n*na+m+1)
38          r(n*na+m+1) = 0
39          DO k=0,na-1
40            tmp0 = qt(n*na+k+1) - qt(m*na+k+1)*tmp
41            CALL PUSHREAL8(qt(n*na+k+1))
42            qt(n*na+k+1) = tmp0
43          ENDDO
44        ENDDO
45        CALL PUSHINTEGER4(ad_from)
46      ENDDO
47      DO ii1=1,na*na
48        rb(ii1) = 0.0
49      ENDDO
50      DO n=0,na-1,1
51        CALL POPINTEGER4(ad_from)
52        DO m=na-1,ad_from,-1
53          tmpb = 0.0
54          DO k=na-1,0,-1
55            CALL POPREAL8(qt(n*na+k+1))
56            tmp0b = qtb(n*na+k+1)
57            qtb(n*na+k+1) = tmp0b
58            qtb(m*na+k+1) = qtb(m*na+k+1) - tmp*tmp0b
59            tmpb = tmpb - qt(m*na+k+1)*tmp0b
60          ENDDO
61          rb(n*na+m+1) = tmpb
62          CALL POPREAL8(tmp)
63        ENDDO
64        tmpb = 0.0
65        DO m=na-1,0,-1
66          CALL POPREAL8(qt(n*na+m+1))
67          CALL POPREAL8(r(n*na+m+1))
68          tmpb = tmpb - r(n*na+m+1)*rb(n*na+m+1)/tmp**2 - qt(n*na+m+1)*
69     +      qtb(n*na+m+1)/tmp**2
70          qtb(n*na+m+1) = qtb(n*na+m+1)/tmp
71          rb(n*na+m+1) = rb(n*na+m+1)/tmp
72        ENDDO
73        CALL POPREAL8(tmp)
74        rb(n*na+n+1) = rb(n*na+n+1) + tmpb
75      ENDDO
76      CALL POPREAL8ARRAY(qt, na**2)
77      CALL POPREAL8ARRAY(r, na**2)
78      CALL QR_B(a, ab, qt, qtb, r, rb, na)
79      END