home *** CD-ROM | disk | FTP | other *** search
/ Geek Gadgets 1 / ADE-1.bin / ade-dist / octave-1.1.1p1-src.tgz / tar.out / fsf / octave / libcruft / minpack / r1mpyq.f < prev    next >
Text File  |  1996-09-28  |  3KB  |  93 lines
  1.       SUBROUTINE R1MPYQ(M,N,A,LDA,V,W)
  2.       INTEGER M,N,LDA
  3.       DOUBLE PRECISION A(LDA,N),V(N),W(N)
  4. C     **********
  5. C
  6. C     SUBROUTINE R1MPYQ
  7. C
  8. C     GIVEN AN M BY N MATRIX A, THIS SUBROUTINE COMPUTES A*Q WHERE
  9. C     Q IS THE PRODUCT OF 2*(N - 1) TRANSFORMATIONS
  10. C
  11. C           GV(N-1)*...*GV(1)*GW(1)*...*GW(N-1)
  12. C
  13. C     AND GV(I), GW(I) ARE GIVENS ROTATIONS IN THE (I,N) PLANE WHICH
  14. C     ELIMINATE ELEMENTS IN THE I-TH AND N-TH PLANES, RESPECTIVELY.
  15. C     Q ITSELF IS NOT GIVEN, RATHER THE INFORMATION TO RECOVER THE
  16. C     GV, GW ROTATIONS IS SUPPLIED.
  17. C
  18. C     THE SUBROUTINE STATEMENT IS
  19. C
  20. C       SUBROUTINE R1MPYQ(M,N,A,LDA,V,W)
  21. C
  22. C     WHERE
  23. C
  24. C       M IS A POSITIVE INTEGER INPUT VARIABLE SET TO THE NUMBER
  25. C         OF ROWS OF A.
  26. C
  27. C       N IS A POSITIVE INTEGER INPUT VARIABLE SET TO THE NUMBER
  28. C         OF COLUMNS OF A.
  29. C
  30. C       A IS AN M BY N ARRAY. ON INPUT A MUST CONTAIN THE MATRIX
  31. C         TO BE POSTMULTIPLIED BY THE ORTHOGONAL MATRIX Q
  32. C         DESCRIBED ABOVE. ON OUTPUT A*Q HAS REPLACED A.
  33. C
  34. C       LDA IS A POSITIVE INTEGER INPUT VARIABLE NOT LESS THAN M
  35. C         WHICH SPECIFIES THE LEADING DIMENSION OF THE ARRAY A.
  36. C
  37. C       V IS AN INPUT ARRAY OF LENGTH N. V(I) MUST CONTAIN THE
  38. C         INFORMATION NECESSARY TO RECOVER THE GIVENS ROTATION GV(I)
  39. C         DESCRIBED ABOVE.
  40. C
  41. C       W IS AN INPUT ARRAY OF LENGTH N. W(I) MUST CONTAIN THE
  42. C         INFORMATION NECESSARY TO RECOVER THE GIVENS ROTATION GW(I)
  43. C         DESCRIBED ABOVE.
  44. C
  45. C     SUBROUTINES CALLED
  46. C
  47. C       FORTRAN-SUPPLIED ... DABS,DSQRT
  48. C
  49. C     MINPACK. VERSION OF DECEMBER 1978.
  50. C     BURTON S. GARBOW, KENNETH E. HILLSTROM, JORGE J. MORE
  51. C
  52. C     **********
  53.       INTEGER I,J,NMJ,NM1
  54.       DOUBLE PRECISION COS,ONE,SIN,TEMP
  55.       DATA ONE /1.0D0/
  56. C
  57. C     APPLY THE FIRST SET OF GIVENS ROTATIONS TO A.
  58. C
  59.       NM1 = N - 1
  60.       IF (NM1 .LT. 1) GO TO 50
  61.       DO 20 NMJ = 1, NM1
  62.          J = N - NMJ
  63.          IF (DABS(V(J)) .GT. ONE) COS = ONE/V(J)
  64.          IF (DABS(V(J)) .GT. ONE) SIN = DSQRT(ONE-COS**2)
  65.          IF (DABS(V(J)) .LE. ONE) SIN = V(J)
  66.          IF (DABS(V(J)) .LE. ONE) COS = DSQRT(ONE-SIN**2)
  67.          DO 10 I = 1, M
  68.             TEMP = COS*A(I,J) - SIN*A(I,N)
  69.             A(I,N) = SIN*A(I,J) + COS*A(I,N)
  70.             A(I,J) = TEMP
  71.    10       CONTINUE
  72.    20    CONTINUE
  73. C
  74. C     APPLY THE SECOND SET OF GIVENS ROTATIONS TO A.
  75. C
  76.       DO 40 J = 1, NM1
  77.          IF (DABS(W(J)) .GT. ONE) COS = ONE/W(J)
  78.          IF (DABS(W(J)) .GT. ONE) SIN = DSQRT(ONE-COS**2)
  79.          IF (DABS(W(J)) .LE. ONE) SIN = W(J)
  80.          IF (DABS(W(J)) .LE. ONE) COS = DSQRT(ONE-SIN**2)
  81.          DO 30 I = 1, M
  82.             TEMP = COS*A(I,J) + SIN*A(I,N)
  83.             A(I,N) = -SIN*A(I,J) + COS*A(I,N)
  84.             A(I,J) = TEMP
  85.    30       CONTINUE
  86.    40    CONTINUE
  87.    50 CONTINUE
  88.       RETURN
  89. C
  90. C     LAST CARD OF SUBROUTINE R1MPYQ.
  91. C
  92.       END
  93.