drot.f 3.0 KB

1234567891011121314151617181920212223242526272829303132333435363738394041424344454647484950515253545556575859606162636465666768697071727374757677787980818283848586878889
  1. *DECK DROT
  2. SUBROUTINE DROT (N, DX, INCX, DY, INCY, DC, DS)
  3. C***BEGIN PROLOGUE DROT
  4. C***PURPOSE Apply a plane Givens rotation.
  5. C***LIBRARY SLATEC (BLAS)
  6. C***CATEGORY D1A8
  7. C***TYPE DOUBLE PRECISION (SROT-S, DROT-D, CSROT-C)
  8. C***KEYWORDS BLAS, GIVENS ROTATION, GIVENS TRANSFORMATION,
  9. C LINEAR ALGEBRA, PLANE ROTATION, VECTOR
  10. C***AUTHOR Lawson, C. L., (JPL)
  11. C Hanson, R. J., (SNLA)
  12. C Kincaid, D. R., (U. of Texas)
  13. C Krogh, F. T., (JPL)
  14. C***DESCRIPTION
  15. C
  16. C B L A S Subprogram
  17. C Description of Parameters
  18. C
  19. C --Input--
  20. C N number of elements in input vector(s)
  21. C DX double precision vector with N elements
  22. C INCX storage spacing between elements of DX
  23. C DY double precision vector with N elements
  24. C INCY storage spacing between elements of DY
  25. C DC D.P. element of rotation matrix
  26. C DS D.P. element of rotation matrix
  27. C
  28. C --Output--
  29. C DX rotated vector DX (unchanged if N .LE. 0)
  30. C DY rotated vector DY (unchanged if N .LE. 0)
  31. C
  32. C Multiply the 2 x 2 matrix ( DC DS) times the 2 x N matrix (DX**T)
  33. C (-DS DC) (DY**T)
  34. C where **T indicates transpose. The elements of DX are in
  35. C DX(LX+I*INCX), I = 0 to N-1, where LX = 1 if INCX .GE. 0, else
  36. C LX = 1+(1-N)*INCX, and similarly for DY using LY and INCY.
  37. C
  38. C***REFERENCES C. L. Lawson, R. J. Hanson, D. R. Kincaid and F. T.
  39. C Krogh, Basic linear algebra subprograms for Fortran
  40. C usage, Algorithm No. 539, Transactions on Mathematical
  41. C Software 5, 3 (September 1979), pp. 308-323.
  42. C***ROUTINES CALLED (NONE)
  43. C***REVISION HISTORY (YYMMDD)
  44. C 791001 DATE WRITTEN
  45. C 861211 REVISION DATE from Version 3.2
  46. C 891214 Prologue converted to Version 4.0 format. (BAB)
  47. C 920310 Corrected definition of LX in DESCRIPTION. (WRB)
  48. C 920501 Reformatted the REFERENCES section. (WRB)
  49. C***END PROLOGUE DROT
  50. DOUBLE PRECISION DX, DY, DC, DS, ZERO, ONE, W, Z
  51. DIMENSION DX(*), DY(*)
  52. SAVE ZERO, ONE
  53. DATA ZERO, ONE /0.0D0, 1.0D0/
  54. C***FIRST EXECUTABLE STATEMENT DROT
  55. IF (N .LE. 0 .OR. (DS .EQ. ZERO .AND. DC .EQ. ONE)) GO TO 40
  56. IF (.NOT. (INCX .EQ. INCY .AND. INCX .GT. 0)) GO TO 20
  57. C
  58. C Code for equal and positive increments.
  59. C
  60. NSTEPS=INCX*N
  61. DO 10 I = 1,NSTEPS,INCX
  62. W=DX(I)
  63. Z=DY(I)
  64. DX(I)=DC*W+DS*Z
  65. DY(I)=-DS*W+DC*Z
  66. 10 CONTINUE
  67. GO TO 40
  68. C
  69. C Code for unequal or nonpositive increments.
  70. C
  71. 20 CONTINUE
  72. KX=1
  73. KY=1
  74. C
  75. IF (INCX .LT. 0) KX = 1-(N-1)*INCX
  76. IF (INCY .LT. 0) KY = 1-(N-1)*INCY
  77. C
  78. DO 30 I = 1,N
  79. W=DX(KX)
  80. Z=DY(KY)
  81. DX(KX)=DC*W+DS*Z
  82. DY(KY)=-DS*W+DC*Z
  83. KX=KX+INCX
  84. KY=KY+INCY
  85. 30 CONTINUE
  86. 40 CONTINUE
  87. C
  88. RETURN
  89. END