LAPACK: SRC/dlasr.f Source File

Go to the documentation of this file.
00001       SUBROUTINE DLASR( SIDE, PIVOT, DIRECT, M, N, C, S, A, LDA )
00002 *
00003 *  -- LAPACK auxiliary routine (version 3.2) --
00004 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
00005 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
00006 *     November 2006
00007 *
00008 *     .. Scalar Arguments ..
00009       CHARACTER          DIRECT, PIVOT, SIDE
00010       INTEGER            LDA, M, N
00011 *     ..
00012 *     .. Array Arguments ..
00013       DOUBLE PRECISION   A( LDA, * ), C( * ), S( * )
00014 *     ..
00015 *
00016 *  Purpose
00017 *  =======
00018 *
00019 *  DLASR applies a sequence of plane rotations to a real matrix A,
00020 *  from either the left or the right.
00021 *  
00022 *  When SIDE = 'L', the transformation takes the form
00023 *  
00024 *     A := P*A
00025 *  
00026 *  and when SIDE = 'R', the transformation takes the form
00027 *  
00028 *     A := A*P**T
00029 *  
00030 *  where P is an orthogonal matrix consisting of a sequence of z plane
00031 *  rotations, with z = M when SIDE = 'L' and z = N when SIDE = 'R',
00032 *  and P**T is the transpose of P.
00033 *  
00034 *  When DIRECT = 'F' (Forward sequence), then
00035 *  
00036 *     P = P(z-1) * ... * P(2) * P(1)
00037 *  
00038 *  and when DIRECT = 'B' (Backward sequence), then
00039 *  
00040 *     P = P(1) * P(2) * ... * P(z-1)
00041 *  
00042 *  where P(k) is a plane rotation matrix defined by the 2-by-2 rotation
00043 *  
00044 *     R(k) = (  c(k)  s(k) )
00045 *          = ( -s(k)  c(k) ).
00046 *  
00047 *  When PIVOT = 'V' (Variable pivot), the rotation is performed
00048 *  for the plane (k,k+1), i.e., P(k) has the form
00049 *  
00050 *     P(k) = (  1                                            )
00051 *            (       ...                                     )
00052 *            (              1                                )
00053 *            (                   c(k)  s(k)                  )
00054 *            (                  -s(k)  c(k)                  )
00055 *            (                                1              )
00056 *            (                                     ...       )
00057 *            (                                            1  )
00058 *  
00059 *  where R(k) appears as a rank-2 modification to the identity matrix in
00060 *  rows and columns k and k+1.
00061 *  
00062 *  When PIVOT = 'T' (Top pivot), the rotation is performed for the
00063 *  plane (1,k+1), so P(k) has the form
00064 *  
00065 *     P(k) = (  c(k)                    s(k)                 )
00066 *            (         1                                     )
00067 *            (              ...                              )
00068 *            (                     1                         )
00069 *            ( -s(k)                    c(k)                 )
00070 *            (                                 1             )
00071 *            (                                      ...      )
00072 *            (                                             1 )
00073 *  
00074 *  where R(k) appears in rows and columns 1 and k+1.
00075 *  
00076 *  Similarly, when PIVOT = 'B' (Bottom pivot), the rotation is
00077 *  performed for the plane (k,z), giving P(k) the form
00078 *  
00079 *     P(k) = ( 1                                             )
00080 *            (      ...                                      )
00081 *            (             1                                 )
00082 *            (                  c(k)                    s(k) )
00083 *            (                         1                     )
00084 *            (                              ...              )
00085 *            (                                     1         )
00086 *            (                 -s(k)                    c(k) )
00087 *  
00088 *  where R(k) appears in rows and columns k and z.  The rotations are
00089 *  performed without ever forming P(k) explicitly.
00090 *
00091 *  Arguments
00092 *  =========
00093 *
00094 *  SIDE    (input) CHARACTER*1
00095 *          Specifies whether the plane rotation matrix P is applied to
00096 *          A on the left or the right.
00097 *          = 'L':  Left, compute A := P*A
00098 *          = 'R':  Right, compute A:= A*P**T
00099 *
00100 *  PIVOT   (input) CHARACTER*1
00101 *          Specifies the plane for which P(k) is a plane rotation
00102 *          matrix.
00103 *          = 'V':  Variable pivot, the plane (k,k+1)
00104 *          = 'T':  Top pivot, the plane (1,k+1)
00105 *          = 'B':  Bottom pivot, the plane (k,z)
00106 *
00107 *  DIRECT  (input) CHARACTER*1
00108 *          Specifies whether P is a forward or backward sequence of
00109 *          plane rotations.
00110 *          = 'F':  Forward, P = P(z-1)*...*P(2)*P(1)
00111 *          = 'B':  Backward, P = P(1)*P(2)*...*P(z-1)
00112 *
00113 *  M       (input) INTEGER
00114 *          The number of rows of the matrix A.  If m <= 1, an immediate
00115 *          return is effected.
00116 *
00117 *  N       (input) INTEGER
00118 *          The number of columns of the matrix A.  If n <= 1, an
00119 *          immediate return is effected.
00120 *
00121 *  C       (input) DOUBLE PRECISION array, dimension
00122 *                  (M-1) if SIDE = 'L'
00123 *                  (N-1) if SIDE = 'R'
00124 *          The cosines c(k) of the plane rotations.
00125 *
00126 *  S       (input) DOUBLE PRECISION array, dimension
00127 *                  (M-1) if SIDE = 'L'
00128 *                  (N-1) if SIDE = 'R'
00129 *          The sines s(k) of the plane rotations.  The 2-by-2 plane
00130 *          rotation part of the matrix P(k), R(k), has the form
00131 *          R(k) = (  c(k)  s(k) )
00132 *                 ( -s(k)  c(k) ).
00133 *
00134 *  A       (input/output) DOUBLE PRECISION array, dimension (LDA,N)
00135 *          The M-by-N matrix A.  On exit, A is overwritten by P*A if
00136 *          SIDE = 'R' or by A*P**T if SIDE = 'L'.
00137 *
00138 *  LDA     (input) INTEGER
00139 *          The leading dimension of the array A.  LDA >= max(1,M).
00140 *
00141 *  =====================================================================
00142 *
00143 *     .. Parameters ..
00144       DOUBLE PRECISION   ONE, ZERO
00145       PARAMETER          ( ONE = 1.0D+0, ZERO = 0.0D+0 )
00146 *     ..
00147 *     .. Local Scalars ..
00148       INTEGER            I, INFO, J
00149       DOUBLE PRECISION   CTEMP, STEMP, TEMP
00150 *     ..
00151 *     .. External Functions ..
00152       LOGICAL            LSAME
00153       EXTERNAL           LSAME
00154 *     ..
00155 *     .. External Subroutines ..
00156       EXTERNAL           XERBLA
00157 *     ..
00158 *     .. Intrinsic Functions ..
00159       INTRINSIC          MAX
00160 *     ..
00161 *     .. Executable Statements ..
00162 *
00163 *     Test the input parameters
00164 *
00165       INFO = 0
00166       IF( .NOT.( LSAME( SIDE, 'L' ) .OR. LSAME( SIDE, 'R' ) ) ) THEN
00167          INFO = 1
00168       ELSE IF( .NOT.( LSAME( PIVOT, 'V' ) .OR. LSAME( PIVOT,
00169      $         'T' ) .OR. LSAME( PIVOT, 'B' ) ) ) THEN
00170          INFO = 2
00171       ELSE IF( .NOT.( LSAME( DIRECT, 'F' ) .OR. LSAME( DIRECT, 'B' ) ) )
00172      $          THEN
00173          INFO = 3
00174       ELSE IF( M.LT.0 ) THEN
00175          INFO = 4
00176       ELSE IF( N.LT.0 ) THEN
00177          INFO = 5
00178       ELSE IF( LDA.LT.MAX( 1, M ) ) THEN
00179          INFO = 9
00180       END IF
00181       IF( INFO.NE.0 ) THEN
00182          CALL XERBLA( 'DLASR ', INFO )
00183          RETURN
00184       END IF
00185 *
00186 *     Quick return if possible
00187 *
00188       IF( ( M.EQ.0 ) .OR. ( N.EQ.0 ) )
00189      $   RETURN
00190       IF( LSAME( SIDE, 'L' ) ) THEN
00191 *
00192 *        Form  P * A
00193 *
00194          IF( LSAME( PIVOT, 'V' ) ) THEN
00195             IF( LSAME( DIRECT, 'F' ) ) THEN
00196                DO 20 J = 1, M - 1
00197                   CTEMP = C( J )
00198                   STEMP = S( J )
00199                   IF( ( CTEMP.NE.ONE ) .OR. ( STEMP.NE.ZERO ) ) THEN
00200                      DO 10 I = 1, N
00201                         TEMP = A( J+1, I )
00202                         A( J+1, I ) = CTEMP*TEMP - STEMP*A( J, I )
00203                         A( J, I ) = STEMP*TEMP + CTEMP*A( J, I )
00204    10                CONTINUE
00205                   END IF
00206    20          CONTINUE
00207             ELSE IF( LSAME( DIRECT, 'B' ) ) THEN
00208                DO 40 J = M - 1, 1, -1
00209                   CTEMP = C( J )
00210                   STEMP = S( J )
00211                   IF( ( CTEMP.NE.ONE ) .OR. ( STEMP.NE.ZERO ) ) THEN
00212                      DO 30 I = 1, N
00213                         TEMP = A( J+1, I )
00214                         A( J+1, I ) = CTEMP*TEMP - STEMP*A( J, I )
00215                         A( J, I ) = STEMP*TEMP + CTEMP*A( J, I )
00216    30                CONTINUE
00217                   END IF
00218    40          CONTINUE
00219             END IF
00220          ELSE IF( LSAME( PIVOT, 'T' ) ) THEN
00221             IF( LSAME( DIRECT, 'F' ) ) THEN
00222                DO 60 J = 2, M
00223                   CTEMP = C( J-1 )
00224                   STEMP = S( J-1 )
00225                   IF( ( CTEMP.NE.ONE ) .OR. ( STEMP.NE.ZERO ) ) THEN
00226                      DO 50 I = 1, N
00227                         TEMP = A( J, I )
00228                         A( J, I ) = CTEMP*TEMP - STEMP*A( 1, I )
00229                         A( 1, I ) = STEMP*TEMP + CTEMP*A( 1, I )
00230    50                CONTINUE
00231                   END IF
00232    60          CONTINUE
00233             ELSE IF( LSAME( DIRECT, 'B' ) ) THEN
00234                DO 80 J = M, 2, -1
00235                   CTEMP = C( J-1 )
00236                   STEMP = S( J-1 )
00237                   IF( ( CTEMP.NE.ONE ) .OR. ( STEMP.NE.ZERO ) ) THEN
00238                      DO 70 I = 1, N
00239                         TEMP = A( J, I )
00240                         A( J, I ) = CTEMP*TEMP - STEMP*A( 1, I )
00241                         A( 1, I ) = STEMP*TEMP + CTEMP*A( 1, I )
00242    70                CONTINUE
00243                   END IF
00244    80          CONTINUE
00245             END IF
00246          ELSE IF( LSAME( PIVOT, 'B' ) ) THEN
00247             IF( LSAME( DIRECT, 'F' ) ) THEN
00248                DO 100 J = 1, M - 1
00249                   CTEMP = C( J )
00250                   STEMP = S( J )
00251                   IF( ( CTEMP.NE.ONE ) .OR. ( STEMP.NE.ZERO ) ) THEN
00252                      DO 90 I = 1, N
00253                         TEMP = A( J, I )
00254                         A( J, I ) = STEMP*A( M, I ) + CTEMP*TEMP
00255                         A( M, I ) = CTEMP*A( M, I ) - STEMP*TEMP
00256    90                CONTINUE
00257                   END IF
00258   100          CONTINUE
00259             ELSE IF( LSAME( DIRECT, 'B' ) ) THEN
00260                DO 120 J = M - 1, 1, -1
00261                   CTEMP = C( J )
00262                   STEMP = S( J )
00263                   IF( ( CTEMP.NE.ONE ) .OR. ( STEMP.NE.ZERO ) ) THEN
00264                      DO 110 I = 1, N
00265                         TEMP = A( J, I )
00266                         A( J, I ) = STEMP*A( M, I ) + CTEMP*TEMP
00267                         A( M, I ) = CTEMP*A( M, I ) - STEMP*TEMP
00268   110                CONTINUE
00269                   END IF
00270   120          CONTINUE
00271             END IF
00272          END IF
00273       ELSE IF( LSAME( SIDE, 'R' ) ) THEN
00274 *
00275 *        Form A * P'
00276 *
00277          IF( LSAME( PIVOT, 'V' ) ) THEN
00278             IF( LSAME( DIRECT, 'F' ) ) THEN
00279                DO 140 J = 1, N - 1
00280                   CTEMP = C( J )
00281                   STEMP = S( J )
00282                   IF( ( CTEMP.NE.ONE ) .OR. ( STEMP.NE.ZERO ) ) THEN
00283                      DO 130 I = 1, M
00284                         TEMP = A( I, J+1 )
00285                         A( I, J+1 ) = CTEMP*TEMP - STEMP*A( I, J )
00286                         A( I, J ) = STEMP*TEMP + CTEMP*A( I, J )
00287   130                CONTINUE
00288                   END IF
00289   140          CONTINUE
00290             ELSE IF( LSAME( DIRECT, 'B' ) ) THEN
00291                DO 160 J = N - 1, 1, -1
00292                   CTEMP = C( J )
00293                   STEMP = S( J )
00294                   IF( ( CTEMP.NE.ONE ) .OR. ( STEMP.NE.ZERO ) ) THEN
00295                      DO 150 I = 1, M
00296                         TEMP = A( I, J+1 )
00297                         A( I, J+1 ) = CTEMP*TEMP - STEMP*A( I, J )
00298                         A( I, J ) = STEMP*TEMP + CTEMP*A( I, J )
00299   150                CONTINUE
00300                   END IF
00301   160          CONTINUE
00302             END IF
00303          ELSE IF( LSAME( PIVOT, 'T' ) ) THEN
00304             IF( LSAME( DIRECT, 'F' ) ) THEN
00305                DO 180 J = 2, N
00306                   CTEMP = C( J-1 )
00307                   STEMP = S( J-1 )
00308                   IF( ( CTEMP.NE.ONE ) .OR. ( STEMP.NE.ZERO ) ) THEN
00309                      DO 170 I = 1, M
00310                         TEMP = A( I, J )
00311                         A( I, J ) = CTEMP*TEMP - STEMP*A( I, 1 )
00312                         A( I, 1 ) = STEMP*TEMP + CTEMP*A( I, 1 )
00313   170                CONTINUE
00314                   END IF
00315   180          CONTINUE
00316             ELSE IF( LSAME( DIRECT, 'B' ) ) THEN
00317                DO 200 J = N, 2, -1
00318                   CTEMP = C( J-1 )
00319                   STEMP = S( J-1 )
00320                   IF( ( CTEMP.NE.ONE ) .OR. ( STEMP.NE.ZERO ) ) THEN
00321                      DO 190 I = 1, M
00322                         TEMP = A( I, J )
00323                         A( I, J ) = CTEMP*TEMP - STEMP*A( I, 1 )
00324                         A( I, 1 ) = STEMP*TEMP + CTEMP*A( I, 1 )
00325   190                CONTINUE
00326                   END IF
00327   200          CONTINUE
00328             END IF
00329          ELSE IF( LSAME( PIVOT, 'B' ) ) THEN
00330             IF( LSAME( DIRECT, 'F' ) ) THEN
00331                DO 220 J = 1, N - 1
00332                   CTEMP = C( J )
00333                   STEMP = S( J )
00334                   IF( ( CTEMP.NE.ONE ) .OR. ( STEMP.NE.ZERO ) ) THEN
00335                      DO 210 I = 1, M
00336                         TEMP = A( I, J )
00337                         A( I, J ) = STEMP*A( I, N ) + CTEMP*TEMP
00338                         A( I, N ) = CTEMP*A( I, N ) - STEMP*TEMP
00339   210                CONTINUE
00340                   END IF
00341   220          CONTINUE
00342             ELSE IF( LSAME( DIRECT, 'B' ) ) THEN
00343                DO 240 J = N - 1, 1, -1
00344                   CTEMP = C( J )
00345                   STEMP = S( J )
00346                   IF( ( CTEMP.NE.ONE ) .OR. ( STEMP.NE.ZERO ) ) THEN
00347                      DO 230 I = 1, M
00348                         TEMP = A( I, J )
00349                         A( I, J ) = STEMP*A( I, N ) + CTEMP*TEMP
00350                         A( I, N ) = CTEMP*A( I, N ) - STEMP*TEMP
00351   230                CONTINUE
00352                   END IF
00353   240          CONTINUE
00354             END IF
00355          END IF
00356       END IF
00357 *
00358       RETURN
00359 *
00360 *     End of DLASR
00361 *
00362       END