|
| 1 | +*> \brief \b DLARF1L applies an elementary reflector to a general rectangular |
| 2 | +* matrix assuming v(lastv) = 1 where lastv is the last non-zero |
| 3 | +* element |
| 4 | +* |
| 5 | +* =========== DOCUMENTATION =========== |
| 6 | +* |
| 7 | +* Online html documentation available at |
| 8 | +* http://www.netlib.org/lapack/explore-html/ |
| 9 | +* |
| 10 | +*> \htmlonly |
| 11 | +*> Download DLARF + dependencies |
| 12 | +*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dlarf.f"> |
| 13 | +*> [TGZ]</a> |
| 14 | +*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dlarf.f"> |
| 15 | +*> [ZIP]</a> |
| 16 | +*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dlarf.f"> |
| 17 | +*> [TXT]</a> |
| 18 | +*> \endhtmlonly |
| 19 | +* |
| 20 | +* Definition: |
| 21 | +* =========== |
| 22 | +* |
| 23 | +* SUBROUTINE DLARF1L( SIDE, M, N, V, INCV, TAU, C, LDC, WORK ) |
| 24 | +* |
| 25 | +* .. Scalar Arguments .. |
| 26 | +* CHARACTER SIDE |
| 27 | +* INTEGER INCV, LDC, M, N |
| 28 | +* DOUBLE PRECISION TAU |
| 29 | +* .. |
| 30 | +* .. Array Arguments .. |
| 31 | +* DOUBLE PRECISION C( LDC, * ), V( * ), WORK( * ) |
| 32 | +* .. |
| 33 | +* |
| 34 | +* |
| 35 | +*> \par Purpose: |
| 36 | +* ============= |
| 37 | +*> |
| 38 | +*> \verbatim |
| 39 | +*> |
| 40 | +*> DLARF1L applies a real elementary reflector H to a real m by n matrix |
| 41 | +*> C, from either the left or the right. H is represented in the form |
| 42 | +*> |
| 43 | +*> H = I - tau * v * v**T |
| 44 | +*> |
| 45 | +*> where tau is a real scalar and v is a real vector. |
| 46 | +*> |
| 47 | +*> If tau = 0, then H is taken to be the unit matrix. |
| 48 | +*> \endverbatim |
| 49 | +* |
| 50 | +* Arguments: |
| 51 | +* ========== |
| 52 | +* |
| 53 | +*> \param[in] SIDE |
| 54 | +*> \verbatim |
| 55 | +*> SIDE is CHARACTER*1 |
| 56 | +*> = 'L': form H * C |
| 57 | +*> = 'R': form C * H |
| 58 | +*> \endverbatim |
| 59 | +*> |
| 60 | +*> \param[in] M |
| 61 | +*> \verbatim |
| 62 | +*> M is INTEGER |
| 63 | +*> The number of rows of the matrix C. |
| 64 | +*> \endverbatim |
| 65 | +*> |
| 66 | +*> \param[in] N |
| 67 | +*> \verbatim |
| 68 | +*> N is INTEGER |
| 69 | +*> The number of columns of the matrix C. |
| 70 | +*> \endverbatim |
| 71 | +*> |
| 72 | +*> \param[in] V |
| 73 | +*> \verbatim |
| 74 | +*> V is DOUBLE PRECISION array, dimension |
| 75 | +*> (1 + (M-1)*abs(INCV)) if SIDE = 'L' |
| 76 | +*> or (1 + (N-1)*abs(INCV)) if SIDE = 'R' |
| 77 | +*> The vector v in the representation of H. V is not used if |
| 78 | +*> TAU = 0. |
| 79 | +*> \endverbatim |
| 80 | +*> |
| 81 | +*> \param[in] INCV |
| 82 | +*> \verbatim |
| 83 | +*> INCV is INTEGER |
| 84 | +*> The increment between elements of v. INCV <> 0. |
| 85 | +*> \endverbatim |
| 86 | +*> |
| 87 | +*> \param[in] TAU |
| 88 | +*> \verbatim |
| 89 | +*> TAU is DOUBLE PRECISION |
| 90 | +*> The value tau in the representation of H. |
| 91 | +*> \endverbatim |
| 92 | +*> |
| 93 | +*> \param[in,out] C |
| 94 | +*> \verbatim |
| 95 | +*> C is DOUBLE PRECISION array, dimension (LDC,N) |
| 96 | +*> On entry, the m by n matrix C. |
| 97 | +*> On exit, C is overwritten by the matrix H * C if SIDE = 'L', |
| 98 | +*> or C * H if SIDE = 'R'. |
| 99 | +*> \endverbatim |
| 100 | +*> |
| 101 | +*> \param[in] LDC |
| 102 | +*> \verbatim |
| 103 | +*> LDC is INTEGER |
| 104 | +*> The leading dimension of the array C. LDC >= max(1,M). |
| 105 | +*> \endverbatim |
| 106 | +*> |
| 107 | +*> \param[out] WORK |
| 108 | +*> \verbatim |
| 109 | +*> WORK is DOUBLE PRECISION array, dimension |
| 110 | +*> (N) if SIDE = 'L' |
| 111 | +*> or (M) if SIDE = 'R' |
| 112 | +*> \endverbatim |
| 113 | +* |
| 114 | +* Authors: |
| 115 | +* ======== |
| 116 | +* |
| 117 | +*> \author Univ. of Tennessee |
| 118 | +*> \author Univ. of California Berkeley |
| 119 | +*> \author Univ. of Colorado Denver |
| 120 | +*> \author NAG Ltd. |
| 121 | +* |
| 122 | +*> \ingroup larf |
| 123 | +* |
| 124 | +* ===================================================================== |
| 125 | + SUBROUTINE DLARF1L( SIDE, M, N, V, INCV, TAU, C, LDC, WORK ) |
| 126 | +* |
| 127 | +* -- LAPACK auxiliary routine -- |
| 128 | +* -- LAPACK is a software package provided by Univ. of Tennessee, -- |
| 129 | +* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- |
| 130 | +* |
| 131 | +* .. Scalar Arguments .. |
| 132 | + CHARACTER SIDE |
| 133 | + INTEGER INCV, LDC, M, N |
| 134 | + DOUBLE PRECISION TAU |
| 135 | +* .. |
| 136 | +* .. Array Arguments .. |
| 137 | + DOUBLE PRECISION C( LDC, * ), V( * ), WORK( * ) |
| 138 | +* .. |
| 139 | +* |
| 140 | +* ===================================================================== |
| 141 | +* |
| 142 | +* .. Parameters .. |
| 143 | + DOUBLE PRECISION ONE, ZERO |
| 144 | + PARAMETER ( ONE = 1.0D+0, ZERO = 0.0D+0 ) |
| 145 | + INTEGER IONE |
| 146 | + PARAMETER ( IONE = 1 ) |
| 147 | +* .. |
| 148 | +* .. Local Scalars .. |
| 149 | + LOGICAL APPLYLEFT |
| 150 | + INTEGER I, LASTV, LASTC, J |
| 151 | +* .. |
| 152 | +* .. External Subroutines .. |
| 153 | + EXTERNAL DGEMV, DGER |
| 154 | +* .. |
| 155 | +* .. External Functions .. |
| 156 | + LOGICAL LSAME |
| 157 | + INTEGER ILADLR, ILADLC |
| 158 | + EXTERNAL LSAME, ILADLR, ILADLC |
| 159 | +* .. |
| 160 | +* .. Executable Statements .. |
| 161 | +* |
| 162 | + APPLYLEFT = LSAME( SIDE, 'L' ) |
| 163 | + LASTV = 0 |
| 164 | + LASTC = 0 |
| 165 | + IF( TAU.NE.ZERO ) THEN |
| 166 | +! Set up variables for scanning V. LASTV begins pointing to the end |
| 167 | +! of V. |
| 168 | + IF( APPLYLEFT ) THEN |
| 169 | + LASTV = M |
| 170 | + ELSE |
| 171 | + LASTV = N |
| 172 | + END IF |
| 173 | + IF( INCV.GT.0 ) THEN |
| 174 | + I = 1 + (LASTV-1) * INCV |
| 175 | + ELSE |
| 176 | + I = 1 |
| 177 | + END IF |
| 178 | +! Look for the last non-zero row in V. |
| 179 | + DO WHILE( LASTV.GT.0 .AND. V( I ).EQ.ZERO ) |
| 180 | + LASTV = LASTV - 1 |
| 181 | + I = I - INCV |
| 182 | + END DO |
| 183 | + IF( APPLYLEFT ) THEN |
| 184 | +! Scan for the last non-zero column in C(1:lastv,:). |
| 185 | + LASTC = ILADLC(LASTV, N, C, LDC) |
| 186 | + ELSE |
| 187 | +! Scan for the last non-zero row in C(:,1:lastv). |
| 188 | + LASTC = ILADLR(M, LASTV, C, LDC) |
| 189 | + END IF |
| 190 | + END IF |
| 191 | +! Note that lastc.eq.0 renders the BLAS operations null; no special |
| 192 | +! case is needed at this level. |
| 193 | + IF( APPLYLEFT ) THEN |
| 194 | +* |
| 195 | +* Form H * C |
| 196 | +* |
| 197 | + IF( LASTV.GT.0 ) THEN |
| 198 | + ! Check if m = 1. This means v = 1, So we just need to compute |
| 199 | + ! C := HC = (1-\tau)C. |
| 200 | + IF( LASTV.EQ.1 ) THEN |
| 201 | + CALL DSCAL(LASTC, ONE - TAU, C, LDC) |
| 202 | + ELSE |
| 203 | +* |
| 204 | +* w(1:lastc,1) := C(1:lastv,1:lastc)**T * v(1:lastv,1) |
| 205 | +* |
| 206 | + ! w(1:lastc,1) := C(1:lastv-1,1:lastc)**T * v(1:lastv-1,1) |
| 207 | + CALL DGEMV( 'Transpose', LASTV-1, LASTC, ONE, C(1,1), |
| 208 | + $ LDC, V(1), INCV, ZERO, WORK, 1) |
| 209 | + ! w(1:lastc,1) += C(lastv,1:lastc)**T * v(lastv,1) = C(lastv,1:lastc)**T |
| 210 | + CALL DAXPY(LASTC, ONE, C(LASTV,1), LDC, WORK, 1) |
| 211 | +* |
| 212 | +* C(1:lastv,1:lastc) := C(...) - tau * v(1:lastv,1) * w(1:lastc,1)**T |
| 213 | +* |
| 214 | + ! C(lastv, 1:lastc) := C(...) - tau * v(lastv,1) * w(1:lastc,1)**T |
| 215 | + ! = C(...) - tau * w(1:lastc,1)**T |
| 216 | + CALL DAXPY(LASTC, -TAU, WORK, 1, C(LASTV,1), LDC) |
| 217 | + ! C(1:lastv-1,1:lastc) := C(...) - tau * v(1:lastv-1,1)*w(1:lastc,1)**T |
| 218 | + CALL DGER(LASTV-1, LASTC, -TAU, V(1), INCV, WORK, 1, |
| 219 | + $ C(1,1), LDC) |
| 220 | + END IF |
| 221 | + END IF |
| 222 | + ELSE |
| 223 | +* |
| 224 | +* Form C * H |
| 225 | +* |
| 226 | + IF( LASTV.GT.0 ) THEN |
| 227 | + ! Check if n = 1. This means v = 1, so we just need to compute |
| 228 | + ! C := CH = C(1-\tau). |
| 229 | + IF( LASTV.EQ.1 ) THEN |
| 230 | + CALL DSCAL(LASTC, ONE - TAU, C, 1) |
| 231 | + ELSE |
| 232 | +* |
| 233 | +* w(1:lastc,1) := C(1:lastc,1:lastv) * v(1:lastv,1) |
| 234 | +* |
| 235 | + ! w(1:lastc,1) := C(1:lastc,1:lastv-1) * v(1:lastv-1,1) |
| 236 | + CALL DGEMV( 'No transpose', LASTC, LASTV-1, ONE, |
| 237 | + $ C(1,1), LDC, V(1), INCV, ZERO, WORK, 1 ) |
| 238 | + ! w(1:lastc,1) += C(1:lastc,lastv) * v(lastv,1) = C(1:lastc,lastv) |
| 239 | + CALL DAXPY(LASTC, ONE, C(1,LASTV), 1, WORK, 1) |
| 240 | +* |
| 241 | +* C(1:lastc,1:lastv) := C(...) - tau * w(1:lastc,1) * v(1:lastv,1)**T |
| 242 | +* |
| 243 | + ! C(1:lastc,lastv) := C(...) - tau * w(1:lastc,1) * v(lastv,1)**T |
| 244 | + ! = C(...) - tau * w(1:lastc,1) |
| 245 | + CALL DAXPY(LASTC, -TAU, WORK, 1, C(1,LASTV), 1) |
| 246 | + ! C(1:lastc,1:lastv-1) := C(...) - tau * w(1:lastc,1) * v(1:lastv-1)**T |
| 247 | + CALL DGER( LASTC, LASTV-1, -TAU, WORK, 1, V(1), |
| 248 | + $ INCV, C(1,1), LDC ) |
| 249 | + END IF |
| 250 | + END IF |
| 251 | + END IF |
| 252 | + RETURN |
| 253 | +* |
| 254 | +* End of DLARF |
| 255 | +* |
| 256 | + END |
0 commit comments