vendor/gonum.org/v1/gonum/lapack/gonum/dsytrd.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 // Dsytrd reduces a symmetric n×n matrix A to symmetric tridiagonal form by an
  13 // orthogonal similarity transformation
  14 //  Q^T * A * Q = T
  15 // where Q is an orthonormal matrix and T is symmetric and tridiagonal.
  16 //
  17 // On entry, a contains the elements of the input matrix in the triangle specified
  18 // by uplo. On exit, the diagonal and sub/super-diagonal are overwritten by the
  19 // corresponding elements of the tridiagonal matrix T. The remaining elements in
  20 // the triangle, along with the array tau, contain the data to construct Q as
  21 // the product of elementary reflectors.
  22 //
  23 // If uplo == blas.Upper, Q is constructed with
  24 //  Q = H_{n-2} * ... * H_1 * H_0
  25 // where
  26 //  H_i = I - tau_i * v * v^T
  27 // v is constructed as v[i+1:n] = 0, v[i] = 1, v[0:i-1] is stored in A[0:i-1, i+1].
  28 // The elements of A are
  29 //  [ d   e  v1  v2  v3]
  30 //  [     d   e  v2  v3]
  31 //  [         d   e  v3]
  32 //  [             d   e]
  33 //  [                 e]
  34 //
  35 // If uplo == blas.Lower, Q is constructed with
  36 //  Q = H_0 * H_1 * ... * H_{n-2}
  37 // where
  38 //  H_i = I - tau_i * v * v^T
  39 // v is constructed as v[0:i+1] = 0, v[i+1] = 1, v[i+2:n] is stored in A[i+2:n, i].
  40 // The elements of A are
  41 //  [ d                ]
  42 //  [ e   d            ]
  43 //  [v0   e   d        ]
  44 //  [v0  v1   e   d    ]
  45 //  [v0  v1  v2   e   d]
  46 //
  47 // d must have length n, and e and tau must have length n-1. Dsytrd will panic if
  48 // these conditions are not met.
  49 //
  50 // work is temporary storage, and lwork specifies the usable memory length. At minimum,
  51 // lwork >= 1, and Dsytrd will panic otherwise. The amount of blocking is
  52 // limited by the usable length.
  53 // If lwork == -1, instead of computing Dsytrd the optimal work length is stored
  54 // into work[0].
  55 //
  56 // Dsytrd is an internal routine. It is exported for testing purposes.
  57 func (impl Implementation) Dsytrd(uplo blas.Uplo, n int, a []float64, lda int, d, e, tau, work []float64, lwork int) {
  58         checkMatrix(n, n, a, lda)
  59         if len(d) < n {
  60                 panic(badD)
  61         }
  62         if len(e) < n-1 {
  63                 panic(badE)
  64         }
  65         if len(tau) < n-1 {
  66                 panic(badTau)
  67         }
  68         if len(work) < lwork {
  69                 panic(shortWork)
  70         }
  71         if lwork != -1 && lwork < 1 {
  72                 panic(badWork)
  73         }
  74
  75         var upper bool
  76         var opts string
  77         switch uplo {
  78         case blas.Upper:
  79                 upper = true
  80                 opts = "U"
  81         case blas.Lower:
  82                 opts = "L"
  83         default:
  84                 panic(badUplo)
  85         }
  86
  87         if n == 0 {
  88                 work[0] = 1
  89                 return
  90         }
  91
  92         nb := impl.Ilaenv(1, "DSYTRD", opts, n, -1, -1, -1)
  93         lworkopt := n * nb
  94         if lwork == -1 {
  95                 work[0] = float64(lworkopt)
  96                 return
  97         }
  98
  99         nx := n
 100         bi := blas64.Implementation()
 101         var ldwork int
 102         if 1 < nb && nb < n {
 103                 // Determine when to cross over from blocked to unblocked code. The last
 104                 // block is always handled by unblocked code.
 105                 opts := "L"
 106                 if upper {
 107                         opts = "U"
 108                 }
 109                 nx = max(nb, impl.Ilaenv(3, "DSYTRD", opts, n, -1, -1, -1))
 110                 if nx < n {
 111                         // Determine if workspace is large enough for blocked code.
 112                         ldwork = nb
 113                         iws := n * ldwork
 114                         if lwork < iws {
 115                                 // Not enough workspace to use optimal nb: determine the minimum
 116                                 // value of nb and reduce nb or force use of unblocked code by
 117                                 // setting nx = n.
 118                                 nb = max(lwork/n, 1)
 119                                 nbmin := impl.Ilaenv(2, "DSYTRD", opts, n, -1, -1, -1)
 120                                 if nb < nbmin {
 121                                         nx = n
 122                                 }
 123                         }
 124                 } else {
 125                         nx = n
 126                 }
 127         } else {
 128                 nb = 1
 129         }
 130         ldwork = nb
 131
 132         if upper {
 133                 // Reduce the upper triangle of A. Columns 0:kk are handled by the
 134                 // unblocked method.
 135                 var i int
 136                 kk := n - ((n-nx+nb-1)/nb)*nb
 137                 for i = n - nb; i >= kk; i -= nb {
 138                         // Reduce columns i:i+nb to tridiagonal form and form the matrix W
 139                         // which is needed to update the unreduced part of the matrix.
 140                         impl.Dlatrd(uplo, i+nb, nb, a, lda, e, tau, work, ldwork)
 141
 142                         // Update the unreduced submatrix A[0:i-1,0:i-1], using an update
 143                         // of the form A = A - V*W^T - W*V^T.
 144                         bi.Dsyr2k(uplo, blas.NoTrans, i, nb, -1, a[i:], lda, work, ldwork, 1, a, lda)
 145
 146                         // Copy superdiagonal elements back into A, and diagonal elements into D.
 147                         for j := i; j < i+nb; j++ {
 148                                 a[(j-1)*lda+j] = e[j-1]
 149                                 d[j] = a[j*lda+j]
 150                         }
 151                 }
 152                 // Use unblocked code to reduce the last or only block
 153                 // check that i == kk.
 154                 impl.Dsytd2(uplo, kk, a, lda, d, e, tau)
 155         } else {
 156                 var i int
 157                 // Reduce the lower triangle of A.
 158                 for i = 0; i < n-nx; i += nb {
 159                         // Reduce columns 0:i+nb to tridiagonal form and form the matrix W
 160                         // which is needed to update the unreduced part of the matrix.
 161                         impl.Dlatrd(uplo, n-i, nb, a[i*lda+i:], lda, e[i:], tau[i:], work, ldwork)
 162
 163                         // Update the unreduced submatrix A[i+ib:n, i+ib:n], using an update
 164                         // of the form A = A + V*W^T - W*V^T.
 165                         bi.Dsyr2k(uplo, blas.NoTrans, n-i-nb, nb, -1, a[(i+nb)*lda+i:], lda,
 166                                 work[nb*ldwork:], ldwork, 1, a[(i+nb)*lda+i+nb:], lda)
 167
 168                         // Copy subdiagonal elements back into A, and diagonal elements into D.
 169                         for j := i; j < i+nb; j++ {
 170                                 a[(j+1)*lda+j] = e[j]
 171                                 d[j] = a[j*lda+j]
 172                         }
 173                 }
 174                 // Use unblocked code to reduce the last or only block.
 175                 impl.Dsytd2(uplo, n-i, a[i*lda+i:], lda, d[i:], e[i:], tau[i:])
 176         }
 177         work[0] = float64(lworkopt)
 178 }