/algorithms_for_tapenade/fortran_code/qr_b.f

  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
 17      INTEGER ad_from
 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
 40        DO m=ad_from,na-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
 70        CALL PUSHINTEGER4(ad_from)
 71      ENDDO
 72      DO n=na-1,0,-1
 73        CALL POPINTEGER4(ad_from)
 74        DO m=na-1,ad_from,-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