Chimie Théorique » Opt'n Path

      SUBROUTINE DLARF( SIDE, M, N, V, INCV, TAU, C, LDC, WORK )

      IMPLICIT NONE

*

*  -- LAPACK auxiliary routine (version 3.2) --

*  -- LAPACK is a software package provided by Univ. of Tennessee,    --

*  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--

*     November 2006

*

*     .. Scalar Arguments ..

      CHARACTER          SIDE

      INTEGER            INCV, LDC, M, N

      DOUBLE PRECISION   TAU

*     ..

*     .. Array Arguments ..

      DOUBLE PRECISION   C( LDC, * ), V( * ), WORK( * )

*     ..

*

*  Purpose

*  =======

*

*  DLARF applies a real elementary reflector H to a real m by n matrix

*  C, from either the left or the right. H is represented in the form

*

*        H = I - tau * v * v'

*

*  where tau is a real scalar and v is a real vector.

*

*  If tau = 0, then H is taken to be the unit matrix.

*

*  Arguments

*  =========

*

*  SIDE    (input) CHARACTER*1

*          = 'L': form  H * C

*          = 'R': form  C * H

*

*  M       (input) INTEGER

*          The number of rows of the matrix C.

*

*  N       (input) INTEGER

*          The number of columns of the matrix C.

*

*  V       (input) DOUBLE PRECISION array, dimension

*                     (1 + (M-1)*abs(INCV)) if SIDE = 'L'

*                  or (1 + (N-1)*abs(INCV)) if SIDE = 'R'

*          The vector v in the representation of H. V is not used if

*          TAU = 0.

*

*  INCV    (input) INTEGER

*          The increment between elements of v. INCV <> 0.

*

*  TAU     (input) DOUBLE PRECISION

*          The value tau in the representation of H.

*

*  C       (input/output) DOUBLE PRECISION array, dimension (LDC,N)

*          On entry, the m by n matrix C.

*          On exit, C is overwritten by the matrix H * C if SIDE = 'L',

*          or C * H if SIDE = 'R'.

*

*  LDC     (input) INTEGER

*          The leading dimension of the array C. LDC >= max(1,M).

*

*  WORK    (workspace) DOUBLE PRECISION array, dimension

*                         (N) if SIDE = 'L'

*                      or (M) if SIDE = 'R'

*

*  =====================================================================

*

*     .. Parameters ..

      DOUBLE PRECISION   ONE, ZERO

      PARAMETER          ( ONE = 1.0D+0, ZERO = 0.0D+0 )

*     ..

*     .. Local Scalars ..

      LOGICAL            APPLYLEFT

      INTEGER            I, LASTV, LASTC

*     ..

*     .. External Subroutines ..

      EXTERNAL           DGEMV, DGER

*     ..

*     .. External Functions ..

      LOGICAL            LSAME

      INTEGER            ILADLR, ILADLC

      EXTERNAL           LSAME, ILADLR, ILADLC

*     ..

*     .. Executable Statements ..

*

      APPLYLEFT = LSAME( SIDE, 'L' )

      LASTV = 0

      LASTC = 0

      IF( TAU.NE.ZERO ) THEN

!     Set up variables for scanning V.  LASTV begins pointing to the end

!     of V.

         IF( APPLYLEFT ) THEN

            LASTV = M

         ELSE

            LASTV = N

         END IF

         IF( INCV.GT.0 ) THEN

            I = 1 + (LASTV-1) * INCV

         ELSE

            I = 1

         END IF

!     Look for the last non-zero row in V.

         DO WHILE( LASTV.GT.0 .AND. V( I ).EQ.ZERO )

            LASTV = LASTV - 1

            I = I - INCV

         END DO

         IF( APPLYLEFT ) THEN

!     Scan for the last non-zero column in C(1:lastv,:).

            LASTC = ILADLC(LASTV, N, C, LDC)

         ELSE

!     Scan for the last non-zero row in C(:,1:lastv).

            LASTC = ILADLR(M, LASTV, C, LDC)

         END IF

      END IF

!     Note that lastc.eq.0 renders the BLAS operations null; no special

!     case is needed at this level.

      IF( APPLYLEFT ) THEN

*

*        Form  H * C

*

         IF( LASTV.GT.0 ) THEN

*

*           w(1:lastc,1) := C(1:lastv,1:lastc)' * v(1:lastv,1)

*

            CALL DGEMV( 'Transpose', LASTV, LASTC, ONE, C, LDC, V, INCV,

     $           ZERO, WORK, 1 )

*

*           C(1:lastv,1:lastc) := C(...) - v(1:lastv,1) * w(1:lastc,1)'

*

            CALL DGER( LASTV, LASTC, -TAU, V, INCV, WORK, 1, C, LDC )

         END IF

      ELSE

*

*        Form  C * H

*

         IF( LASTV.GT.0 ) THEN

*

*           w(1:lastc,1) := C(1:lastc,1:lastv) * v(1:lastv,1)

*

            CALL DGEMV( 'No transpose', LASTC, LASTV, ONE, C, LDC,

     $           V, INCV, ZERO, WORK, 1 )

*

*           C(1:lastc,1:lastv) := C(...) - w(1:lastc,1) * v(1:lastv,1)'

*

            CALL DGER( LASTC, LASTV, -TAU, WORK, 1, V, INCV, C, LDC )

         END IF

      END IF

      RETURN

*

*     End of DLARF

*

      END