http://github.com/b45ch1/hpsc_hanoi_2009_walter
FORTRAN Legacy | 135 lines | 95 code | 0 blank | 40 comment | 0 complexity | 24f491681f37aa519c8f2dacbe6ea13d MD5 | raw file
```  1C        Generated by TAPENADE     (INRIA, Tropics team)
2C  Tapenade 3.4 (r3375) - 10 Feb 2010 15:08
3C
4C  Differentiation of qr in reverse (adjoint) mode:
5C   gradient     of useful results: r qt
6C   with respect to varying inputs: qt a
7      SUBROUTINE QR_B(a, ab, qt, qtb, r, rb, na)
8      IMPLICIT NONE
9      INTEGER*4 na
10      REAL*8 a(na*na), qt(na*na), r(na*na)
11      REAL*8 ab(na*na), qtb(na*na), rb(na*na)
12      REAL*8 tmp, at, bt, rt, c, s, rnk, qtnk
13      REAL*8 atb, btb, rtb, cb, sb, rnkb, qtnkb
14      INTEGER*4 n, m, k
15      REAL*8 tmp0
16      REAL*8 tmp1
18      REAL*8 tmp0b
19      REAL*8 tempb
20      INTEGER ii1
21      INTRINSIC SQRT
22      REAL*8 tmp1b
23C     PREPARE Q AND R
24      DO n=0,na-1
25        DO m=0,na-1
26          IF (n .EQ. m) THEN
27            tmp = 1.0
28          ELSE
29            tmp = 0.0
30          END IF
31          qt(n*na+m+1) = tmp
32        ENDDO
33      ENDDO
34      DO n=1,na*na
35        r(n) = a(n)
36      ENDDO
37C	MAIN ALGORITHM
38      DO n=0,na-1
39        ad_from = n + 1
41          CALL PUSHREAL8(at)
42          at = r(n*na+n+1)
43          CALL PUSHREAL8(bt)
44          bt = r(m*na+n+1)
45          CALL PUSHREAL8(rt)
46          rt = SQRT(at*at + bt*bt)
47          CALL PUSHREAL8(c)
48          c = at/rt
49          CALL PUSHREAL8(s)
50          s = bt/rt
51          DO k=1,na
52            CALL PUSHREAL8(rnk)
53C				UPDATE R
54            rnk = r(n*na+k)
55            tmp0 = c*rnk + s*r(m*na+k)
56            CALL PUSHREAL8(r(n*na+k))
57            r(n*na+k) = tmp0
58            CALL PUSHREAL8(r(m*na+k))
59            r(m*na+k) = -(s*rnk) + c*r(m*na+k)
60            CALL PUSHREAL8(qtnk)
61C				UPDATE Q
62            qtnk = qt(n*na+k)
63            tmp1 = c*qtnk + s*qt(m*na+k)
64            CALL PUSHREAL8(qt(n*na+k))
65            qt(n*na+k) = tmp1
66            CALL PUSHREAL8(qt(m*na+k))
67            qt(m*na+k) = -(s*qtnk) + c*qt(m*na+k)
68          ENDDO
69        ENDDO
71      ENDDO
72      DO n=na-1,0,-1
75          sb = 0.0
76          cb = 0.0
77          DO k=na,1,-1
78            CALL POPREAL8(qt(m*na+k))
79            cb = cb + qt(m*na+k)*qtb(m*na+k)
80            sb = sb - qtnk*qtb(m*na+k)
81            qtnkb = -(s*qtb(m*na+k))
82            qtb(m*na+k) = c*qtb(m*na+k)
83            CALL POPREAL8(qt(n*na+k))
84            tmp1b = qtb(n*na+k)
85            qtb(n*na+k) = 0.0
86            cb = cb + qtnk*tmp1b
87            qtnkb = qtnkb + c*tmp1b
88            sb = sb + qt(m*na+k)*tmp1b - rnk*rb(m*na+k)
89            qtb(m*na+k) = qtb(m*na+k) + s*tmp1b
90            CALL POPREAL8(qtnk)
91            qtb(n*na+k) = qtb(n*na+k) + qtnkb
92            CALL POPREAL8(r(m*na+k))
93            cb = cb + r(m*na+k)*rb(m*na+k)
94            rnkb = -(s*rb(m*na+k))
95            rb(m*na+k) = c*rb(m*na+k)
96            CALL POPREAL8(r(n*na+k))
97            tmp0b = rb(n*na+k)
98            rb(n*na+k) = 0.0
99            cb = cb + rnk*tmp0b
100            rnkb = rnkb + c*tmp0b
101            sb = sb + r(m*na+k)*tmp0b
102            rb(m*na+k) = rb(m*na+k) + s*tmp0b
103            CALL POPREAL8(rnk)
104            rb(n*na+k) = rb(n*na+k) + rnkb
105          ENDDO
106          CALL POPREAL8(s)
107          rtb = -(at*cb/rt**2) - bt*sb/rt**2
108          IF (at**2 + bt**2 .EQ. 0.0) THEN
109            tempb = 0.0
110          ELSE
111            tempb = rtb/(2.0*SQRT(at**2+bt**2))
112          END IF
113          btb = 2*bt*tempb + sb/rt
114          CALL POPREAL8(c)
115          atb = 2*at*tempb + cb/rt
116          CALL POPREAL8(rt)
117          CALL POPREAL8(bt)
118          rb(m*na+n+1) = rb(m*na+n+1) + btb
119          CALL POPREAL8(at)
120          rb(n*na+n+1) = rb(n*na+n+1) + atb
121        ENDDO
122      ENDDO
123      DO ii1=1,na*na
124        ab(ii1) = 0.0
125      ENDDO
126      DO n=na*na,1,-1
127        ab(n) = ab(n) + rb(n)
128        rb(n) = 0.0
129      ENDDO
130      DO n=na-1,0,-1
131        DO m=na-1,0,-1
132          qtb(n*na+m+1) = 0.0
133        ENDDO
134      ENDDO
135      END
```