vendor/gonum.org/v1/gonum/lapack/gonum/dsytd2.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 // Dsytd2 reduces a symmetric n×n matrix A to symmetric tridiagonal form T by an
  13 // orthogonal similarity transformation
  14 //  Q^T * A * Q = T
  15 // On entry, the matrix is contained in the specified triangle of a. On exit,
  16 // if uplo == blas.Upper, the diagonal and first super-diagonal of a are
  17 // overwritten with the elements of T. The elements above the first super-diagonal
  18 // are overwritten with the the elementary reflectors that are used with the
  19 // elements written to tau in order to construct Q. If uplo == blas.Lower, the
  20 // elements are written in the lower triangular region.
  21 //
  22 // d must have length at least n. e and tau must have length at least n-1. Dsytd2
  23 // will panic if these sizes are not met.
  24 //
  25 // Q is represented as a product of elementary reflectors.
  26 // If uplo == blas.Upper
  27 //  Q = H_{n-2} * ... * H_1 * H_0
  28 // and if uplo == blas.Lower
  29 //  Q = H_0 * H_1 * ... * H_{n-2}
  30 // where
  31 //  H_i = I - tau * v * v^T
  32 // where tau is stored in tau[i], and v is stored in a.
  33 //
  34 // If uplo == blas.Upper, v[0:i-1] is stored in A[0:i-1,i+1], v[i] = 1, and
  35 // v[i+1:] = 0. The elements of a are
  36 //  [ d   e  v2  v3  v4]
  37 //  [     d   e  v3  v4]
  38 //  [         d   e  v4]
  39 //  [             d   e]
  40 //  [                 d]
  41 // If uplo == blas.Lower, v[0:i+1] = 0, v[i+1] = 1, and v[i+2:] is stored in
  42 // A[i+2:n,i].
  43 // The elements of a are
  44 //  [ d                ]
  45 //  [ e   d            ]
  46 //  [v1   e   d        ]
  47 //  [v1  v2   e   d    ]
  48 //  [v1  v2  v3   e   d]
  49 //
  50 // Dsytd2 is an internal routine. It is exported for testing purposes.
  51 func (impl Implementation) Dsytd2(uplo blas.Uplo, n int, a []float64, lda int, d, e, tau []float64) {
  52         checkMatrix(n, n, a, lda)
  53         if len(d) < n {
  54                 panic(badD)
  55         }
  56         if len(e) < n-1 {
  57                 panic(badE)
  58         }
  59         if len(tau) < n-1 {
  60                 panic(badTau)
  61         }
  62         if n <= 0 {
  63                 return
  64         }
  65         bi := blas64.Implementation()
  66         if uplo == blas.Upper {
  67                 // Reduce the upper triangle of A.
  68                 for i := n - 2; i >= 0; i-- {
  69                         // Generate elementary reflector H_i = I - tau * v * v^T to
  70                         // annihilate A[i:i-1, i+1].
  71                         var taui float64
  72                         a[i*lda+i+1], taui = impl.Dlarfg(i+1, a[i*lda+i+1], a[i+1:], lda)
  73                         e[i] = a[i*lda+i+1]
  74                         if taui != 0 {
  75                                 // Apply H_i from both sides to A[0:i,0:i].
  76                                 a[i*lda+i+1] = 1
  77
  78                                 // Compute x := tau * A * v storing x in tau[0:i].
  79                                 bi.Dsymv(uplo, i+1, taui, a, lda, a[i+1:], lda, 0, tau, 1)
  80
  81                                 // Compute w := x - 1/2 * tau * (x^T * v) * v.
  82                                 alpha := -0.5 * taui * bi.Ddot(i+1, tau, 1, a[i+1:], lda)
  83                                 bi.Daxpy(i+1, alpha, a[i+1:], lda, tau, 1)
  84
  85                                 // Apply the transformation as a rank-2 update
  86                                 // A = A - v * w^T - w * v^T.
  87                                 bi.Dsyr2(uplo, i+1, -1, a[i+1:], lda, tau, 1, a, lda)
  88                                 a[i*lda+i+1] = e[i]
  89                         }
  90                         d[i+1] = a[(i+1)*lda+i+1]
  91                         tau[i] = taui
  92                 }
  93                 d[0] = a[0]
  94                 return
  95         }
  96         // Reduce the lower triangle of A.
  97         for i := 0; i < n-1; i++ {
  98                 // Generate elementary reflector H_i = I - tau * v * v^T to
  99                 // annihilate A[i+2:n, i].
 100                 var taui float64
 101                 a[(i+1)*lda+i], taui = impl.Dlarfg(n-i-1, a[(i+1)*lda+i], a[min(i+2, n-1)*lda+i:], lda)
 102                 e[i] = a[(i+1)*lda+i]
 103                 if taui != 0 {
 104                         // Apply H_i from both sides to A[i+1:n, i+1:n].
 105                         a[(i+1)*lda+i] = 1
 106
 107                         // Compute x := tau * A * v, storing y in tau[i:n-1].
 108                         bi.Dsymv(uplo, n-i-1, taui, a[(i+1)*lda+i+1:], lda, a[(i+1)*lda+i:], lda, 0, tau[i:], 1)
 109
 110                         // Compute w := x - 1/2 * tau * (x^T * v) * v.
 111                         alpha := -0.5 * taui * bi.Ddot(n-i-1, tau[i:], 1, a[(i+1)*lda+i:], lda)
 112                         bi.Daxpy(n-i-1, alpha, a[(i+1)*lda+i:], lda, tau[i:], 1)
 113
 114                         // Apply the transformation as a rank-2 update
 115                         // A = A - v * w^T - w * v^T.
 116                         bi.Dsyr2(uplo, n-i-1, -1, a[(i+1)*lda+i:], lda, tau[i:], 1, a[(i+1)*lda+i+1:], lda)
 117                         a[(i+1)*lda+i] = e[i]
 118                 }
 119                 d[i] = a[i*lda+i]
 120                 tau[i] = taui
 121         }
 122         d[n-1] = a[(n-1)*lda+n-1]
 123 }