vendor/gonum.org/v1/gonum/lapack/internal/testdata/netlib/dlahr2.f

   1 *> \brief \b DLAHR2 reduces the specified number of first columns of a general rectangular matrix A so that elements below the specified subdiagonal are zero, and returns auxiliary matrices which are needed to apply the transformation to the unreduced part of A.
   2 *
   3 *  =========== DOCUMENTATION ===========
   4 *
   5 * Online html documentation available at
   6 *            http://www.netlib.org/lapack/explore-html/
   7 *
   8 *> \htmlonly
   9 *> Download DLAHR2 + dependencies
  10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dlahr2.f">
  11 *> [TGZ]</a>
  12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dlahr2.f">
  13 *> [ZIP]</a>
  14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dlahr2.f">
  15 *> [TXT]</a>
  16 *> \endhtmlonly
  17 *
  18 *  Definition:
  19 *  ===========
  20 *
  21 *       SUBROUTINE DLAHR2( N, K, NB, A, LDA, TAU, T, LDT, Y, LDY )
  22 *
  23 *       .. Scalar Arguments ..
  24 *       INTEGER            K, LDA, LDT, LDY, N, NB
  25 *       ..
  26 *       .. Array Arguments ..
  27 *       DOUBLE PRECISION  A( LDA, * ), T( LDT, NB ), TAU( NB ),
  28 *      $                   Y( LDY, NB )
  29 *       ..
  30 *
  31 *
  32 *> \par Purpose:
  33 *  =============
  34 *>
  35 *> \verbatim
  36 *>
  37 *> DLAHR2 reduces the first NB columns of A real general n-BY-(n-k+1)
  38 *> matrix A so that elements below the k-th subdiagonal are zero. The
  39 *> reduction is performed by an orthogonal similarity transformation
  40 *> Q**T * A * Q. The routine returns the matrices V and T which determine
  41 *> Q as a block reflector I - V*T*V**T, and also the matrix Y = A * V * T.
  42 *>
  43 *> This is an auxiliary routine called by DGEHRD.
  44 *> \endverbatim
  45 *
  46 *  Arguments:
  47 *  ==========
  48 *
  49 *> \param[in] N
  50 *> \verbatim
  51 *>          N is INTEGER
  52 *>          The order of the matrix A.
  53 *> \endverbatim
  54 *>
  55 *> \param[in] K
  56 *> \verbatim
  57 *>          K is INTEGER
  58 *>          The offset for the reduction. Elements below the k-th
  59 *>          subdiagonal in the first NB columns are reduced to zero.
  60 *>          K < N.
  61 *> \endverbatim
  62 *>
  63 *> \param[in] NB
  64 *> \verbatim
  65 *>          NB is INTEGER
  66 *>          The number of columns to be reduced.
  67 *> \endverbatim
  68 *>
  69 *> \param[in,out] A
  70 *> \verbatim
  71 *>          A is DOUBLE PRECISION array, dimension (LDA,N-K+1)
  72 *>          On entry, the n-by-(n-k+1) general matrix A.
  73 *>          On exit, the elements on and above the k-th subdiagonal in
  74 *>          the first NB columns are overwritten with the corresponding
  75 *>          elements of the reduced matrix; the elements below the k-th
  76 *>          subdiagonal, with the array TAU, represent the matrix Q as a
  77 *>          product of elementary reflectors. The other columns of A are
  78 *>          unchanged. See Further Details.
  79 *> \endverbatim
  80 *>
  81 *> \param[in] LDA
  82 *> \verbatim
  83 *>          LDA is INTEGER
  84 *>          The leading dimension of the array A.  LDA >= max(1,N).
  85 *> \endverbatim
  86 *>
  87 *> \param[out] TAU
  88 *> \verbatim
  89 *>          TAU is DOUBLE PRECISION array, dimension (NB)
  90 *>          The scalar factors of the elementary reflectors. See Further
  91 *>          Details.
  92 *> \endverbatim
  93 *>
  94 *> \param[out] T
  95 *> \verbatim
  96 *>          T is DOUBLE PRECISION array, dimension (LDT,NB)
  97 *>          The upper triangular matrix T.
  98 *> \endverbatim
  99 *>
 100 *> \param[in] LDT
 101 *> \verbatim
 102 *>          LDT is INTEGER
 103 *>          The leading dimension of the array T.  LDT >= NB.
 104 *> \endverbatim
 105 *>
 106 *> \param[out] Y
 107 *> \verbatim
 108 *>          Y is DOUBLE PRECISION array, dimension (LDY,NB)
 109 *>          The n-by-nb matrix Y.
 110 *> \endverbatim
 111 *>
 112 *> \param[in] LDY
 113 *> \verbatim
 114 *>          LDY is INTEGER
 115 *>          The leading dimension of the array Y. LDY >= N.
 116 *> \endverbatim
 117 *
 118 *  Authors:
 119 *  ========
 120 *
 121 *> \author Univ. of Tennessee
 122 *> \author Univ. of California Berkeley
 123 *> \author Univ. of Colorado Denver
 124 *> \author NAG Ltd.
 125 *
 126 *> \date September 2012
 127 *
 128 *> \ingroup doubleOTHERauxiliary
 129 *
 130 *> \par Further Details:
 131 *  =====================
 132 *>
 133 *> \verbatim
 134 *>
 135 *>  The matrix Q is represented as a product of nb elementary reflectors
 136 *>
 137 *>     Q = H(1) H(2) . . . H(nb).
 138 *>
 139 *>  Each H(i) has the form
 140 *>
 141 *>     H(i) = I - tau * v * v**T
 142 *>
 143 *>  where tau is a real scalar, and v is a real vector with
 144 *>  v(1:i+k-1) = 0, v(i+k) = 1; v(i+k+1:n) is stored on exit in
 145 *>  A(i+k+1:n,i), and tau in TAU(i).
 146 *>
 147 *>  The elements of the vectors v together form the (n-k+1)-by-nb matrix
 148 *>  V which is needed, with T and Y, to apply the transformation to the
 149 *>  unreduced part of the matrix, using an update of the form:
 150 *>  A := (I - V*T*V**T) * (A - Y*V**T).
 151 *>
 152 *>  The contents of A on exit are illustrated by the following example
 153 *>  with n = 7, k = 3 and nb = 2:
 154 *>
 155 *>     ( a   a   a   a   a )
 156 *>     ( a   a   a   a   a )
 157 *>     ( a   a   a   a   a )
 158 *>     ( h   h   a   a   a )
 159 *>     ( v1  h   a   a   a )
 160 *>     ( v1  v2  a   a   a )
 161 *>     ( v1  v2  a   a   a )
 162 *>
 163 *>  where a denotes an element of the original matrix A, h denotes a
 164 *>  modified element of the upper Hessenberg matrix H, and vi denotes an
 165 *>  element of the vector defining H(i).
 166 *>
 167 *>  This subroutine is a slight modification of LAPACK-3.0's DLAHRD
 168 *>  incorporating improvements proposed by Quintana-Orti and Van de
 169 *>  Gejin. Note that the entries of A(1:K,2:NB) differ from those
 170 *>  returned by the original LAPACK-3.0's DLAHRD routine. (This
 171 *>  subroutine is not backward compatible with LAPACK-3.0's DLAHRD.)
 172 *> \endverbatim
 173 *
 174 *> \par References:
 175 *  ================
 176 *>
 177 *>  Gregorio Quintana-Orti and Robert van de Geijn, "Improving the
 178 *>  performance of reduction to Hessenberg form," ACM Transactions on
 179 *>  Mathematical Software, 32(2):180-194, June 2006.
 180 *>
 181 *  =====================================================================
 182       SUBROUTINE DLAHR2( N, K, NB, A, LDA, TAU, T, LDT, Y, LDY )
 183 *
 184 *  -- LAPACK auxiliary routine (version 3.4.2) --
 185 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
 186 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
 187 *     September 2012
 188 *
 189 *     .. Scalar Arguments ..
 190       INTEGER            K, LDA, LDT, LDY, N, NB
 191 *     ..
 192 *     .. Array Arguments ..
 193       DOUBLE PRECISION  A( LDA, * ), T( LDT, NB ), TAU( NB ),
 194      $                   Y( LDY, NB )
 195 *     ..
 196 *
 197 *  =====================================================================
 198 *
 199 *     .. Parameters ..
 200       DOUBLE PRECISION  ZERO, ONE
 201       PARAMETER          ( ZERO = 0.0D+0,
 202      $                     ONE = 1.0D+0 )
 203 *     ..
 204 *     .. Local Scalars ..
 205       INTEGER            I
 206       DOUBLE PRECISION  EI
 207 *     ..
 208 *     .. External Subroutines ..
 209       EXTERNAL           DAXPY, DCOPY, DGEMM, DGEMV, DLACPY,
 210      $                   DLARFG, DSCAL, DTRMM, DTRMV
 211 *     ..
 212 *     .. Intrinsic Functions ..
 213       INTRINSIC          MIN
 214 *     ..
 215 *     .. Executable Statements ..
 216 *
 217 *     Quick return if possible
 218 *
 219       IF( N.LE.1 )
 220      $   RETURN
 221 *
 222       DO 10 I = 1, NB
 223          IF( I.GT.1 ) THEN
 224 *
 225 *           Update A(K+1:N,I)
 226 *
 227 *           Update I-th column of A - Y * V**T
 228 *
 229             CALL DGEMV( 'NO TRANSPOSE', N-K, I-1, -ONE, Y(K+1,1), LDY,
 230      $                  A( K+I-1, 1 ), LDA, ONE, A( K+1, I ), 1 )
 231 *
 232 *           Apply I - V * T**T * V**T to this column (call it b) from the
 233 *           left, using the last column of T as workspace
 234 *
 235 *           Let  V = ( V1 )   and   b = ( b1 )   (first I-1 rows)
 236 *                    ( V2 )             ( b2 )
 237 *
 238 *           where V1 is unit lower triangular
 239 *
 240 *           w := V1**T * b1
 241 *
 242             CALL DCOPY( I-1, A( K+1, I ), 1, T( 1, NB ), 1 )
 243             CALL DTRMV( 'Lower', 'Transpose', 'UNIT',
 244      $                  I-1, A( K+1, 1 ),
 245      $                  LDA, T( 1, NB ), 1 )
 246 *
 247 *           w := w + V2**T * b2
 248 *
 249             CALL DGEMV( 'Transpose', N-K-I+1, I-1,
 250      $                  ONE, A( K+I, 1 ),
 251      $                  LDA, A( K+I, I ), 1, ONE, T( 1, NB ), 1 )
 252 *
 253 *           w := T**T * w
 254 *
 255             CALL DTRMV( 'Upper', 'Transpose', 'NON-UNIT',
 256      $                  I-1, T, LDT,
 257      $                  T( 1, NB ), 1 )
 258 *
 259 *           b2 := b2 - V2*w
 260 *
 261             CALL DGEMV( 'NO TRANSPOSE', N-K-I+1, I-1, -ONE,
 262      $                  A( K+I, 1 ),
 263      $                  LDA, T( 1, NB ), 1, ONE, A( K+I, I ), 1 )
 264 *
 265 *           b1 := b1 - V1*w
 266 *
 267             CALL DTRMV( 'Lower', 'NO TRANSPOSE',
 268      $                  'UNIT', I-1,
 269      $                  A( K+1, 1 ), LDA, T( 1, NB ), 1 )
 270             CALL DAXPY( I-1, -ONE, T( 1, NB ), 1, A( K+1, I ), 1 )
 271 *
 272             A( K+I-1, I-1 ) = EI
 273          END IF
 274 *
 275 *        Generate the elementary reflector H(I) to annihilate
 276 *        A(K+I+1:N,I)
 277 *
 278          CALL DLARFG( N-K-I+1, A( K+I, I ), A( MIN( K+I+1, N ), I ), 1,
 279      $                TAU( I ) )
 280          EI = A( K+I, I )
 281          A( K+I, I ) = ONE
 282 *
 283 *        Compute  Y(K+1:N,I)
 284 *
 285          CALL DGEMV( 'NO TRANSPOSE', N-K, N-K-I+1,
 286      $               ONE, A( K+1, I+1 ),
 287      $               LDA, A( K+I, I ), 1, ZERO, Y( K+1, I ), 1 )
 288          CALL DGEMV( 'Transpose', N-K-I+1, I-1,
 289      $               ONE, A( K+I, 1 ), LDA,
 290      $               A( K+I, I ), 1, ZERO, T( 1, I ), 1 )
 291          CALL DGEMV( 'NO TRANSPOSE', N-K, I-1, -ONE,
 292      $               Y( K+1, 1 ), LDY,
 293      $               T( 1, I ), 1, ONE, Y( K+1, I ), 1 )
 294          CALL DSCAL( N-K, TAU( I ), Y( K+1, I ), 1 )
 295 *
 296 *        Compute T(1:I,I)
 297 *
 298          CALL DSCAL( I-1, -TAU( I ), T( 1, I ), 1 )
 299          CALL DTRMV( 'Upper', 'No Transpose', 'NON-UNIT',
 300      $               I-1, T, LDT,
 301      $               T( 1, I ), 1 )
 302          T( I, I ) = TAU( I )
 303 *
 304    10 CONTINUE
 305       A( K+NB, NB ) = EI
 306 *
 307 *     Compute Y(1:K,1:NB)
 308 *
 309       CALL DLACPY( 'ALL', K, NB, A( 1, 2 ), LDA, Y, LDY )
 310       CALL DTRMM( 'RIGHT', 'Lower', 'NO TRANSPOSE',
 311      $            'UNIT', K, NB,
 312      $            ONE, A( K+1, 1 ), LDA, Y, LDY )
 313       IF( N.GT.K+NB )
 314      $   CALL DGEMM( 'NO TRANSPOSE', 'NO TRANSPOSE', K,
 315      $               NB, N-K-NB, ONE,
 316      $               A( 1, 2+NB ), LDA, A( K+1+NB, 1 ), LDA, ONE, Y,
 317      $               LDY )
 318       CALL DTRMM( 'RIGHT', 'Upper', 'NO TRANSPOSE',
 319      $            'NON-UNIT', K, NB,
 320      $            ONE, T, LDT, Y, LDY )
 321 *
 322       RETURN
 323 *
 324 *     End of DLAHR2
 325 *
 326       END