blob: 8a4fe2b5ff4aabef089733c86c8dd46f4de1aff9 [file] [log] [blame]
 // Copyright ©2015 The Gonum Authors. All rights reserved. // Use of this source code is governed by a BSD-style // license that can be found in the LICENSE file. package gonum import "gonum.org/v1/gonum/lapack" // Dorgbr generates one of the matrices Q or Pᵀ computed by Dgebrd // computed from the decomposition Dgebrd. See Dgebd2 for the description of // Q and Pᵀ. // // If vect == lapack.GenerateQ, then a is assumed to have been an m×k matrix and // Q is of order m. If m >= k, then Dorgbr returns the first n columns of Q // where m >= n >= k. If m < k, then Dorgbr returns Q as an m×m matrix. // // If vect == lapack.GeneratePT, then A is assumed to have been a k×n matrix, and // Pᵀ is of order n. If k < n, then Dorgbr returns the first m rows of Pᵀ, // where n >= m >= k. If k >= n, then Dorgbr returns Pᵀ as an n×n matrix. // // Dorgbr is an internal routine. It is exported for testing purposes. func (impl Implementation) Dorgbr(vect lapack.GenOrtho, m, n, k int, a []float64, lda int, tau, work []float64, lwork int) { wantq := vect == lapack.GenerateQ mn := min(m, n) switch { case vect != lapack.GenerateQ && vect != lapack.GeneratePT: panic(badGenOrtho) case m < 0: panic(mLT0) case n < 0: panic(nLT0) case k < 0: panic(kLT0) case wantq && n > m: panic(nGTM) case wantq && n < min(m, k): panic("lapack: n < min(m,k)") case !wantq && m > n: panic(mGTN) case !wantq && m < min(n, k): panic("lapack: m < min(n,k)") case lda < max(1, n) && lwork != -1: // Normally, we follow the reference and require the leading // dimension to be always valid, even in case of workspace // queries. However, if a caller provided a placeholder value // for lda (and a) when doing a workspace query that didn't // fulfill the condition here, it would cause a panic. This is // exactly what Dgesvd does. panic(badLdA) case lwork < max(1, mn) && lwork != -1: panic(badLWork) case len(work) < max(1, lwork): panic(shortWork) } // Quick return if possible. work[0] = 1 if m == 0 || n == 0 { return } if wantq { if m >= k { impl.Dorgqr(m, n, k, a, lda, tau, work, -1) } else if m > 1 { impl.Dorgqr(m-1, m-1, m-1, a[lda+1:], lda, tau, work, -1) } } else { if k < n { impl.Dorglq(m, n, k, a, lda, tau, work, -1) } else if n > 1 { impl.Dorglq(n-1, n-1, n-1, a[lda+1:], lda, tau, work, -1) } } lworkopt := int(work[0]) lworkopt = max(lworkopt, mn) if lwork == -1 { work[0] = float64(lworkopt) return } switch { case len(a) < (m-1)*lda+n: panic(shortA) case wantq && len(tau) < min(m, k): panic(shortTau) case !wantq && len(tau) < min(n, k): panic(shortTau) } if wantq { // Form Q, determined by a call to Dgebrd to reduce an m×k matrix. if m >= k { impl.Dorgqr(m, n, k, a, lda, tau, work, lwork) } else { // Shift the vectors which define the elementary reflectors one // column to the right, and set the first row and column of Q to // those of the unit matrix. for j := m - 1; j >= 1; j-- { a[j] = 0 for i := j + 1; i < m; i++ { a[i*lda+j] = a[i*lda+j-1] } } a[0] = 1 for i := 1; i < m; i++ { a[i*lda] = 0 } if m > 1 { // Form Q[1:m-1, 1:m-1] impl.Dorgqr(m-1, m-1, m-1, a[lda+1:], lda, tau, work, lwork) } } } else { // Form Pᵀ, determined by a call to Dgebrd to reduce a k×n matrix. if k < n { impl.Dorglq(m, n, k, a, lda, tau, work, lwork) } else { // Shift the vectors which define the elementary reflectors one // row downward, and set the first row and column of Pᵀ to // those of the unit matrix. a[0] = 1 for i := 1; i < n; i++ { a[i*lda] = 0 } for j := 1; j < n; j++ { for i := j - 1; i >= 1; i-- { a[i*lda+j] = a[(i-1)*lda+j] } a[j] = 0 } if n > 1 { impl.Dorglq(n-1, n-1, n-1, a[lda+1:], lda, tau, work, lwork) } } } work[0] = float64(lworkopt) }