vendor/gonum.org/v1/gonum/lapack/gonum/dlahr2.go

   1 // Copyright ©2016 The Gonum Authors. All rights reserved.
   2 // Use of this source code is governed by a BSD-style
   3 // license that can be found in the LICENSE file.
   4
   5 package gonum
   6
   7 import (
   8         "gonum.org/v1/gonum/blas"
   9         "gonum.org/v1/gonum/blas/blas64"
  10 )
  11
  12 // Dlahr2 reduces the first nb columns of a real general n×(n-k+1) matrix A so
  13 // that elements below the k-th subdiagonal are zero. The reduction is performed
  14 // by an orthogonal similarity transformation Q^T * A * Q. Dlahr2 returns the
  15 // matrices V and T which determine Q as a block reflector I - V*T*V^T, and
  16 // also the matrix Y = A * V * T.
  17 //
  18 // The matrix Q is represented as a product of nb elementary reflectors
  19 //  Q = H_0 * H_1 * ... * H_{nb-1}.
  20 // Each H_i has the form
  21 //  H_i = I - tau[i] * v * v^T,
  22 // where v is a real vector with v[0:i+k-1] = 0 and v[i+k-1] = 1. v[i+k:n] is
  23 // stored on exit in A[i+k+1:n,i].
  24 //
  25 // The elements of the vectors v together form the (n-k+1)×nb matrix
  26 // V which is needed, with T and Y, to apply the transformation to the
  27 // unreduced part of the matrix, using an update of the form
  28 //  A = (I - V*T*V^T) * (A - Y*V^T).
  29 //
  30 // On entry, a contains the n×(n-k+1) general matrix A. On return, the elements
  31 // on and above the k-th subdiagonal in the first nb columns are overwritten
  32 // with the corresponding elements of the reduced matrix; the elements below the
  33 // k-th subdiagonal, with the slice tau, represent the matrix Q as a product of
  34 // elementary reflectors. The other columns of A are unchanged.
  35 //
  36 // The contents of A on exit are illustrated by the following example
  37 // with n = 7, k = 3 and nb = 2:
  38 //  [ a   a   a   a   a ]
  39 //  [ a   a   a   a   a ]
  40 //  [ a   a   a   a   a ]
  41 //  [ h   h   a   a   a ]
  42 //  [ v0  h   a   a   a ]
  43 //  [ v0  v1  a   a   a ]
  44 //  [ v0  v1  a   a   a ]
  45 // where a denotes an element of the original matrix A, h denotes a
  46 // modified element of the upper Hessenberg matrix H, and vi denotes an
  47 // element of the vector defining H_i.
  48 //
  49 // k is the offset for the reduction. Elements below the k-th subdiagonal in the
  50 // first nb columns are reduced to zero.
  51 //
  52 // nb is the number of columns to be reduced.
  53 //
  54 // On entry, a represents the n×(n-k+1) matrix A. On return, the elements on and
  55 // above the k-th subdiagonal in the first nb columns are overwritten with the
  56 // corresponding elements of the reduced matrix. The elements below the k-th
  57 // subdiagonal, with the slice tau, represent the matrix Q as a product of
  58 // elementary reflectors. The other columns of A are unchanged.
  59 //
  60 // tau will contain the scalar factors of the elementary reflectors. It must
  61 // have length at least nb.
  62 //
  63 // t and ldt represent the nb×nb upper triangular matrix T, and y and ldy
  64 // represent the n×nb matrix Y.
  65 //
  66 // Dlahr2 is an internal routine. It is exported for testing purposes.
  67 func (impl Implementation) Dlahr2(n, k, nb int, a []float64, lda int, tau, t []float64, ldt int, y []float64, ldy int) {
  68         checkMatrix(n, n-k+1, a, lda)
  69         if len(tau) < nb {
  70                 panic(badTau)
  71         }
  72         checkMatrix(nb, nb, t, ldt)
  73         checkMatrix(n, nb, y, ldy)
  74
  75         // Quick return if possible.
  76         if n <= 1 {
  77                 return
  78         }
  79
  80         bi := blas64.Implementation()
  81         var ei float64
  82         for i := 0; i < nb; i++ {
  83                 if i > 0 {
  84                         // Update A[k:n,i].
  85
  86                         // Update i-th column of A - Y * V^T.
  87                         bi.Dgemv(blas.NoTrans, n-k, i,
  88                                 -1, y[k*ldy:], ldy,
  89                                 a[(k+i-1)*lda:], 1,
  90                                 1, a[k*lda+i:], lda)
  91
  92                         // Apply I - V * T^T * V^T to this column (call it b)
  93                         // from the left, using the last column of T as
  94                         // workspace.
  95                         // Let V = [ V1 ]   and   b = [ b1 ]   (first i rows)
  96                         //         [ V2 ]             [ b2 ]
  97                         // where V1 is unit lower triangular.
  98                         //
  99                         // w := V1^T * b1.
 100                         bi.Dcopy(i, a[k*lda+i:], lda, t[nb-1:], ldt)
 101                         bi.Dtrmv(blas.Lower, blas.Trans, blas.Unit, i,
 102                                 a[k*lda:], lda, t[nb-1:], ldt)
 103
 104                         // w := w + V2^T * b2.
 105                         bi.Dgemv(blas.Trans, n-k-i, i,
 106                                 1, a[(k+i)*lda:], lda,
 107                                 a[(k+i)*lda+i:], lda,
 108                                 1, t[nb-1:], ldt)
 109
 110                         // w := T^T * w.
 111                         bi.Dtrmv(blas.Upper, blas.Trans, blas.NonUnit, i,
 112                                 t, ldt, t[nb-1:], ldt)
 113
 114                         // b2 := b2 - V2*w.
 115                         bi.Dgemv(blas.NoTrans, n-k-i, i,
 116                                 -1, a[(k+i)*lda:], lda,
 117                                 t[nb-1:], ldt,
 118                                 1, a[(k+i)*lda+i:], lda)
 119
 120                         // b1 := b1 - V1*w.
 121                         bi.Dtrmv(blas.Lower, blas.NoTrans, blas.Unit, i,
 122                                 a[k*lda:], lda, t[nb-1:], ldt)
 123                         bi.Daxpy(i, -1, t[nb-1:], ldt, a[k*lda+i:], lda)
 124
 125                         a[(k+i-1)*lda+i-1] = ei
 126                 }
 127
 128                 // Generate the elementary reflector H_i to annihilate
 129                 // A[k+i+1:n,i].
 130                 ei, tau[i] = impl.Dlarfg(n-k-i, a[(k+i)*lda+i], a[min(k+i+1, n-1)*lda+i:], lda)
 131                 a[(k+i)*lda+i] = 1
 132
 133                 // Compute Y[k:n,i].
 134                 bi.Dgemv(blas.NoTrans, n-k, n-k-i,
 135                         1, a[k*lda+i+1:], lda,
 136                         a[(k+i)*lda+i:], lda,
 137                         0, y[k*ldy+i:], ldy)
 138                 bi.Dgemv(blas.Trans, n-k-i, i,
 139                         1, a[(k+i)*lda:], lda,
 140                         a[(k+i)*lda+i:], lda,
 141                         0, t[i:], ldt)
 142                 bi.Dgemv(blas.NoTrans, n-k, i,
 143                         -1, y[k*ldy:], ldy,
 144                         t[i:], ldt,
 145                         1, y[k*ldy+i:], ldy)
 146                 bi.Dscal(n-k, tau[i], y[k*ldy+i:], ldy)
 147
 148                 // Compute T[0:i,i].
 149                 bi.Dscal(i, -tau[i], t[i:], ldt)
 150                 bi.Dtrmv(blas.Upper, blas.NoTrans, blas.NonUnit, i,
 151                         t, ldt, t[i:], ldt)
 152
 153                 t[i*ldt+i] = tau[i]
 154         }
 155         a[(k+nb-1)*lda+nb-1] = ei
 156
 157         // Compute Y[0:k,0:nb].
 158         impl.Dlacpy(blas.All, k, nb, a[1:], lda, y, ldy)
 159         bi.Dtrmm(blas.Right, blas.Lower, blas.NoTrans, blas.Unit, k, nb,
 160                 1, a[k*lda:], lda, y, ldy)
 161         if n > k+nb {
 162                 bi.Dgemm(blas.NoTrans, blas.NoTrans, k, nb, n-k-nb,
 163                         1, a[1+nb:], lda,
 164                         a[(k+nb)*lda:], lda,
 165                         1, y, ldy)
 166         }
 167         bi.Dtrmm(blas.Right, blas.Upper, blas.NoTrans, blas.NonUnit, k, nb,
 168                 1, t, ldt, y, ldy)
 169 }