/algorithms_for_tapenade/fortran_code/qr_b.f
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 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