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 / lapack / dlasr.f < prev    next >
Text File  |  1996-09-28  |  12KB  |  326 lines
  1.       SUBROUTINE DLASR( SIDE, PIVOT, DIRECT, M, N, C, S, A, LDA )
  2. *
  3. *  -- LAPACK auxiliary routine (version 2.0) --
  4. *     Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
  5. *     Courant Institute, Argonne National Lab, and Rice University
  6. *     October 31, 1992
  7. *
  8. *     .. Scalar Arguments ..
  9.       CHARACTER          DIRECT, PIVOT, SIDE
  10.       INTEGER            LDA, M, N
  11. *     ..
  12. *     .. Array Arguments ..
  13.       DOUBLE PRECISION   A( LDA, * ), C( * ), S( * )
  14. *     ..
  15. *
  16. *  Purpose
  17. *  =======
  18. *
  19. *  DLASR   performs the transformation
  20. *
  21. *     A := P*A,   when SIDE = 'L' or 'l'  (  Left-hand side )
  22. *
  23. *     A := A*P',  when SIDE = 'R' or 'r'  ( Right-hand side )
  24. *
  25. *  where A is an m by n real matrix and P is an orthogonal matrix,
  26. *  consisting of a sequence of plane rotations determined by the
  27. *  parameters PIVOT and DIRECT as follows ( z = m when SIDE = 'L' or 'l'
  28. *  and z = n when SIDE = 'R' or 'r' ):
  29. *
  30. *  When  DIRECT = 'F' or 'f'  ( Forward sequence ) then
  31. *
  32. *     P = P( z - 1 )*...*P( 2 )*P( 1 ),
  33. *
  34. *  and when DIRECT = 'B' or 'b'  ( Backward sequence ) then
  35. *
  36. *     P = P( 1 )*P( 2 )*...*P( z - 1 ),
  37. *
  38. *  where  P( k ) is a plane rotation matrix for the following planes:
  39. *
  40. *     when  PIVOT = 'V' or 'v'  ( Variable pivot ),
  41. *        the plane ( k, k + 1 )
  42. *
  43. *     when  PIVOT = 'T' or 't'  ( Top pivot ),
  44. *        the plane ( 1, k + 1 )
  45. *
  46. *     when  PIVOT = 'B' or 'b'  ( Bottom pivot ),
  47. *        the plane ( k, z )
  48. *
  49. *  c( k ) and s( k )  must contain the  cosine and sine that define the
  50. *  matrix  P( k ).  The two by two plane rotation part of the matrix
  51. *  P( k ), R( k ), is assumed to be of the form
  52. *
  53. *     R( k ) = (  c( k )  s( k ) ).
  54. *              ( -s( k )  c( k ) )
  55. *
  56. *  This version vectorises across rows of the array A when SIDE = 'L'.
  57. *
  58. *  Arguments
  59. *  =========
  60. *
  61. *  SIDE    (input) CHARACTER*1
  62. *          Specifies whether the plane rotation matrix P is applied to
  63. *          A on the left or the right.
  64. *          = 'L':  Left, compute A := P*A
  65. *          = 'R':  Right, compute A:= A*P'
  66. *
  67. *  DIRECT  (input) CHARACTER*1
  68. *          Specifies whether P is a forward or backward sequence of
  69. *          plane rotations.
  70. *          = 'F':  Forward, P = P( z - 1 )*...*P( 2 )*P( 1 )
  71. *          = 'B':  Backward, P = P( 1 )*P( 2 )*...*P( z - 1 )
  72. *
  73. *  PIVOT   (input) CHARACTER*1
  74. *          Specifies the plane for which P(k) is a plane rotation
  75. *          matrix.
  76. *          = 'V':  Variable pivot, the plane (k,k+1)
  77. *          = 'T':  Top pivot, the plane (1,k+1)
  78. *          = 'B':  Bottom pivot, the plane (k,z)
  79. *
  80. *  M       (input) INTEGER
  81. *          The number of rows of the matrix A.  If m <= 1, an immediate
  82. *          return is effected.
  83. *
  84. *  N       (input) INTEGER
  85. *          The number of columns of the matrix A.  If n <= 1, an
  86. *          immediate return is effected.
  87. *
  88. *  C, S    (input) DOUBLE PRECISION arrays, dimension
  89. *                  (M-1) if SIDE = 'L'
  90. *                  (N-1) if SIDE = 'R'
  91. *          c(k) and s(k) contain the cosine and sine that define the
  92. *          matrix P(k).  The two by two plane rotation part of the
  93. *          matrix P(k), R(k), is assumed to be of the form
  94. *          R( k ) = (  c( k )  s( k ) ).
  95. *                   ( -s( k )  c( k ) )
  96. *
  97. *  A       (input/output) DOUBLE PRECISION array, dimension (LDA,N)
  98. *          The m by n matrix A.  On exit, A is overwritten by P*A if
  99. *          SIDE = 'R' or by A*P' if SIDE = 'L'.
  100. *
  101. *  LDA     (input) INTEGER
  102. *          The leading dimension of the array A.  LDA >= max(1,M).
  103. *
  104. *  =====================================================================
  105. *
  106. *     .. Parameters ..
  107.       DOUBLE PRECISION   ONE, ZERO
  108.       PARAMETER          ( ONE = 1.0D+0, ZERO = 0.0D+0 )
  109. *     ..
  110. *     .. Local Scalars ..
  111.       INTEGER            I, INFO, J
  112.       DOUBLE PRECISION   CTEMP, STEMP, TEMP
  113. *     ..
  114. *     .. External Functions ..
  115.       LOGICAL            LSAME
  116.       EXTERNAL           LSAME
  117. *     ..
  118. *     .. External Subroutines ..
  119.       EXTERNAL           XERBLA
  120. *     ..
  121. *     .. Intrinsic Functions ..
  122.       INTRINSIC          MAX
  123. *     ..
  124. *     .. Executable Statements ..
  125. *
  126. *     Test the input parameters
  127. *
  128.       INFO = 0
  129.       IF( .NOT.( LSAME( SIDE, 'L' ) .OR. LSAME( SIDE, 'R' ) ) ) THEN
  130.          INFO = 1
  131.       ELSE IF( .NOT.( LSAME( PIVOT, 'V' ) .OR. LSAME( PIVOT,
  132.      $         'T' ) .OR. LSAME( PIVOT, 'B' ) ) ) THEN
  133.          INFO = 2
  134.       ELSE IF( .NOT.( LSAME( DIRECT, 'F' ) .OR. LSAME( DIRECT, 'B' ) ) )
  135.      $          THEN
  136.          INFO = 3
  137.       ELSE IF( M.LT.0 ) THEN
  138.          INFO = 4
  139.       ELSE IF( N.LT.0 ) THEN
  140.          INFO = 5
  141.       ELSE IF( LDA.LT.MAX( 1, M ) ) THEN
  142.          INFO = 9
  143.       END IF
  144.       IF( INFO.NE.0 ) THEN
  145.          CALL XERBLA( 'DLASR ', INFO )
  146.          RETURN
  147.       END IF
  148. *
  149. *     Quick return if possible
  150. *
  151.       IF( ( M.EQ.0 ) .OR. ( N.EQ.0 ) )
  152.      $   RETURN
  153.       IF( LSAME( SIDE, 'L' ) ) THEN
  154. *
  155. *        Form  P * A
  156. *
  157.          IF( LSAME( PIVOT, 'V' ) ) THEN
  158.             IF( LSAME( DIRECT, 'F' ) ) THEN
  159.                DO 20 J = 1, M - 1
  160.                   CTEMP = C( J )
  161.                   STEMP = S( J )
  162.                   IF( ( CTEMP.NE.ONE ) .OR. ( STEMP.NE.ZERO ) ) THEN
  163.                      DO 10 I = 1, N
  164.                         TEMP = A( J+1, I )
  165.                         A( J+1, I ) = CTEMP*TEMP - STEMP*A( J, I )
  166.                         A( J, I ) = STEMP*TEMP + CTEMP*A( J, I )
  167.    10                CONTINUE
  168.                   END IF
  169.    20          CONTINUE
  170.             ELSE IF( LSAME( DIRECT, 'B' ) ) THEN
  171.                DO 40 J = M - 1, 1, -1
  172.                   CTEMP = C( J )
  173.                   STEMP = S( J )
  174.                   IF( ( CTEMP.NE.ONE ) .OR. ( STEMP.NE.ZERO ) ) THEN
  175.                      DO 30 I = 1, N
  176.                         TEMP = A( J+1, I )
  177.                         A( J+1, I ) = CTEMP*TEMP - STEMP*A( J, I )
  178.                         A( J, I ) = STEMP*TEMP + CTEMP*A( J, I )
  179.    30                CONTINUE
  180.                   END IF
  181.    40          CONTINUE
  182.             END IF
  183.          ELSE IF( LSAME( PIVOT, 'T' ) ) THEN
  184.             IF( LSAME( DIRECT, 'F' ) ) THEN
  185.                DO 60 J = 2, M
  186.                   CTEMP = C( J-1 )
  187.                   STEMP = S( J-1 )
  188.                   IF( ( CTEMP.NE.ONE ) .OR. ( STEMP.NE.ZERO ) ) THEN
  189.                      DO 50 I = 1, N
  190.                         TEMP = A( J, I )
  191.                         A( J, I ) = CTEMP*TEMP - STEMP*A( 1, I )
  192.                         A( 1, I ) = STEMP*TEMP + CTEMP*A( 1, I )
  193.    50                CONTINUE
  194.                   END IF
  195.    60          CONTINUE
  196.             ELSE IF( LSAME( DIRECT, 'B' ) ) THEN
  197.                DO 80 J = M, 2, -1
  198.                   CTEMP = C( J-1 )
  199.                   STEMP = S( J-1 )
  200.                   IF( ( CTEMP.NE.ONE ) .OR. ( STEMP.NE.ZERO ) ) THEN
  201.                      DO 70 I = 1, N
  202.                         TEMP = A( J, I )
  203.                         A( J, I ) = CTEMP*TEMP - STEMP*A( 1, I )
  204.                         A( 1, I ) = STEMP*TEMP + CTEMP*A( 1, I )
  205.    70                CONTINUE
  206.                   END IF
  207.    80          CONTINUE
  208.             END IF
  209.          ELSE IF( LSAME( PIVOT, 'B' ) ) THEN
  210.             IF( LSAME( DIRECT, 'F' ) ) THEN
  211.                DO 100 J = 1, M - 1
  212.                   CTEMP = C( J )
  213.                   STEMP = S( J )
  214.                   IF( ( CTEMP.NE.ONE ) .OR. ( STEMP.NE.ZERO ) ) THEN
  215.                      DO 90 I = 1, N
  216.                         TEMP = A( J, I )
  217.                         A( J, I ) = STEMP*A( M, I ) + CTEMP*TEMP
  218.                         A( M, I ) = CTEMP*A( M, I ) - STEMP*TEMP
  219.    90                CONTINUE
  220.                   END IF
  221.   100          CONTINUE
  222.             ELSE IF( LSAME( DIRECT, 'B' ) ) THEN
  223.                DO 120 J = M - 1, 1, -1
  224.                   CTEMP = C( J )
  225.                   STEMP = S( J )
  226.                   IF( ( CTEMP.NE.ONE ) .OR. ( STEMP.NE.ZERO ) ) THEN
  227.                      DO 110 I = 1, N
  228.                         TEMP = A( J, I )
  229.                         A( J, I ) = STEMP*A( M, I ) + CTEMP*TEMP
  230.                         A( M, I ) = CTEMP*A( M, I ) - STEMP*TEMP
  231.   110                CONTINUE
  232.                   END IF
  233.   120          CONTINUE
  234.             END IF
  235.          END IF
  236.       ELSE IF( LSAME( SIDE, 'R' ) ) THEN
  237. *
  238. *        Form A * P'
  239. *
  240.          IF( LSAME( PIVOT, 'V' ) ) THEN
  241.             IF( LSAME( DIRECT, 'F' ) ) THEN
  242.                DO 140 J = 1, N - 1
  243.                   CTEMP = C( J )
  244.                   STEMP = S( J )
  245.                   IF( ( CTEMP.NE.ONE ) .OR. ( STEMP.NE.ZERO ) ) THEN
  246.                      DO 130 I = 1, M
  247.                         TEMP = A( I, J+1 )
  248.                         A( I, J+1 ) = CTEMP*TEMP - STEMP*A( I, J )
  249.                         A( I, J ) = STEMP*TEMP + CTEMP*A( I, J )
  250.   130                CONTINUE
  251.                   END IF
  252.   140          CONTINUE
  253.             ELSE IF( LSAME( DIRECT, 'B' ) ) THEN
  254.                DO 160 J = N - 1, 1, -1
  255.                   CTEMP = C( J )
  256.                   STEMP = S( J )
  257.                   IF( ( CTEMP.NE.ONE ) .OR. ( STEMP.NE.ZERO ) ) THEN
  258.                      DO 150 I = 1, M
  259.                         TEMP = A( I, J+1 )
  260.                         A( I, J+1 ) = CTEMP*TEMP - STEMP*A( I, J )
  261.                         A( I, J ) = STEMP*TEMP + CTEMP*A( I, J )
  262.   150                CONTINUE
  263.                   END IF
  264.   160          CONTINUE
  265.             END IF
  266.          ELSE IF( LSAME( PIVOT, 'T' ) ) THEN
  267.             IF( LSAME( DIRECT, 'F' ) ) THEN
  268.                DO 180 J = 2, N
  269.                   CTEMP = C( J-1 )
  270.                   STEMP = S( J-1 )
  271.                   IF( ( CTEMP.NE.ONE ) .OR. ( STEMP.NE.ZERO ) ) THEN
  272.                      DO 170 I = 1, M
  273.                         TEMP = A( I, J )
  274.                         A( I, J ) = CTEMP*TEMP - STEMP*A( I, 1 )
  275.                         A( I, 1 ) = STEMP*TEMP + CTEMP*A( I, 1 )
  276.   170                CONTINUE
  277.                   END IF
  278.   180          CONTINUE
  279.             ELSE IF( LSAME( DIRECT, 'B' ) ) THEN
  280.                DO 200 J = N, 2, -1
  281.                   CTEMP = C( J-1 )
  282.                   STEMP = S( J-1 )
  283.                   IF( ( CTEMP.NE.ONE ) .OR. ( STEMP.NE.ZERO ) ) THEN
  284.                      DO 190 I = 1, M
  285.                         TEMP = A( I, J )
  286.                         A( I, J ) = CTEMP*TEMP - STEMP*A( I, 1 )
  287.                         A( I, 1 ) = STEMP*TEMP + CTEMP*A( I, 1 )
  288.   190                CONTINUE
  289.                   END IF
  290.   200          CONTINUE
  291.             END IF
  292.          ELSE IF( LSAME( PIVOT, 'B' ) ) THEN
  293.             IF( LSAME( DIRECT, 'F' ) ) THEN
  294.                DO 220 J = 1, N - 1
  295.                   CTEMP = C( J )
  296.                   STEMP = S( J )
  297.                   IF( ( CTEMP.NE.ONE ) .OR. ( STEMP.NE.ZERO ) ) THEN
  298.                      DO 210 I = 1, M
  299.                         TEMP = A( I, J )
  300.                         A( I, J ) = STEMP*A( I, N ) + CTEMP*TEMP
  301.                         A( I, N ) = CTEMP*A( I, N ) - STEMP*TEMP
  302.   210                CONTINUE
  303.                   END IF
  304.   220          CONTINUE
  305.             ELSE IF( LSAME( DIRECT, 'B' ) ) THEN
  306.                DO 240 J = N - 1, 1, -1
  307.                   CTEMP = C( J )
  308.                   STEMP = S( J )
  309.                   IF( ( CTEMP.NE.ONE ) .OR. ( STEMP.NE.ZERO ) ) THEN
  310.                      DO 230 I = 1, M
  311.                         TEMP = A( I, J )
  312.                         A( I, J ) = STEMP*A( I, N ) + CTEMP*TEMP
  313.                         A( I, N ) = CTEMP*A( I, N ) - STEMP*TEMP
  314.   230                CONTINUE
  315.                   END IF
  316.   240          CONTINUE
  317.             END IF
  318.          END IF
  319.       END IF
  320. *
  321.       RETURN
  322. *
  323. *     End of DLASR
  324. *
  325.       END
  326.