You can not select more than 25 topics Topics must start with a letter or number, can include dashes ('-') and can be up to 35 characters long.

206 lines
4.8 KiB

  1. SUBROUTINE DROTMG(DD1,DD2,DX1,DY1,DPARAM)
  2. * .. Scalar Arguments ..
  3. DOUBLE PRECISION DD1,DD2,DX1,DY1
  4. * ..
  5. * .. Array Arguments ..
  6. DOUBLE PRECISION DPARAM(5)
  7. * ..
  8. *
  9. * Purpose
  10. * =======
  11. *
  12. * CONSTRUCT THE MODIFIED GIVENS TRANSFORMATION MATRIX H WHICH ZEROS
  13. * THE SECOND COMPONENT OF THE 2-VECTOR (DSQRT(DD1)*DX1,DSQRT(DD2)*
  14. * DY2)**T.
  15. * WITH DPARAM(1)=DFLAG, H HAS ONE OF THE FOLLOWING FORMS..
  16. *
  17. * DFLAG=-1.D0 DFLAG=0.D0 DFLAG=1.D0 DFLAG=-2.D0
  18. *
  19. * (DH11 DH12) (1.D0 DH12) (DH11 1.D0) (1.D0 0.D0)
  20. * H=( ) ( ) ( ) ( )
  21. * (DH21 DH22), (DH21 1.D0), (-1.D0 DH22), (0.D0 1.D0).
  22. * LOCATIONS 2-4 OF DPARAM CONTAIN DH11, DH21, DH12, AND DH22
  23. * RESPECTIVELY. (VALUES OF 1.D0, -1.D0, OR 0.D0 IMPLIED BY THE
  24. * VALUE OF DPARAM(1) ARE NOT STORED IN DPARAM.)
  25. *
  26. * THE VALUES OF GAMSQ AND RGAMSQ SET IN THE DATA STATEMENT MAY BE
  27. * INEXACT. THIS IS OK AS THEY ARE ONLY USED FOR TESTING THE SIZE
  28. * OF DD1 AND DD2. ALL ACTUAL SCALING OF DATA IS DONE USING GAM.
  29. *
  30. *
  31. * Arguments
  32. * =========
  33. *
  34. * DD1 (input/output) DOUBLE PRECISION
  35. *
  36. * DD2 (input/output) DOUBLE PRECISION
  37. *
  38. * DX1 (input/output) DOUBLE PRECISION
  39. *
  40. * DY1 (input) DOUBLE PRECISION
  41. *
  42. * DPARAM (input/output) DOUBLE PRECISION array, dimension 5
  43. * DPARAM(1)=DFLAG
  44. * DPARAM(2)=DH11
  45. * DPARAM(3)=DH21
  46. * DPARAM(4)=DH12
  47. * DPARAM(5)=DH22
  48. *
  49. * =====================================================================
  50. *
  51. * .. Local Scalars ..
  52. DOUBLE PRECISION DFLAG,DH11,DH12,DH21,DH22,DP1,DP2,DQ1,DQ2,DTEMP,
  53. + DU,GAM,GAMSQ,ONE,RGAMSQ,TWO,ZERO
  54. INTEGER IGO
  55. * ..
  56. * .. Intrinsic Functions ..
  57. INTRINSIC DABS
  58. * ..
  59. * .. Data statements ..
  60. *
  61. DATA ZERO,ONE,TWO/0.D0,1.D0,2.D0/
  62. DATA GAM,GAMSQ,RGAMSQ/4096.D0,16777216.D0,5.9604645D-8/
  63. * ..
  64. IF (.NOT.DD1.LT.ZERO) GO TO 10
  65. * GO ZERO-H-D-AND-DX1..
  66. GO TO 60
  67. 10 CONTINUE
  68. * CASE-DD1-NONNEGATIVE
  69. DP2 = DD2*DY1
  70. IF (.NOT.DP2.EQ.ZERO) GO TO 20
  71. DFLAG = -TWO
  72. GO TO 260
  73. * REGULAR-CASE..
  74. 20 CONTINUE
  75. DP1 = DD1*DX1
  76. DQ2 = DP2*DY1
  77. DQ1 = DP1*DX1
  78. *
  79. IF (.NOT.DABS(DQ1).GT.DABS(DQ2)) GO TO 40
  80. DH21 = -DY1/DX1
  81. DH12 = DP2/DP1
  82. *
  83. DU = ONE - DH12*DH21
  84. *
  85. IF (.NOT.DU.LE.ZERO) GO TO 30
  86. * GO ZERO-H-D-AND-DX1..
  87. GO TO 60
  88. 30 CONTINUE
  89. DFLAG = ZERO
  90. DD1 = DD1/DU
  91. DD2 = DD2/DU
  92. DX1 = DX1*DU
  93. * GO SCALE-CHECK..
  94. GO TO 100
  95. 40 CONTINUE
  96. IF (.NOT.DQ2.LT.ZERO) GO TO 50
  97. * GO ZERO-H-D-AND-DX1..
  98. GO TO 60
  99. 50 CONTINUE
  100. DFLAG = ONE
  101. DH11 = DP1/DP2
  102. DH22 = DX1/DY1
  103. DU = ONE + DH11*DH22
  104. DTEMP = DD2/DU
  105. DD2 = DD1/DU
  106. DD1 = DTEMP
  107. DX1 = DY1*DU
  108. * GO SCALE-CHECK
  109. GO TO 100
  110. * PROCEDURE..ZERO-H-D-AND-DX1..
  111. 60 CONTINUE
  112. DFLAG = -ONE
  113. DH11 = ZERO
  114. DH12 = ZERO
  115. DH21 = ZERO
  116. DH22 = ZERO
  117. *
  118. DD1 = ZERO
  119. DD2 = ZERO
  120. DX1 = ZERO
  121. * RETURN..
  122. GO TO 220
  123. * PROCEDURE..FIX-H..
  124. 70 CONTINUE
  125. IF (.NOT.DFLAG.GE.ZERO) GO TO 90
  126. *
  127. IF (.NOT.DFLAG.EQ.ZERO) GO TO 80
  128. DH11 = ONE
  129. DH22 = ONE
  130. DFLAG = -ONE
  131. GO TO 90
  132. 80 CONTINUE
  133. DH21 = -ONE
  134. DH12 = ONE
  135. DFLAG = -ONE
  136. 90 CONTINUE
  137. GO TO IGO(120,150,180,210)
  138. * PROCEDURE..SCALE-CHECK
  139. 100 CONTINUE
  140. 110 CONTINUE
  141. IF (.NOT.DD1.LE.RGAMSQ) GO TO 130
  142. IF (DD1.EQ.ZERO) GO TO 160
  143. ASSIGN 120 TO IGO
  144. * FIX-H..
  145. GO TO 70
  146. 120 CONTINUE
  147. DD1 = DD1*GAM**2
  148. DX1 = DX1/GAM
  149. DH11 = DH11/GAM
  150. DH12 = DH12/GAM
  151. GO TO 110
  152. 130 CONTINUE
  153. 140 CONTINUE
  154. IF (.NOT.DD1.GE.GAMSQ) GO TO 160
  155. ASSIGN 150 TO IGO
  156. * FIX-H..
  157. GO TO 70
  158. 150 CONTINUE
  159. DD1 = DD1/GAM**2
  160. DX1 = DX1*GAM
  161. DH11 = DH11*GAM
  162. DH12 = DH12*GAM
  163. GO TO 140
  164. 160 CONTINUE
  165. 170 CONTINUE
  166. IF (.NOT.DABS(DD2).LE.RGAMSQ) GO TO 190
  167. IF (DD2.EQ.ZERO) GO TO 220
  168. ASSIGN 180 TO IGO
  169. * FIX-H..
  170. GO TO 70
  171. 180 CONTINUE
  172. DD2 = DD2*GAM**2
  173. DH21 = DH21/GAM
  174. DH22 = DH22/GAM
  175. GO TO 170
  176. 190 CONTINUE
  177. 200 CONTINUE
  178. IF (.NOT.DABS(DD2).GE.GAMSQ) GO TO 220
  179. ASSIGN 210 TO IGO
  180. * FIX-H..
  181. GO TO 70
  182. 210 CONTINUE
  183. DD2 = DD2/GAM**2
  184. DH21 = DH21*GAM
  185. DH22 = DH22*GAM
  186. GO TO 200
  187. 220 CONTINUE
  188. IF (DFLAG) 250,230,240
  189. 230 CONTINUE
  190. DPARAM(3) = DH21
  191. DPARAM(4) = DH12
  192. GO TO 260
  193. 240 CONTINUE
  194. DPARAM(2) = DH11
  195. DPARAM(5) = DH22
  196. GO TO 260
  197. 250 CONTINUE
  198. DPARAM(2) = DH11
  199. DPARAM(3) = DH21
  200. DPARAM(4) = DH12
  201. DPARAM(5) = DH22
  202. 260 CONTINUE
  203. DPARAM(1) = DFLAG
  204. RETURN
  205. END