LAPACK 3.3.0
|
00001 SUBROUTINE DLARF( SIDE, M, N, V, INCV, TAU, C, LDC, WORK ) 00002 IMPLICIT NONE 00003 * 00004 * -- LAPACK auxiliary routine (version 3.2) -- 00005 * -- LAPACK is a software package provided by Univ. of Tennessee, -- 00006 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 00007 * November 2006 00008 * 00009 * .. Scalar Arguments .. 00010 CHARACTER SIDE 00011 INTEGER INCV, LDC, M, N 00012 DOUBLE PRECISION TAU 00013 * .. 00014 * .. Array Arguments .. 00015 DOUBLE PRECISION C( LDC, * ), V( * ), WORK( * ) 00016 * .. 00017 * 00018 * Purpose 00019 * ======= 00020 * 00021 * DLARF applies a real elementary reflector H to a real m by n matrix 00022 * C, from either the left or the right. H is represented in the form 00023 * 00024 * H = I - tau * v * v' 00025 * 00026 * where tau is a real scalar and v is a real vector. 00027 * 00028 * If tau = 0, then H is taken to be the unit matrix. 00029 * 00030 * Arguments 00031 * ========= 00032 * 00033 * SIDE (input) CHARACTER*1 00034 * = 'L': form H * C 00035 * = 'R': form C * H 00036 * 00037 * M (input) INTEGER 00038 * The number of rows of the matrix C. 00039 * 00040 * N (input) INTEGER 00041 * The number of columns of the matrix C. 00042 * 00043 * V (input) DOUBLE PRECISION array, dimension 00044 * (1 + (M-1)*abs(INCV)) if SIDE = 'L' 00045 * or (1 + (N-1)*abs(INCV)) if SIDE = 'R' 00046 * The vector v in the representation of H. V is not used if 00047 * TAU = 0. 00048 * 00049 * INCV (input) INTEGER 00050 * The increment between elements of v. INCV <> 0. 00051 * 00052 * TAU (input) DOUBLE PRECISION 00053 * The value tau in the representation of H. 00054 * 00055 * C (input/output) DOUBLE PRECISION array, dimension (LDC,N) 00056 * On entry, the m by n matrix C. 00057 * On exit, C is overwritten by the matrix H * C if SIDE = 'L', 00058 * or C * H if SIDE = 'R'. 00059 * 00060 * LDC (input) INTEGER 00061 * The leading dimension of the array C. LDC >= max(1,M). 00062 * 00063 * WORK (workspace) DOUBLE PRECISION array, dimension 00064 * (N) if SIDE = 'L' 00065 * or (M) if SIDE = 'R' 00066 * 00067 * ===================================================================== 00068 * 00069 * .. Parameters .. 00070 DOUBLE PRECISION ONE, ZERO 00071 PARAMETER ( ONE = 1.0D+0, ZERO = 0.0D+0 ) 00072 * .. 00073 * .. Local Scalars .. 00074 LOGICAL APPLYLEFT 00075 INTEGER I, LASTV, LASTC 00076 * .. 00077 * .. External Subroutines .. 00078 EXTERNAL DGEMV, DGER 00079 * .. 00080 * .. External Functions .. 00081 LOGICAL LSAME 00082 INTEGER ILADLR, ILADLC 00083 EXTERNAL LSAME, ILADLR, ILADLC 00084 * .. 00085 * .. Executable Statements .. 00086 * 00087 APPLYLEFT = LSAME( SIDE, 'L' ) 00088 LASTV = 0 00089 LASTC = 0 00090 IF( TAU.NE.ZERO ) THEN 00091 ! Set up variables for scanning V. LASTV begins pointing to the end 00092 ! of V. 00093 IF( APPLYLEFT ) THEN 00094 LASTV = M 00095 ELSE 00096 LASTV = N 00097 END IF 00098 IF( INCV.GT.0 ) THEN 00099 I = 1 + (LASTV-1) * INCV 00100 ELSE 00101 I = 1 00102 END IF 00103 ! Look for the last non-zero row in V. 00104 DO WHILE( LASTV.GT.0 .AND. V( I ).EQ.ZERO ) 00105 LASTV = LASTV - 1 00106 I = I - INCV 00107 END DO 00108 IF( APPLYLEFT ) THEN 00109 ! Scan for the last non-zero column in C(1:lastv,:). 00110 LASTC = ILADLC(LASTV, N, C, LDC) 00111 ELSE 00112 ! Scan for the last non-zero row in C(:,1:lastv). 00113 LASTC = ILADLR(M, LASTV, C, LDC) 00114 END IF 00115 END IF 00116 ! Note that lastc.eq.0 renders the BLAS operations null; no special 00117 ! case is needed at this level. 00118 IF( APPLYLEFT ) THEN 00119 * 00120 * Form H * C 00121 * 00122 IF( LASTV.GT.0 ) THEN 00123 * 00124 * w(1:lastc,1) := C(1:lastv,1:lastc)' * v(1:lastv,1) 00125 * 00126 CALL DGEMV( 'Transpose', LASTV, LASTC, ONE, C, LDC, V, INCV, 00127 $ ZERO, WORK, 1 ) 00128 * 00129 * C(1:lastv,1:lastc) := C(...) - v(1:lastv,1) * w(1:lastc,1)' 00130 * 00131 CALL DGER( LASTV, LASTC, -TAU, V, INCV, WORK, 1, C, LDC ) 00132 END IF 00133 ELSE 00134 * 00135 * Form C * H 00136 * 00137 IF( LASTV.GT.0 ) THEN 00138 * 00139 * w(1:lastc,1) := C(1:lastc,1:lastv) * v(1:lastv,1) 00140 * 00141 CALL DGEMV( 'No transpose', LASTC, LASTV, ONE, C, LDC, 00142 $ V, INCV, ZERO, WORK, 1 ) 00143 * 00144 * C(1:lastc,1:lastv) := C(...) - w(1:lastc,1) * v(1:lastv,1)' 00145 * 00146 CALL DGER( LASTC, LASTV, -TAU, WORK, 1, V, INCV, C, LDC ) 00147 END IF 00148 END IF 00149 RETURN 00150 * 00151 * End of DLARF 00152 * 00153 END