File:  [local] / rpl / lapack / lapack / dlasr.f
Revision 1.1: download - view: text, annotated - select for diffs - revision graph
Tue Jan 26 15:22:46 2010 UTC (14 years, 3 months ago) by bertrand
Branches: MAIN
CVS tags: HEAD
Initial revision

    1:       SUBROUTINE DLASR( SIDE, PIVOT, DIRECT, M, N, C, S, A, LDA )
    2: *
    3: *  -- LAPACK auxiliary routine (version 3.2) --
    4: *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
    5: *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
    6: *     November 2006
    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 applies a sequence of plane rotations to a real matrix A,
   20: *  from either the left or the right.
   21: *  
   22: *  When SIDE = 'L', the transformation takes the form
   23: *  
   24: *     A := P*A
   25: *  
   26: *  and when SIDE = 'R', the transformation takes the form
   27: *  
   28: *     A := A*P**T
   29: *  
   30: *  where P is an orthogonal matrix consisting of a sequence of z plane
   31: *  rotations, with z = M when SIDE = 'L' and z = N when SIDE = 'R',
   32: *  and P**T is the transpose of P.
   33: *  
   34: *  When DIRECT = 'F' (Forward sequence), then
   35: *  
   36: *     P = P(z-1) * ... * P(2) * P(1)
   37: *  
   38: *  and when DIRECT = 'B' (Backward sequence), then
   39: *  
   40: *     P = P(1) * P(2) * ... * P(z-1)
   41: *  
   42: *  where P(k) is a plane rotation matrix defined by the 2-by-2 rotation
   43: *  
   44: *     R(k) = (  c(k)  s(k) )
   45: *          = ( -s(k)  c(k) ).
   46: *  
   47: *  When PIVOT = 'V' (Variable pivot), the rotation is performed
   48: *  for the plane (k,k+1), i.e., P(k) has the form
   49: *  
   50: *     P(k) = (  1                                            )
   51: *            (       ...                                     )
   52: *            (              1                                )
   53: *            (                   c(k)  s(k)                  )
   54: *            (                  -s(k)  c(k)                  )
   55: *            (                                1              )
   56: *            (                                     ...       )
   57: *            (                                            1  )
   58: *  
   59: *  where R(k) appears as a rank-2 modification to the identity matrix in
   60: *  rows and columns k and k+1.
   61: *  
   62: *  When PIVOT = 'T' (Top pivot), the rotation is performed for the
   63: *  plane (1,k+1), so P(k) has the form
   64: *  
   65: *     P(k) = (  c(k)                    s(k)                 )
   66: *            (         1                                     )
   67: *            (              ...                              )
   68: *            (                     1                         )
   69: *            ( -s(k)                    c(k)                 )
   70: *            (                                 1             )
   71: *            (                                      ...      )
   72: *            (                                             1 )
   73: *  
   74: *  where R(k) appears in rows and columns 1 and k+1.
   75: *  
   76: *  Similarly, when PIVOT = 'B' (Bottom pivot), the rotation is
   77: *  performed for the plane (k,z), giving P(k) the form
   78: *  
   79: *     P(k) = ( 1                                             )
   80: *            (      ...                                      )
   81: *            (             1                                 )
   82: *            (                  c(k)                    s(k) )
   83: *            (                         1                     )
   84: *            (                              ...              )
   85: *            (                                     1         )
   86: *            (                 -s(k)                    c(k) )
   87: *  
   88: *  where R(k) appears in rows and columns k and z.  The rotations are
   89: *  performed without ever forming P(k) explicitly.
   90: *
   91: *  Arguments
   92: *  =========
   93: *
   94: *  SIDE    (input) CHARACTER*1
   95: *          Specifies whether the plane rotation matrix P is applied to
   96: *          A on the left or the right.
   97: *          = 'L':  Left, compute A := P*A
   98: *          = 'R':  Right, compute A:= A*P**T
   99: *
  100: *  PIVOT   (input) CHARACTER*1
  101: *          Specifies the plane for which P(k) is a plane rotation
  102: *          matrix.
  103: *          = 'V':  Variable pivot, the plane (k,k+1)
  104: *          = 'T':  Top pivot, the plane (1,k+1)
  105: *          = 'B':  Bottom pivot, the plane (k,z)
  106: *
  107: *  DIRECT  (input) CHARACTER*1
  108: *          Specifies whether P is a forward or backward sequence of
  109: *          plane rotations.
  110: *          = 'F':  Forward, P = P(z-1)*...*P(2)*P(1)
  111: *          = 'B':  Backward, P = P(1)*P(2)*...*P(z-1)
  112: *
  113: *  M       (input) INTEGER
  114: *          The number of rows of the matrix A.  If m <= 1, an immediate
  115: *          return is effected.
  116: *
  117: *  N       (input) INTEGER
  118: *          The number of columns of the matrix A.  If n <= 1, an
  119: *          immediate return is effected.
  120: *
  121: *  C       (input) DOUBLE PRECISION array, dimension
  122: *                  (M-1) if SIDE = 'L'
  123: *                  (N-1) if SIDE = 'R'
  124: *          The cosines c(k) of the plane rotations.
  125: *
  126: *  S       (input) DOUBLE PRECISION array, dimension
  127: *                  (M-1) if SIDE = 'L'
  128: *                  (N-1) if SIDE = 'R'
  129: *          The sines s(k) of the plane rotations.  The 2-by-2 plane
  130: *          rotation part of the matrix P(k), R(k), has the form
  131: *          R(k) = (  c(k)  s(k) )
  132: *                 ( -s(k)  c(k) ).
  133: *
  134: *  A       (input/output) DOUBLE PRECISION array, dimension (LDA,N)
  135: *          The M-by-N matrix A.  On exit, A is overwritten by P*A if
  136: *          SIDE = 'R' or by A*P**T if SIDE = 'L'.
  137: *
  138: *  LDA     (input) INTEGER
  139: *          The leading dimension of the array A.  LDA >= max(1,M).
  140: *
  141: *  =====================================================================
  142: *
  143: *     .. Parameters ..
  144:       DOUBLE PRECISION   ONE, ZERO
  145:       PARAMETER          ( ONE = 1.0D+0, ZERO = 0.0D+0 )
  146: *     ..
  147: *     .. Local Scalars ..
  148:       INTEGER            I, INFO, J
  149:       DOUBLE PRECISION   CTEMP, STEMP, TEMP
  150: *     ..
  151: *     .. External Functions ..
  152:       LOGICAL            LSAME
  153:       EXTERNAL           LSAME
  154: *     ..
  155: *     .. External Subroutines ..
  156:       EXTERNAL           XERBLA
  157: *     ..
  158: *     .. Intrinsic Functions ..
  159:       INTRINSIC          MAX
  160: *     ..
  161: *     .. Executable Statements ..
  162: *
  163: *     Test the input parameters
  164: *
  165:       INFO = 0
  166:       IF( .NOT.( LSAME( SIDE, 'L' ) .OR. LSAME( SIDE, 'R' ) ) ) THEN
  167:          INFO = 1
  168:       ELSE IF( .NOT.( LSAME( PIVOT, 'V' ) .OR. LSAME( PIVOT,
  169:      $         'T' ) .OR. LSAME( PIVOT, 'B' ) ) ) THEN
  170:          INFO = 2
  171:       ELSE IF( .NOT.( LSAME( DIRECT, 'F' ) .OR. LSAME( DIRECT, 'B' ) ) )
  172:      $          THEN
  173:          INFO = 3
  174:       ELSE IF( M.LT.0 ) THEN
  175:          INFO = 4
  176:       ELSE IF( N.LT.0 ) THEN
  177:          INFO = 5
  178:       ELSE IF( LDA.LT.MAX( 1, M ) ) THEN
  179:          INFO = 9
  180:       END IF
  181:       IF( INFO.NE.0 ) THEN
  182:          CALL XERBLA( 'DLASR ', INFO )
  183:          RETURN
  184:       END IF
  185: *
  186: *     Quick return if possible
  187: *
  188:       IF( ( M.EQ.0 ) .OR. ( N.EQ.0 ) )
  189:      $   RETURN
  190:       IF( LSAME( SIDE, 'L' ) ) THEN
  191: *
  192: *        Form  P * A
  193: *
  194:          IF( LSAME( PIVOT, 'V' ) ) THEN
  195:             IF( LSAME( DIRECT, 'F' ) ) THEN
  196:                DO 20 J = 1, M - 1
  197:                   CTEMP = C( J )
  198:                   STEMP = S( J )
  199:                   IF( ( CTEMP.NE.ONE ) .OR. ( STEMP.NE.ZERO ) ) THEN
  200:                      DO 10 I = 1, N
  201:                         TEMP = A( J+1, I )
  202:                         A( J+1, I ) = CTEMP*TEMP - STEMP*A( J, I )
  203:                         A( J, I ) = STEMP*TEMP + CTEMP*A( J, I )
  204:    10                CONTINUE
  205:                   END IF
  206:    20          CONTINUE
  207:             ELSE IF( LSAME( DIRECT, 'B' ) ) THEN
  208:                DO 40 J = M - 1, 1, -1
  209:                   CTEMP = C( J )
  210:                   STEMP = S( J )
  211:                   IF( ( CTEMP.NE.ONE ) .OR. ( STEMP.NE.ZERO ) ) THEN
  212:                      DO 30 I = 1, N
  213:                         TEMP = A( J+1, I )
  214:                         A( J+1, I ) = CTEMP*TEMP - STEMP*A( J, I )
  215:                         A( J, I ) = STEMP*TEMP + CTEMP*A( J, I )
  216:    30                CONTINUE
  217:                   END IF
  218:    40          CONTINUE
  219:             END IF
  220:          ELSE IF( LSAME( PIVOT, 'T' ) ) THEN
  221:             IF( LSAME( DIRECT, 'F' ) ) THEN
  222:                DO 60 J = 2, M
  223:                   CTEMP = C( J-1 )
  224:                   STEMP = S( J-1 )
  225:                   IF( ( CTEMP.NE.ONE ) .OR. ( STEMP.NE.ZERO ) ) THEN
  226:                      DO 50 I = 1, N
  227:                         TEMP = A( J, I )
  228:                         A( J, I ) = CTEMP*TEMP - STEMP*A( 1, I )
  229:                         A( 1, I ) = STEMP*TEMP + CTEMP*A( 1, I )
  230:    50                CONTINUE
  231:                   END IF
  232:    60          CONTINUE
  233:             ELSE IF( LSAME( DIRECT, 'B' ) ) THEN
  234:                DO 80 J = M, 2, -1
  235:                   CTEMP = C( J-1 )
  236:                   STEMP = S( J-1 )
  237:                   IF( ( CTEMP.NE.ONE ) .OR. ( STEMP.NE.ZERO ) ) THEN
  238:                      DO 70 I = 1, N
  239:                         TEMP = A( J, I )
  240:                         A( J, I ) = CTEMP*TEMP - STEMP*A( 1, I )
  241:                         A( 1, I ) = STEMP*TEMP + CTEMP*A( 1, I )
  242:    70                CONTINUE
  243:                   END IF
  244:    80          CONTINUE
  245:             END IF
  246:          ELSE IF( LSAME( PIVOT, 'B' ) ) THEN
  247:             IF( LSAME( DIRECT, 'F' ) ) THEN
  248:                DO 100 J = 1, M - 1
  249:                   CTEMP = C( J )
  250:                   STEMP = S( J )
  251:                   IF( ( CTEMP.NE.ONE ) .OR. ( STEMP.NE.ZERO ) ) THEN
  252:                      DO 90 I = 1, N
  253:                         TEMP = A( J, I )
  254:                         A( J, I ) = STEMP*A( M, I ) + CTEMP*TEMP
  255:                         A( M, I ) = CTEMP*A( M, I ) - STEMP*TEMP
  256:    90                CONTINUE
  257:                   END IF
  258:   100          CONTINUE
  259:             ELSE IF( LSAME( DIRECT, 'B' ) ) THEN
  260:                DO 120 J = M - 1, 1, -1
  261:                   CTEMP = C( J )
  262:                   STEMP = S( J )
  263:                   IF( ( CTEMP.NE.ONE ) .OR. ( STEMP.NE.ZERO ) ) THEN
  264:                      DO 110 I = 1, N
  265:                         TEMP = A( J, I )
  266:                         A( J, I ) = STEMP*A( M, I ) + CTEMP*TEMP
  267:                         A( M, I ) = CTEMP*A( M, I ) - STEMP*TEMP
  268:   110                CONTINUE
  269:                   END IF
  270:   120          CONTINUE
  271:             END IF
  272:          END IF
  273:       ELSE IF( LSAME( SIDE, 'R' ) ) THEN
  274: *
  275: *        Form A * P'
  276: *
  277:          IF( LSAME( PIVOT, 'V' ) ) THEN
  278:             IF( LSAME( DIRECT, 'F' ) ) THEN
  279:                DO 140 J = 1, N - 1
  280:                   CTEMP = C( J )
  281:                   STEMP = S( J )
  282:                   IF( ( CTEMP.NE.ONE ) .OR. ( STEMP.NE.ZERO ) ) THEN
  283:                      DO 130 I = 1, M
  284:                         TEMP = A( I, J+1 )
  285:                         A( I, J+1 ) = CTEMP*TEMP - STEMP*A( I, J )
  286:                         A( I, J ) = STEMP*TEMP + CTEMP*A( I, J )
  287:   130                CONTINUE
  288:                   END IF
  289:   140          CONTINUE
  290:             ELSE IF( LSAME( DIRECT, 'B' ) ) THEN
  291:                DO 160 J = N - 1, 1, -1
  292:                   CTEMP = C( J )
  293:                   STEMP = S( J )
  294:                   IF( ( CTEMP.NE.ONE ) .OR. ( STEMP.NE.ZERO ) ) THEN
  295:                      DO 150 I = 1, M
  296:                         TEMP = A( I, J+1 )
  297:                         A( I, J+1 ) = CTEMP*TEMP - STEMP*A( I, J )
  298:                         A( I, J ) = STEMP*TEMP + CTEMP*A( I, J )
  299:   150                CONTINUE
  300:                   END IF
  301:   160          CONTINUE
  302:             END IF
  303:          ELSE IF( LSAME( PIVOT, 'T' ) ) THEN
  304:             IF( LSAME( DIRECT, 'F' ) ) THEN
  305:                DO 180 J = 2, N
  306:                   CTEMP = C( J-1 )
  307:                   STEMP = S( J-1 )
  308:                   IF( ( CTEMP.NE.ONE ) .OR. ( STEMP.NE.ZERO ) ) THEN
  309:                      DO 170 I = 1, M
  310:                         TEMP = A( I, J )
  311:                         A( I, J ) = CTEMP*TEMP - STEMP*A( I, 1 )
  312:                         A( I, 1 ) = STEMP*TEMP + CTEMP*A( I, 1 )
  313:   170                CONTINUE
  314:                   END IF
  315:   180          CONTINUE
  316:             ELSE IF( LSAME( DIRECT, 'B' ) ) THEN
  317:                DO 200 J = N, 2, -1
  318:                   CTEMP = C( J-1 )
  319:                   STEMP = S( J-1 )
  320:                   IF( ( CTEMP.NE.ONE ) .OR. ( STEMP.NE.ZERO ) ) THEN
  321:                      DO 190 I = 1, M
  322:                         TEMP = A( I, J )
  323:                         A( I, J ) = CTEMP*TEMP - STEMP*A( I, 1 )
  324:                         A( I, 1 ) = STEMP*TEMP + CTEMP*A( I, 1 )
  325:   190                CONTINUE
  326:                   END IF
  327:   200          CONTINUE
  328:             END IF
  329:          ELSE IF( LSAME( PIVOT, 'B' ) ) THEN
  330:             IF( LSAME( DIRECT, 'F' ) ) THEN
  331:                DO 220 J = 1, N - 1
  332:                   CTEMP = C( J )
  333:                   STEMP = S( J )
  334:                   IF( ( CTEMP.NE.ONE ) .OR. ( STEMP.NE.ZERO ) ) THEN
  335:                      DO 210 I = 1, M
  336:                         TEMP = A( I, J )
  337:                         A( I, J ) = STEMP*A( I, N ) + CTEMP*TEMP
  338:                         A( I, N ) = CTEMP*A( I, N ) - STEMP*TEMP
  339:   210                CONTINUE
  340:                   END IF
  341:   220          CONTINUE
  342:             ELSE IF( LSAME( DIRECT, 'B' ) ) THEN
  343:                DO 240 J = N - 1, 1, -1
  344:                   CTEMP = C( J )
  345:                   STEMP = S( J )
  346:                   IF( ( CTEMP.NE.ONE ) .OR. ( STEMP.NE.ZERO ) ) THEN
  347:                      DO 230 I = 1, M
  348:                         TEMP = A( I, J )
  349:                         A( I, J ) = STEMP*A( I, N ) + CTEMP*TEMP
  350:                         A( I, N ) = CTEMP*A( I, N ) - STEMP*TEMP
  351:   230                CONTINUE
  352:                   END IF
  353:   240          CONTINUE
  354:             END IF
  355:          END IF
  356:       END IF
  357: *
  358:       RETURN
  359: *
  360: *     End of DLASR
  361: *
  362:       END

CVSweb interface <joel.bertrand@systella.fr>