vendor/gonum.org/v1/gonum/lapack/gonum/dlatrd.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 // Dlatrd reduces nb rows and columns of a real n×n symmetric matrix A to symmetric
  13 // tridiagonal form. It computes the orthonormal similarity transformation
  14 //  Q^T * A * Q
  15 // and returns the matrices V and W to apply to the unreduced part of A. If
  16 // uplo == blas.Upper, the upper triangle is supplied and the last nb rows are
  17 // reduced. If uplo == blas.Lower, the lower triangle is supplied and the first
  18 // nb rows are reduced.
  19 //
  20 // a contains the symmetric matrix on entry with active triangular half specified
  21 // by uplo. On exit, the nb columns have been reduced to tridiagonal form. The
  22 // diagonal contains the diagonal of the reduced matrix, the off-diagonal is
  23 // set to 1, and the remaining elements contain the data to construct Q.
  24 //
  25 // If uplo == blas.Upper, with n = 5 and nb = 2 on exit a is
  26 //  [ a   a   a  v4  v5]
  27 //  [     a   a  v4  v5]
  28 //  [         a   1  v5]
  29 //  [             d   1]
  30 //  [                 d]
  31 //
  32 // If uplo == blas.Lower, with n = 5 and nb = 2, on exit a is
  33 //  [ d                ]
  34 //  [ 1   d            ]
  35 //  [v1   1   a        ]
  36 //  [v1  v2   a   a    ]
  37 //  [v1  v2   a   a   a]
  38 //
  39 // e contains the superdiagonal elements of the reduced matrix. If uplo == blas.Upper,
  40 // e[n-nb:n-1] contains the last nb columns of the reduced matrix, while if
  41 // uplo == blas.Lower, e[:nb] contains the first nb columns of the reduced matrix.
  42 // e must have length at least n-1, and Dlatrd will panic otherwise.
  43 //
  44 // tau contains the scalar factors of the elementary reflectors needed to construct Q.
  45 // The reflectors are stored in tau[n-nb:n-1] if uplo == blas.Upper, and in
  46 // tau[:nb] if uplo == blas.Lower. tau must have length n-1, and Dlatrd will panic
  47 // otherwise.
  48 //
  49 // w is an n×nb matrix. On exit it contains the data to update the unreduced part
  50 // of A.
  51 //
  52 // The matrix Q is represented as a product of elementary reflectors. Each reflector
  53 // H has the form
  54 //  I - tau * v * v^T
  55 // If uplo == blas.Upper,
  56 //  Q = H_{n-1} * H_{n-2} * ... * H_{n-nb}
  57 // where v[:i-1] is stored in A[:i-1,i], v[i-1] = 1, and v[i:n] = 0.
  58 //
  59 // If uplo == blas.Lower,
  60 //  Q = H_0 * H_1 * ... * H_{nb-1}
  61 // where v[:i+1] = 0, v[i+1] = 1, and v[i+2:n] is stored in A[i+2:n,i].
  62 //
  63 // The vectors v form the n×nb matrix V which is used with W to apply a
  64 // symmetric rank-2 update to the unreduced part of A
  65 //  A = A - V * W^T - W * V^T
  66 //
  67 // Dlatrd is an internal routine. It is exported for testing purposes.
  68 func (impl Implementation) Dlatrd(uplo blas.Uplo, n, nb int, a []float64, lda int, e, tau, w []float64, ldw int) {
  69         checkMatrix(n, n, a, lda)
  70         checkMatrix(n, nb, w, ldw)
  71         if len(e) < n-1 {
  72                 panic(badE)
  73         }
  74         if len(tau) < n-1 {
  75                 panic(badTau)
  76         }
  77         if n <= 0 {
  78                 return
  79         }
  80         bi := blas64.Implementation()
  81         if uplo == blas.Upper {
  82                 for i := n - 1; i >= n-nb; i-- {
  83                         iw := i - n + nb
  84                         if i < n-1 {
  85                                 // Update A(0:i, i).
  86                                 bi.Dgemv(blas.NoTrans, i+1, n-i-1, -1, a[i+1:], lda,
  87                                         w[i*ldw+iw+1:], 1, 1, a[i:], lda)
  88                                 bi.Dgemv(blas.NoTrans, i+1, n-i-1, -1, w[iw+1:], ldw,
  89                                         a[i*lda+i+1:], 1, 1, a[i:], lda)
  90                         }
  91                         if i > 0 {
  92                                 // Generate elementary reflector H_i to annihilate A(0:i-2,i).
  93                                 e[i-1], tau[i-1] = impl.Dlarfg(i, a[(i-1)*lda+i], a[i:], lda)
  94                                 a[(i-1)*lda+i] = 1
  95
  96                                 // Compute W(0:i-1, i).
  97                                 bi.Dsymv(blas.Upper, i, 1, a, lda, a[i:], lda, 0, w[iw:], ldw)
  98                                 if i < n-1 {
  99                                         bi.Dgemv(blas.Trans, i, n-i-1, 1, w[iw+1:], ldw,
 100                                                 a[i:], lda, 0, w[(i+1)*ldw+iw:], ldw)
 101                                         bi.Dgemv(blas.NoTrans, i, n-i-1, -1, a[i+1:], lda,
 102                                                 w[(i+1)*ldw+iw:], ldw, 1, w[iw:], ldw)
 103                                         bi.Dgemv(blas.Trans, i, n-i-1, 1, a[i+1:], lda,
 104                                                 a[i:], lda, 0, w[(i+1)*ldw+iw:], ldw)
 105                                         bi.Dgemv(blas.NoTrans, i, n-i-1, -1, w[iw+1:], ldw,
 106                                                 w[(i+1)*ldw+iw:], ldw, 1, w[iw:], ldw)
 107                                 }
 108                                 bi.Dscal(i, tau[i-1], w[iw:], ldw)
 109                                 alpha := -0.5 * tau[i-1] * bi.Ddot(i, w[iw:], ldw, a[i:], lda)
 110                                 bi.Daxpy(i, alpha, a[i:], lda, w[iw:], ldw)
 111                         }
 112                 }
 113         } else {
 114                 // Reduce first nb columns of lower triangle.
 115                 for i := 0; i < nb; i++ {
 116                         // Update A(i:n, i)
 117                         bi.Dgemv(blas.NoTrans, n-i, i, -1, a[i*lda:], lda,
 118                                 w[i*ldw:], 1, 1, a[i*lda+i:], lda)
 119                         bi.Dgemv(blas.NoTrans, n-i, i, -1, w[i*ldw:], ldw,
 120                                 a[i*lda:], 1, 1, a[i*lda+i:], lda)
 121                         if i < n-1 {
 122                                 // Generate elementary reflector H_i to annihilate A(i+2:n,i).
 123                                 e[i], tau[i] = impl.Dlarfg(n-i-1, a[(i+1)*lda+i], a[min(i+2, n-1)*lda+i:], lda)
 124                                 a[(i+1)*lda+i] = 1
 125
 126                                 // Compute W(i+1:n,i).
 127                                 bi.Dsymv(blas.Lower, n-i-1, 1, a[(i+1)*lda+i+1:], lda,
 128                                         a[(i+1)*lda+i:], lda, 0, w[(i+1)*ldw+i:], ldw)
 129                                 bi.Dgemv(blas.Trans, n-i-1, i, 1, w[(i+1)*ldw:], ldw,
 130                                         a[(i+1)*lda+i:], lda, 0, w[i:], ldw)
 131                                 bi.Dgemv(blas.NoTrans, n-i-1, i, -1, a[(i+1)*lda:], lda,
 132                                         w[i:], ldw, 1, w[(i+1)*ldw+i:], ldw)
 133                                 bi.Dgemv(blas.Trans, n-i-1, i, 1, a[(i+1)*lda:], lda,
 134                                         a[(i+1)*lda+i:], lda, 0, w[i:], ldw)
 135                                 bi.Dgemv(blas.NoTrans, n-i-1, i, -1, w[(i+1)*ldw:], ldw,
 136                                         w[i:], ldw, 1, w[(i+1)*ldw+i:], ldw)
 137                                 bi.Dscal(n-i-1, tau[i], w[(i+1)*ldw+i:], ldw)
 138                                 alpha := -0.5 * tau[i] * bi.Ddot(n-i-1, w[(i+1)*ldw+i:], ldw,
 139                                         a[(i+1)*lda+i:], lda)
 140                                 bi.Daxpy(n-i-1, alpha, a[(i+1)*lda+i:], lda,
 141                                         w[(i+1)*ldw+i:], ldw)
 142                         }
 143                 }
 144         }
 145 }