8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
24
25 implicit none
26
27
28
29 integer, intent(in) :: innLoop
30
31 real(r8), dimension(innLoop,innLoop), intent(inout) :: a, tl_a
32 real(r8), dimension(innLoop), intent(inout) :: tau, tl_tau
33 real(r8), dimension(innLoop), intent(inout) :: y, tl_y
34
35
36
37 integer :: i, j, ii, jj
38
39 real(r8) :: znorm, zbeta, zaii, ztemp
40 real(r8) :: tl_znorm, tl_zbeta, tl_zaii, tl_ztemp
41
42
43
44
45
46 DO i=1,innloop
47 tau(i)=0.0_r8
48 tl_tau(i)=0.0_r8
49 END DO
50
51 DO ii=1,innloop-1
52
53
54
55 znorm=0.0_r8
56 tl_znorm=0.0_r8
57 DO j=ii+1,innloop
58 znorm=znorm+a(ii,j)*a(ii,j)
59 tl_znorm=tl_znorm+2.0_r8*a(ii,j)*tl_a(ii,j)
60 END DO
61 znorm=sqrt(znorm)
62 IF (znorm.NE.0.0_r8) THEN
63 tl_znorm=0.5_r8*tl_znorm/znorm
64 ELSE
65 tl_znorm=0.0_r8
66 END IF
67 zbeta=sqrt(znorm*znorm+a(ii,ii)*a(ii,ii))
68 IF (zbeta.NE.0.0_r8) THEN
69 tl_zbeta=(tl_znorm*znorm+tl_a(ii,ii)*a(ii,ii))/zbeta
70 ELSE
71 tl_zbeta=0.0_r8
72 END IF
73
74
75 tl_zbeta=-sign(1.0_r8,zbeta)*sign(1.0_r8,a(ii,ii))*tl_zbeta
76 zbeta=-sign(zbeta,a(ii,ii))
77 tau(ii)=(zbeta-a(ii,ii))/zbeta
78 tl_tau(ii)=(tl_zbeta-tl_a(ii,ii))/zbeta-tl_zbeta*tau(ii)/zbeta
79 DO j=ii+1,innloop
80 a(ii,j)=a(ii,j)/(a(ii,ii)-zbeta)
81 tl_a(ii,j)=tl_a(ii,j)/(a(ii,ii)-zbeta)- &
82 & (tl_a(ii,ii)-tl_zbeta)*a(ii,j)/(a(ii,ii)-zbeta)
83 END DO
84 a(ii,ii)=zbeta
85 tl_a(ii,ii)=tl_zbeta
86
87
88
89 zaii=a(ii,ii)
90 tl_zaii=tl_a(ii,ii)
91 a(ii,ii)=1.0_r8
92 tl_a(ii,ii)=0.0_r8
93 IF (tau(ii).ne.0.0_r8) THEN
94 DO i=1,innloop
95 y(i)=0.0_r8
96 tl_y(i)=0.0_r8
97 END DO
98 jj=ii-1
99 DO j=1,innloop-jj
100 ztemp=a(ii,j+jj)
101 tl_ztemp=tl_a(ii,j+jj)
102 DO i=1,innloop-ii
103 y(i)=y(i)+ztemp*a(ii+i,j+jj)
104 tl_y(i)=tl_y(i)+tl_ztemp*a(ii+i,j+jj)+ &
105 & ztemp*tl_a(ii+i,j+jj)
106 END DO
107 END DO
108
109 DO j=1,innloop-jj
110 ztemp=-tau(ii)*a(ii,j+jj)
111 tl_ztemp=-tl_tau(ii)*a(ii,j+jj)-tau(ii)*tl_a(ii,j+jj)
112 DO i=1,innloop-ii
113 a(ii+i,j+jj)=a(ii+i,j+jj)+y(i)*ztemp
114 tl_a(ii+i,j+jj)=tl_a(ii+i,j+jj)+tl_y(i)*ztemp+ &
115 & y(i)*tl_ztemp
116 END DO
117 END DO
118 END IF
119 a(ii,ii)=zaii
120 tl_a(ii,ii)=tl_zaii
121 END DO
122
123 RETURN