/algorithms_for_tapenade/fortran_code/inv_b.f

http://github.com/b45ch1/hpsc_hanoi_2009_walter · FORTRAN Legacy · 79 lines · 53 code · 0 blank · 26 comment · 0 complexity · 63f57b025270c75988d7c079c9bac788 MD5 · raw file 

    

  1. C        Generated by TAPENADE     (INRIA, Tropics team)
    2. 

  2. C  Tapenade 3.4 (r3375) - 10 Feb 2010 15:08
    3. 

  3. C
    4. 

  4. C  Differentiation of inv in reverse (adjoint) mode:
    5. 

  5. C   gradient     of useful results: qt
    6. 

  6. C   with respect to varying inputs: qt a
    7. 

  7. C   RW status of diff variables: qt:in-out a:out
    8. 

  8.       SUBROUTINE INV_B(a, ab, qt, qtb, r, na)
    9. 

  9.       IMPLICIT NONE
    10. 

  10. C  Hint: SIZE1OFrbINinv should be the value of na*na
    11. 

  11.       INTEGER*4 na
    12. 

  12.       REAL*8 a(na*na), qt(na*na), r(na*na)
    13. 

  13.       REAL*8 ab(na*na), qtb(na*na)
    14. 

  14.       INTEGER*4 n, m, k
    15. 

  15.       REAL*8 tmp
    16. 

  16.       REAL*8 tmpb
    17. 

  17.       REAL*8 tmp0
    18. 

  18.       INTEGER ad_from
    19. 

  19.       REAL*8 tmp0b
    20. 

  20.       REAL*8 rb(na*na)
    21. 

  21.       INTEGER ii1
    22. 

  22.       CALL PUSHREAL8ARRAY(r, na**2)
    23. 

  23.       CALL PUSHREAL8ARRAY(qt, na**2)
    24. 

  24.       CALL QR(a, qt, r, na)
    25. 

  25.       DO n=na-1,0,-1
    26. 

  26.         CALL PUSHREAL8(tmp)
    27. 

  27.         tmp = r(n*na+n+1)
    28. 

  28.         DO m=0,na-1
    29. 

  29.           CALL PUSHREAL8(r(n*na+m+1))
    30. 

  30.           r(n*na+m+1) = r(n*na+m+1)/tmp
    31. 

  31.           CALL PUSHREAL8(qt(n*na+m+1))
    32. 

  32.           qt(n*na+m+1) = qt(n*na+m+1)/tmp
    33. 

  33.         ENDDO
    34. 

  34.         ad_from = n + 1
    35. 

  35.         DO m=ad_from,na-1
    36. 

  36.           CALL PUSHREAL8(tmp)
    37. 

  37.           tmp = r(n*na+m+1)
    38. 

  38.           r(n*na+m+1) = 0
    39. 

  39.           DO k=0,na-1
    40. 

  40.             tmp0 = qt(n*na+k+1) - qt(m*na+k+1)*tmp
    41. 

  41.             CALL PUSHREAL8(qt(n*na+k+1))
    42. 

  42.             qt(n*na+k+1) = tmp0
    43. 

  43.           ENDDO
    44. 

  44.         ENDDO
    45. 

  45.         CALL PUSHINTEGER4(ad_from)
    46. 

  46.       ENDDO
    47. 

  47.       DO ii1=1,na*na
    48. 

  48.         rb(ii1) = 0.0
    49. 

  49.       ENDDO
    50. 

  50.       DO n=0,na-1,1
    51. 

  51.         CALL POPINTEGER4(ad_from)
    52. 

  52.         DO m=na-1,ad_from,-1
    53. 

  53.           tmpb = 0.0
    54. 

  54.           DO k=na-1,0,-1
    55. 

  55.             CALL POPREAL8(qt(n*na+k+1))
    56. 

  56.             tmp0b = qtb(n*na+k+1)
    57. 

  57.             qtb(n*na+k+1) = tmp0b
    58. 

  58.             qtb(m*na+k+1) = qtb(m*na+k+1) - tmp*tmp0b
    59. 

  59.             tmpb = tmpb - qt(m*na+k+1)*tmp0b
    60. 

  60.           ENDDO
    61. 

  61.           rb(n*na+m+1) = tmpb
    62. 

  62.           CALL POPREAL8(tmp)
    63. 

  63.         ENDDO
    64. 

  64.         tmpb = 0.0
    65. 

  65.         DO m=na-1,0,-1
    66. 

  66.           CALL POPREAL8(qt(n*na+m+1))
    67. 

  67.           CALL POPREAL8(r(n*na+m+1))
    68. 

  68.           tmpb = tmpb - r(n*na+m+1)*rb(n*na+m+1)/tmp**2 - qt(n*na+m+1)*
    69. 

  69.      +      qtb(n*na+m+1)/tmp**2
    70. 

  70.           qtb(n*na+m+1) = qtb(n*na+m+1)/tmp
    71. 

  71.           rb(n*na+m+1) = rb(n*na+m+1)/tmp
    72. 

  72.         ENDDO
    73. 

  73.         CALL POPREAL8(tmp)
    74. 

  74.         rb(n*na+n+1) = rb(n*na+n+1) + tmpb
    75. 

  75.       ENDDO
    76. 

  76.       CALL POPREAL8ARRAY(qt, na**2)
    77. 

  77.       CALL POPREAL8ARRAY(r, na**2)
    78. 

  78.       CALL QR_B(a, ab, qt, qtb, r, rb, na)
    79. 

  79.       END
    80.