forked from LaconicNetwork/kompose
Upgrade OpenShift and its dependencies.
OpenShift version 1.4.0-alpha.0
This commit is contained in:
Tomas Kral
+44
@@ -0,0 +1,44 @@
|
||||
// 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 native
|
||||
|
||||
import "github.com/gonum/blas"
|
||||
|
||||
// Dgelq2 computes the LQ factorization of the m×n matrix a.
|
||||
//
|
||||
// During Dgelq2, a is modified to contain the information to construct Q and L.
|
||||
// The lower triangle of a contains the matrix L. The upper triangular elements
|
||||
// (not including the diagonal) contain the elementary reflectors. Tau is modified
|
||||
// to contain the reflector scales. Tau must have length of at least k = min(m,n)
|
||||
// and this function will panic otherwise.
|
||||
//
|
||||
// See Dgeqr2 for a description of the elementary reflectors and orthonormal
|
||||
// matrix Q. Q is constructed as a product of these elementary reflectors,
|
||||
// Q = H_k ... H_2*H_1.
|
||||
//
|
||||
// Work is temporary storage of length at least m and this function will panic otherwise.
|
||||
func (impl Implementation) Dgelq2(m, n int, a []float64, lda int, tau, work []float64) {
|
||||
checkMatrix(m, n, a, lda)
|
||||
k := min(m, n)
|
||||
if len(tau) < k {
|
||||
panic(badTau)
|
||||
}
|
||||
if len(work) < m {
|
||||
panic(badWork)
|
||||
}
|
||||
for i := 0; i < k; i++ {
|
||||
a[i*lda+i], tau[i] = impl.Dlarfg(n-i, a[i*lda+i], a[i*lda+min(i+1, n-1):], 1)
|
||||
if i < m-1 {
|
||||
aii := a[i*lda+i]
|
||||
a[i*lda+i] = 1
|
||||
impl.Dlarf(blas.Right, m-i-1, n-i,
|
||||
a[i*lda+i:], 1,
|
||||
tau[i],
|
||||
a[(i+1)*lda+i:], lda,
|
||||
work)
|
||||
a[i*lda+i] = aii
|
||||
}
|
||||
}
|
||||
}
|
||||
+84
@@ -0,0 +1,84 @@
|
||||
// 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 native
|
||||
|
||||
import (
|
||||
"github.com/gonum/blas"
|
||||
"github.com/gonum/lapack"
|
||||
)
|
||||
|
||||
// Dgelqf computes the LQ factorization of the m×n matrix a using a blocked
|
||||
// algorithm. Please see the documentation for Dgelq2 for a description of the
|
||||
// parameters at entry and exit.
|
||||
//
|
||||
// Work is temporary storage, and lwork specifies the usable memory length.
|
||||
// At minimum, lwork >= m, and this function will panic otherwise.
|
||||
// Dgelqf is a blocked LQ factorization, but the block size is limited
|
||||
// by the temporary space available. If lwork == -1, instead of performing Dgelqf,
|
||||
// the optimal work length will be stored into work[0].
|
||||
//
|
||||
// tau must have length at least min(m,n), and this function will panic otherwise.
|
||||
func (impl Implementation) Dgelqf(m, n int, a []float64, lda int, tau, work []float64, lwork int) {
|
||||
nb := impl.Ilaenv(1, "DGELQF", " ", m, n, -1, -1)
|
||||
lworkopt := m * max(nb, 1)
|
||||
if lwork == -1 {
|
||||
work[0] = float64(lworkopt)
|
||||
return
|
||||
}
|
||||
checkMatrix(m, n, a, lda)
|
||||
if len(work) < lwork {
|
||||
panic(shortWork)
|
||||
}
|
||||
if lwork < m {
|
||||
panic(badWork)
|
||||
}
|
||||
k := min(m, n)
|
||||
if len(tau) < k {
|
||||
panic(badTau)
|
||||
}
|
||||
if k == 0 {
|
||||
return
|
||||
}
|
||||
// Find the optimal blocking size based on the size of available memory
|
||||
// and optimal machine parameters.
|
||||
nbmin := 2
|
||||
var nx int
|
||||
iws := m
|
||||
ldwork := nb
|
||||
if nb > 1 && k > nb {
|
||||
nx = max(0, impl.Ilaenv(3, "DGELQF", " ", m, n, -1, -1))
|
||||
if nx < k {
|
||||
iws = m * nb
|
||||
if lwork < iws {
|
||||
nb = lwork / m
|
||||
nbmin = max(2, impl.Ilaenv(2, "DGELQF", " ", m, n, -1, -1))
|
||||
}
|
||||
}
|
||||
}
|
||||
// Computed blocked LQ factorization.
|
||||
var i int
|
||||
if nb >= nbmin && nb < k && nx < k {
|
||||
for i = 0; i < k-nx; i += nb {
|
||||
ib := min(k-i, nb)
|
||||
impl.Dgelq2(ib, n-i, a[i*lda+i:], lda, tau[i:], work)
|
||||
if i+ib < m {
|
||||
impl.Dlarft(lapack.Forward, lapack.RowWise, n-i, ib,
|
||||
a[i*lda+i:], lda,
|
||||
tau[i:],
|
||||
work, ldwork)
|
||||
impl.Dlarfb(blas.Right, blas.NoTrans, lapack.Forward, lapack.RowWise,
|
||||
m-i-ib, n-i, ib,
|
||||
a[i*lda+i:], lda,
|
||||
work, ldwork,
|
||||
a[(i+ib)*lda+i:], lda,
|
||||
work[ib*ldwork:], ldwork)
|
||||
}
|
||||
}
|
||||
}
|
||||
// Perform unblocked LQ factorization on the remainder.
|
||||
if i < k {
|
||||
impl.Dgelq2(m-i, n-i, a[i*lda+i:], lda, tau[i:], work)
|
||||
}
|
||||
}
|
||||
+200
@@ -0,0 +1,200 @@
|
||||
// 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 native
|
||||
|
||||
import (
|
||||
"github.com/gonum/blas"
|
||||
"github.com/gonum/lapack"
|
||||
)
|
||||
|
||||
// Dgels finds a minimum-norm solution based on the matrices a and b using the
|
||||
// QR or LQ factorization. Dgels returns false if the matrix
|
||||
// A is singular, and true if this solution was successfully found.
|
||||
//
|
||||
// The minimization problem solved depends on the input parameters.
|
||||
//
|
||||
// 1. If m >= n and trans == blas.NoTrans, Dgels finds X such that || A*X - B||_2
|
||||
// is minimized.
|
||||
// 2. If m < n and trans == blas.NoTrans, Dgels finds the minimum norm solution of
|
||||
// A * X = B.
|
||||
// 3. If m >= n and trans == blas.Trans, Dgels finds the minimum norm solution of
|
||||
// A^T * X = B.
|
||||
// 4. If m < n and trans == blas.Trans, Dgels finds X such that || A*X - B||_2
|
||||
// is minimized.
|
||||
// Note that the least-squares solutions (cases 1 and 3) perform the minimization
|
||||
// per column of B. This is not the same as finding the minimum-norm matrix.
|
||||
//
|
||||
// The matrix a is a general matrix of size m×n and is modified during this call.
|
||||
// The input matrix b is of size max(m,n)×nrhs, and serves two purposes. On entry,
|
||||
// the elements of b specify the input matrix B. B has size m×nrhs if
|
||||
// trans == blas.NoTrans, and n×nrhs if trans == blas.Trans. On exit, the
|
||||
// leading submatrix of b contains the solution vectors X. If trans == blas.NoTrans,
|
||||
// this submatrix is of size n×nrhs, and of size m×nrhs otherwise.
|
||||
//
|
||||
// Work is temporary storage, and lwork specifies the usable memory length.
|
||||
// At minimum, lwork >= max(m,n) + max(m,n,nrhs), and this function will panic
|
||||
// otherwise. A longer work will enable blocked algorithms to be called.
|
||||
// In the special case that lwork == -1, work[0] will be set to the optimal working
|
||||
// length.
|
||||
func (impl Implementation) Dgels(trans blas.Transpose, m, n, nrhs int, a []float64, lda int, b []float64, ldb int, work []float64, lwork int) bool {
|
||||
notran := trans == blas.NoTrans
|
||||
checkMatrix(m, n, a, lda)
|
||||
mn := min(m, n)
|
||||
checkMatrix(mn, nrhs, b, ldb)
|
||||
|
||||
// Find optimal block size.
|
||||
tpsd := true
|
||||
if notran {
|
||||
tpsd = false
|
||||
}
|
||||
var nb int
|
||||
if m >= n {
|
||||
nb = impl.Ilaenv(1, "DGEQRF", " ", m, n, -1, -1)
|
||||
if tpsd {
|
||||
nb = max(nb, impl.Ilaenv(1, "DORMQR", "LN", m, nrhs, n, -1))
|
||||
} else {
|
||||
nb = max(nb, impl.Ilaenv(1, "DORMQR", "LT", m, nrhs, n, -1))
|
||||
}
|
||||
} else {
|
||||
nb = impl.Ilaenv(1, "DGELQF", " ", m, n, -1, -1)
|
||||
if tpsd {
|
||||
nb = max(nb, impl.Ilaenv(1, "DORMLQ", "LT", n, nrhs, m, -1))
|
||||
} else {
|
||||
nb = max(nb, impl.Ilaenv(1, "DORMLQ", "LN", n, nrhs, m, -1))
|
||||
}
|
||||
}
|
||||
if lwork == -1 {
|
||||
work[0] = float64(max(1, mn+max(mn, nrhs)*nb))
|
||||
return true
|
||||
}
|
||||
|
||||
if len(work) < lwork {
|
||||
panic(shortWork)
|
||||
}
|
||||
if lwork < mn+max(mn, nrhs) {
|
||||
panic(badWork)
|
||||
}
|
||||
if m == 0 || n == 0 || nrhs == 0 {
|
||||
impl.Dlaset(blas.All, max(m, n), nrhs, 0, 0, b, ldb)
|
||||
return true
|
||||
}
|
||||
|
||||
// Scale the input matrices if they contain extreme values.
|
||||
smlnum := dlamchS / dlamchP
|
||||
bignum := 1 / smlnum
|
||||
anrm := impl.Dlange(lapack.MaxAbs, m, n, a, lda, nil)
|
||||
var iascl int
|
||||
if anrm > 0 && anrm < smlnum {
|
||||
impl.Dlascl(lapack.General, 0, 0, anrm, smlnum, m, n, a, lda)
|
||||
iascl = 1
|
||||
} else if anrm > bignum {
|
||||
impl.Dlascl(lapack.General, 0, 0, anrm, bignum, m, n, a, lda)
|
||||
} else if anrm == 0 {
|
||||
// Matrix all zeros
|
||||
impl.Dlaset(blas.All, max(m, n), nrhs, 0, 0, b, ldb)
|
||||
return true
|
||||
}
|
||||
brow := m
|
||||
if tpsd {
|
||||
brow = n
|
||||
}
|
||||
bnrm := impl.Dlange(lapack.MaxAbs, brow, nrhs, b, ldb, nil)
|
||||
ibscl := 0
|
||||
if bnrm > 0 && bnrm < smlnum {
|
||||
impl.Dlascl(lapack.General, 0, 0, bnrm, smlnum, brow, nrhs, b, ldb)
|
||||
ibscl = 1
|
||||
} else if bnrm > bignum {
|
||||
impl.Dlascl(lapack.General, 0, 0, bnrm, bignum, brow, nrhs, b, ldb)
|
||||
ibscl = 2
|
||||
}
|
||||
|
||||
// Solve the minimization problem using a QR or an LQ decomposition.
|
||||
var scllen int
|
||||
if m >= n {
|
||||
impl.Dgeqrf(m, n, a, lda, work, work[mn:], lwork-mn)
|
||||
if !tpsd {
|
||||
impl.Dormqr(blas.Left, blas.Trans, m, nrhs, n,
|
||||
a, lda,
|
||||
work,
|
||||
b, ldb,
|
||||
work[mn:], lwork-mn)
|
||||
ok := impl.Dtrtrs(blas.Upper, blas.NoTrans, blas.NonUnit, n, nrhs,
|
||||
a, lda,
|
||||
b, ldb)
|
||||
if !ok {
|
||||
return false
|
||||
}
|
||||
scllen = n
|
||||
} else {
|
||||
ok := impl.Dtrtrs(blas.Upper, blas.Trans, blas.NonUnit, n, nrhs,
|
||||
a, lda,
|
||||
b, ldb)
|
||||
if !ok {
|
||||
return false
|
||||
}
|
||||
for i := n; i < m; i++ {
|
||||
for j := 0; j < nrhs; j++ {
|
||||
b[i*ldb+j] = 0
|
||||
}
|
||||
}
|
||||
impl.Dormqr(blas.Left, blas.NoTrans, m, nrhs, n,
|
||||
a, lda,
|
||||
work,
|
||||
b, ldb,
|
||||
work[mn:], lwork-mn)
|
||||
scllen = m
|
||||
}
|
||||
} else {
|
||||
impl.Dgelqf(m, n, a, lda, work, work[mn:], lwork-mn)
|
||||
if !tpsd {
|
||||
ok := impl.Dtrtrs(blas.Lower, blas.NoTrans, blas.NonUnit,
|
||||
m, nrhs,
|
||||
a, lda,
|
||||
b, ldb)
|
||||
if !ok {
|
||||
return false
|
||||
}
|
||||
for i := m; i < n; i++ {
|
||||
for j := 0; j < nrhs; j++ {
|
||||
b[i*ldb+j] = 0
|
||||
}
|
||||
}
|
||||
impl.Dormlq(blas.Left, blas.Trans, n, nrhs, m,
|
||||
a, lda,
|
||||
work,
|
||||
b, ldb,
|
||||
work[mn:], lwork-mn)
|
||||
scllen = n
|
||||
} else {
|
||||
impl.Dormlq(blas.Left, blas.NoTrans, n, nrhs, m,
|
||||
a, lda,
|
||||
work,
|
||||
b, ldb,
|
||||
work[mn:], lwork-mn)
|
||||
ok := impl.Dtrtrs(blas.Lower, blas.Trans, blas.NonUnit,
|
||||
m, nrhs,
|
||||
a, lda,
|
||||
b, ldb)
|
||||
if !ok {
|
||||
return false
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
// Adjust answer vector based on scaling.
|
||||
if iascl == 1 {
|
||||
impl.Dlascl(lapack.General, 0, 0, anrm, smlnum, scllen, nrhs, b, ldb)
|
||||
}
|
||||
if iascl == 2 {
|
||||
impl.Dlascl(lapack.General, 0, 0, anrm, bignum, scllen, nrhs, b, ldb)
|
||||
}
|
||||
if ibscl == 1 {
|
||||
impl.Dlascl(lapack.General, 0, 0, smlnum, bnrm, scllen, nrhs, b, ldb)
|
||||
}
|
||||
if ibscl == 2 {
|
||||
impl.Dlascl(lapack.General, 0, 0, bignum, bnrm, scllen, nrhs, b, ldb)
|
||||
}
|
||||
return true
|
||||
}
|
||||
+57
@@ -0,0 +1,57 @@
|
||||
// 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 native
|
||||
|
||||
import "github.com/gonum/blas"
|
||||
|
||||
// Dgeqr2 computes a QR factorization of the m×n matrix a.
|
||||
//
|
||||
// In a QR factorization, Q is an m×m orthonormal matrix, and R is an
|
||||
// upper triangular m×n matrix.
|
||||
//
|
||||
// During Dgeqr2, a is modified to contain the information to construct Q and R.
|
||||
// The upper triangle of a contains the matrix R. The lower triangular elements
|
||||
// (not including the diagonal) contain the elementary reflectors. Tau is modified
|
||||
// to contain the reflector scales. Tau must have length at least k = min(m,n), and
|
||||
// this function will panic otherwise.
|
||||
//
|
||||
// The ith elementary reflector can be explicitly constructed by first extracting
|
||||
// the
|
||||
// v[j] = 0 j < i
|
||||
// v[j] = i j == i
|
||||
// v[j] = a[i*lda+j] j > i
|
||||
// and computing h_i = I - tau[i] * v * v^T.
|
||||
//
|
||||
// The orthonormal matrix Q can be constucted from a product of these elementary
|
||||
// reflectors, Q = H_1*H_2 ... H_k, where k = min(m,n).
|
||||
//
|
||||
// Work is temporary storage of length at least n and this function will panic otherwise.
|
||||
func (impl Implementation) Dgeqr2(m, n int, a []float64, lda int, tau, work []float64) {
|
||||
// TODO(btracey): This is oriented such that columns of a are eliminated.
|
||||
// This likely could be re-arranged to take better advantage of row-major
|
||||
// storage.
|
||||
checkMatrix(m, n, a, lda)
|
||||
if len(work) < n {
|
||||
panic(badWork)
|
||||
}
|
||||
k := min(m, n)
|
||||
if len(tau) < k {
|
||||
panic(badTau)
|
||||
}
|
||||
for i := 0; i < k; i++ {
|
||||
// Generate elementary reflector H(i).
|
||||
a[i*lda+i], tau[i] = impl.Dlarfg(m-i, a[i*lda+i], a[min((i+1), m-1)*lda+i:], lda)
|
||||
if i < n-1 {
|
||||
aii := a[i*lda+i]
|
||||
a[i*lda+i] = 1
|
||||
impl.Dlarf(blas.Left, m-i, n-i-1,
|
||||
a[i*lda+i:], lda,
|
||||
tau[i],
|
||||
a[i*lda+i+1:], lda,
|
||||
work)
|
||||
a[i*lda+i] = aii
|
||||
}
|
||||
}
|
||||
}
|
||||
+98
@@ -0,0 +1,98 @@
|
||||
// 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 native
|
||||
|
||||
import (
|
||||
"github.com/gonum/blas"
|
||||
"github.com/gonum/lapack"
|
||||
)
|
||||
|
||||
// Dgeqrf computes the QR factorization of the m×n matrix a using a blocked
|
||||
// algorithm. Please see the documentation for Dgeqr2 for a description of the
|
||||
// parameters at entry and exit.
|
||||
//
|
||||
// Work is temporary storage, and lwork specifies the usable memory length.
|
||||
// At minimum, lwork >= m and this function will panic otherwise.
|
||||
// Dgeqrf is a blocked LQ factorization, but the block size is limited
|
||||
// by the temporary space available. If lwork == -1, instead of performing Dgelqf,
|
||||
// the optimal work length will be stored into work[0].
|
||||
//
|
||||
// tau must be at least len min(m,n), and this function will panic otherwise.
|
||||
func (impl Implementation) Dgeqrf(m, n int, a []float64, lda int, tau, work []float64, lwork int) {
|
||||
// TODO(btracey): This algorithm is oriented for column-major storage.
|
||||
// Consider modifying the algorithm to better suit row-major storage.
|
||||
|
||||
// nb is the optimal blocksize, i.e. the number of columns transformed at a time.
|
||||
nb := impl.Ilaenv(1, "DGEQRF", " ", m, n, -1, -1)
|
||||
lworkopt := n * max(nb, 1)
|
||||
lworkopt = max(n, lworkopt)
|
||||
if lwork == -1 {
|
||||
work[0] = float64(lworkopt)
|
||||
return
|
||||
}
|
||||
checkMatrix(m, n, a, lda)
|
||||
if len(work) < lwork {
|
||||
panic(shortWork)
|
||||
}
|
||||
if lwork < n {
|
||||
panic(badWork)
|
||||
}
|
||||
k := min(m, n)
|
||||
if len(tau) < k {
|
||||
panic(badTau)
|
||||
}
|
||||
if k == 0 {
|
||||
return
|
||||
}
|
||||
nbmin := 2 // Minimal number of blocks
|
||||
var nx int // Use unblocked (unless changed in the next for loop)
|
||||
iws := n
|
||||
ldwork := nb
|
||||
// Only consider blocked if the suggested number of blocks is > 1 and the
|
||||
// number of columns is sufficiently large.
|
||||
if nb > 1 && k > nb {
|
||||
// nx is the crossover point. Above this value the blocked routine should be used.
|
||||
nx = max(0, impl.Ilaenv(3, "DGEQRF", " ", m, n, -1, -1))
|
||||
if k > nx {
|
||||
iws = ldwork * n
|
||||
if lwork < iws {
|
||||
// Not enough workspace to use the optimal number of blocks. Instead,
|
||||
// get the maximum allowable number of blocks.
|
||||
nb = lwork / n
|
||||
nbmin = max(2, impl.Ilaenv(2, "DGEQRF", " ", m, n, -1, -1))
|
||||
}
|
||||
}
|
||||
}
|
||||
for i := range work {
|
||||
work[i] = 0
|
||||
}
|
||||
// Compute QR using a blocked algorithm.
|
||||
var i int
|
||||
if nb >= nbmin && nb < k && nx < k {
|
||||
for i = 0; i < k-nx; i += nb {
|
||||
ib := min(k-i, nb)
|
||||
// Compute the QR factorization of the current block.
|
||||
impl.Dgeqr2(m-i, ib, a[i*lda+i:], lda, tau[i:], work)
|
||||
if i+ib < n {
|
||||
// Form the triangular factor of the block reflector and apply H^T
|
||||
// In Dlarft, work becomes the T matrix.
|
||||
impl.Dlarft(lapack.Forward, lapack.ColumnWise, m-i, ib,
|
||||
a[i*lda+i:], lda,
|
||||
tau[i:],
|
||||
work, ldwork)
|
||||
impl.Dlarfb(blas.Left, blas.Trans, lapack.Forward, lapack.ColumnWise,
|
||||
m-i, n-i-ib, ib,
|
||||
a[i*lda+i:], lda,
|
||||
work, ldwork,
|
||||
a[i*lda+i+ib:], lda,
|
||||
work[ib*ldwork:], ldwork)
|
||||
}
|
||||
}
|
||||
}
|
||||
// Call unblocked code on the remaining columns.
|
||||
if i < k {
|
||||
impl.Dgeqr2(m-i, n-i, a[i*lda+i:], lda, tau[i:], work)
|
||||
}
|
||||
}
|
||||
+76
@@ -0,0 +1,76 @@
|
||||
// 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 native
|
||||
|
||||
import (
|
||||
"math"
|
||||
|
||||
"github.com/gonum/lapack"
|
||||
)
|
||||
|
||||
// Dlange computes the matrix norm of the general m×n matrix a. The input norm
|
||||
// specifies the norm computed.
|
||||
// lapack.MaxAbs: the maximum absolute value of an element.
|
||||
// lapack.MaxColumnSum: the maximum column sum of the absolute values of the entries.
|
||||
// lapack.MaxRowSum: the maximum row sum of the absolute values of the entries.
|
||||
// lapack.Frobenius: the square root of the sum of the squares of the entries.
|
||||
// If norm == lapack.MaxColumnSum, work must be of length n, and this function will panic otherwise.
|
||||
// There are no restrictions on work for the other matrix norms.
|
||||
func (impl Implementation) Dlange(norm lapack.MatrixNorm, m, n int, a []float64, lda int, work []float64) float64 {
|
||||
// TODO(btracey): These should probably be refactored to use BLAS calls.
|
||||
checkMatrix(m, n, a, lda)
|
||||
if m == 0 && n == 0 {
|
||||
return 0
|
||||
}
|
||||
if norm == lapack.MaxAbs {
|
||||
var value float64
|
||||
for i := 0; i < m; i++ {
|
||||
for j := 0; j < n; j++ {
|
||||
value = math.Max(value, math.Abs(a[i*lda+j]))
|
||||
}
|
||||
}
|
||||
return value
|
||||
}
|
||||
if norm == lapack.MaxColumnSum {
|
||||
if len(work) < n {
|
||||
panic(badWork)
|
||||
}
|
||||
for i := 0; i < n; i++ {
|
||||
work[i] = 0
|
||||
}
|
||||
for i := 0; i < m; i++ {
|
||||
for j := 0; j < n; j++ {
|
||||
work[j] += math.Abs(a[i*lda+j])
|
||||
}
|
||||
}
|
||||
var value float64
|
||||
for i := 0; i < n; i++ {
|
||||
value = math.Max(value, work[i])
|
||||
}
|
||||
return value
|
||||
}
|
||||
if norm == lapack.MaxRowSum {
|
||||
var value float64
|
||||
for i := 0; i < m; i++ {
|
||||
var sum float64
|
||||
for j := 0; j < n; j++ {
|
||||
sum += math.Abs(a[i*lda+j])
|
||||
}
|
||||
value = math.Max(value, sum)
|
||||
}
|
||||
return value
|
||||
}
|
||||
if norm == lapack.NormFrob {
|
||||
var value float64
|
||||
scale := 0.0
|
||||
sum := 1.0
|
||||
for i := 0; i < m; i++ {
|
||||
scale, sum = impl.Dlassq(n, a[i*lda:], 1, scale, sum)
|
||||
}
|
||||
value = scale * math.Sqrt(sum)
|
||||
return value
|
||||
}
|
||||
panic("lapack: bad matrix norm")
|
||||
}
|
||||
+12
@@ -0,0 +1,12 @@
|
||||
// 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 native
|
||||
|
||||
import "math"
|
||||
|
||||
// Dlapy2 is the LAPACK version of math.Hypot.
|
||||
func (Implementation) Dlapy2(x, y float64) float64 {
|
||||
return math.Hypot(x, y)
|
||||
}
|
||||
+82
@@ -0,0 +1,82 @@
|
||||
// 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 native
|
||||
|
||||
import (
|
||||
"github.com/gonum/blas"
|
||||
"github.com/gonum/blas/blas64"
|
||||
)
|
||||
|
||||
// Dlarf applies an elementary reflector to a general rectangular matrix c.
|
||||
// This computes
|
||||
// c = h * c if side == Left
|
||||
// c = c * h if side == right
|
||||
// where
|
||||
// h = 1 - tau * v * v^T
|
||||
// and c is an m * n matrix.
|
||||
//
|
||||
|
||||
// Work is temporary storage of length at least m if side == Left and at least
|
||||
// n if side == Right. This function will panic if this length requirement is not met.
|
||||
func (impl Implementation) Dlarf(side blas.Side, m, n int, v []float64, incv int, tau float64, c []float64, ldc int, work []float64) {
|
||||
applyleft := side == blas.Left
|
||||
if (applyleft && len(work) < n) || (!applyleft && len(work) < m) {
|
||||
panic(badWork)
|
||||
}
|
||||
checkMatrix(m, n, c, ldc)
|
||||
|
||||
// v has length m if applyleft and n otherwise.
|
||||
lenV := n
|
||||
if applyleft {
|
||||
lenV = m
|
||||
}
|
||||
|
||||
checkVector(lenV, v, incv)
|
||||
|
||||
lastv := 0 // last non-zero element of v
|
||||
lastc := 0 // last non-zero row/column of c
|
||||
if tau != 0 {
|
||||
var i int
|
||||
if applyleft {
|
||||
lastv = m - 1
|
||||
} else {
|
||||
lastv = n - 1
|
||||
}
|
||||
if incv > 0 {
|
||||
i = lastv * incv
|
||||
}
|
||||
|
||||
// Look for the last non-zero row in v.
|
||||
for lastv >= 0 && v[i] == 0 {
|
||||
lastv--
|
||||
i -= incv
|
||||
}
|
||||
if applyleft {
|
||||
// Scan for the last non-zero column in C[0:lastv, :]
|
||||
lastc = impl.Iladlc(lastv+1, n, c, ldc)
|
||||
} else {
|
||||
// Scan for the last non-zero row in C[:, 0:lastv]
|
||||
lastc = impl.Iladlr(m, lastv+1, c, ldc)
|
||||
}
|
||||
}
|
||||
if lastv == -1 || lastc == -1 {
|
||||
return
|
||||
}
|
||||
// Sometimes 1-indexing is nicer ...
|
||||
bi := blas64.Implementation()
|
||||
if applyleft {
|
||||
// Form H * C
|
||||
// w[0:lastc+1] = c[1:lastv+1, 1:lastc+1]^T * v[1:lastv+1,1]
|
||||
bi.Dgemv(blas.Trans, lastv+1, lastc+1, 1, c, ldc, v, incv, 0, work, 1)
|
||||
// c[0: lastv, 0: lastc] = c[...] - w[0:lastv, 1] * v[1:lastc, 1]^T
|
||||
bi.Dger(lastv+1, lastc+1, -tau, v, incv, work, 1, c, ldc)
|
||||
return
|
||||
}
|
||||
// Form C*H
|
||||
// w[0:lastc+1,1] := c[0:lastc+1,0:lastv+1] * v[0:lastv+1,1]
|
||||
bi.Dgemv(blas.NoTrans, lastc+1, lastv+1, 1, c, ldc, v, incv, 0, work, 1)
|
||||
// c[0:lastc+1,0:lastv+1] = c[...] - w[0:lastc+1,0] * v[0:lastv+1,0]^T
|
||||
bi.Dger(lastc+1, lastv+1, -tau, work, 1, v, incv, c, ldc)
|
||||
}
|
||||
+424
@@ -0,0 +1,424 @@
|
||||
// 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 native
|
||||
|
||||
import (
|
||||
"github.com/gonum/blas"
|
||||
"github.com/gonum/blas/blas64"
|
||||
"github.com/gonum/lapack"
|
||||
)
|
||||
|
||||
// Dlarfb applies a block reflector to a matrix.
|
||||
//
|
||||
// In the call to Dlarfb, the mxn c is multiplied by the implicitly defined matrix h as follows:
|
||||
// c = h * c if side == Left and trans == NoTrans
|
||||
// c = c * h if side == Right and trans == NoTrans
|
||||
// c = h^T * c if side == Left and trans == Trans
|
||||
// c = c * h^t if side == Right and trans == Trans
|
||||
// h is a product of elementary reflectors. direct sets the direction of multiplication
|
||||
// h = h_1 * h_2 * ... * h_k if direct == Forward
|
||||
// h = h_k * h_k-1 * ... * h_1 if direct == Backward
|
||||
// The combination of direct and store defines the orientation of the elementary
|
||||
// reflectors. In all cases the ones on the diagonal are implicitly represented.
|
||||
//
|
||||
// If direct == lapack.Forward and store == lapack.ColumnWise
|
||||
// V = ( 1 )
|
||||
// ( v1 1 )
|
||||
// ( v1 v2 1 )
|
||||
// ( v1 v2 v3 )
|
||||
// ( v1 v2 v3 )
|
||||
// If direct == lapack.Forward and store == lapack.RowWise
|
||||
// V = ( 1 v1 v1 v1 v1 )
|
||||
// ( 1 v2 v2 v2 )
|
||||
// ( 1 v3 v3 )
|
||||
// If direct == lapack.Backward and store == lapack.ColumnWise
|
||||
// V = ( v1 v2 v3 )
|
||||
// ( v1 v2 v3 )
|
||||
// ( 1 v2 v3 )
|
||||
// ( 1 v3 )
|
||||
// ( 1 )
|
||||
// If direct == lapack.Backward and store == lapack.RowWise
|
||||
// V = ( v1 v1 1 )
|
||||
// ( v2 v2 v2 1 )
|
||||
// ( v3 v3 v3 v3 1 )
|
||||
// An elementary reflector can be explicitly constructed by extracting the
|
||||
// corresponding elements of v, placing a 1 where the diagonal would be, and
|
||||
// placing zeros in the remaining elements.
|
||||
//
|
||||
// t is a k×k matrix containing the block reflector, and this function will panic
|
||||
// if t is not of sufficient size. See Dlarft for more information.
|
||||
//
|
||||
// Work is a temporary storage matrix with stride ldwork.
|
||||
// Work must be of size at least n×k side == Left and m×k if side == Right, and
|
||||
// this function will panic if this size is not met.
|
||||
func (Implementation) Dlarfb(side blas.Side, trans blas.Transpose, direct lapack.Direct,
|
||||
store lapack.StoreV, m, n, k int, v []float64, ldv int, t []float64, ldt int,
|
||||
c []float64, ldc int, work []float64, ldwork int) {
|
||||
|
||||
checkMatrix(m, n, c, ldc)
|
||||
if m == 0 || n == 0 {
|
||||
return
|
||||
}
|
||||
if k < 0 {
|
||||
panic("lapack: negative number of transforms")
|
||||
}
|
||||
if side != blas.Left && side != blas.Right {
|
||||
panic(badSide)
|
||||
}
|
||||
if trans != blas.Trans && trans != blas.NoTrans {
|
||||
panic(badTrans)
|
||||
}
|
||||
if direct != lapack.Forward && direct != lapack.Backward {
|
||||
panic(badDirect)
|
||||
}
|
||||
if store != lapack.ColumnWise && store != lapack.RowWise {
|
||||
panic(badStore)
|
||||
}
|
||||
|
||||
rowsWork := n
|
||||
if side == blas.Right {
|
||||
rowsWork = m
|
||||
}
|
||||
checkMatrix(rowsWork, k, work, ldwork)
|
||||
|
||||
bi := blas64.Implementation()
|
||||
|
||||
transt := blas.Trans
|
||||
if trans == blas.Trans {
|
||||
transt = blas.NoTrans
|
||||
}
|
||||
// TODO(btracey): This follows the original Lapack code where the
|
||||
// elements are copied into the columns of the working array. The
|
||||
// loops should go in the other direction so the data is written
|
||||
// into the rows of work so the copy is not strided. A bigger change
|
||||
// would be to replace work with work^T, but benchmarks would be
|
||||
// needed to see if the change is merited.
|
||||
if store == lapack.ColumnWise {
|
||||
if direct == lapack.Forward {
|
||||
// V1 is the first k rows of C. V2 is the remaining rows.
|
||||
if side == blas.Left {
|
||||
// W = C^T V = C1^T V1 + C2^T V2 (stored in work).
|
||||
|
||||
// W = C1.
|
||||
for j := 0; j < k; j++ {
|
||||
bi.Dcopy(n, c[j*ldc:], 1, work[j:], ldwork)
|
||||
}
|
||||
// W = W * V1.
|
||||
bi.Dtrmm(blas.Right, blas.Lower, blas.NoTrans, blas.Unit,
|
||||
n, k, 1,
|
||||
v, ldv,
|
||||
work, ldwork)
|
||||
if m > k {
|
||||
// W = W + C2^T V2.
|
||||
bi.Dgemm(blas.Trans, blas.NoTrans, n, k, m-k,
|
||||
1, c[k*ldc:], ldc, v[k*ldv:], ldv,
|
||||
1, work, ldwork)
|
||||
}
|
||||
// W = W * T^T or W * T.
|
||||
bi.Dtrmm(blas.Right, blas.Upper, transt, blas.NonUnit, n, k,
|
||||
1, t, ldt,
|
||||
work, ldwork)
|
||||
// C -= V * W^T.
|
||||
if m > k {
|
||||
// C2 -= V2 * W^T.
|
||||
bi.Dgemm(blas.NoTrans, blas.Trans, m-k, n, k,
|
||||
-1, v[k*ldv:], ldv, work, ldwork,
|
||||
1, c[k*ldc:], ldc)
|
||||
}
|
||||
// W *= V1^T.
|
||||
bi.Dtrmm(blas.Right, blas.Lower, blas.Trans, blas.Unit, n, k,
|
||||
1, v, ldv,
|
||||
work, ldwork)
|
||||
// C1 -= W^T.
|
||||
// TODO(btracey): This should use blas.Axpy.
|
||||
for i := 0; i < n; i++ {
|
||||
for j := 0; j < k; j++ {
|
||||
c[j*ldc+i] -= work[i*ldwork+j]
|
||||
}
|
||||
}
|
||||
return
|
||||
}
|
||||
// Form C = C * H or C * H^T, where C = (C1 C2).
|
||||
|
||||
// W = C1.
|
||||
for i := 0; i < k; i++ {
|
||||
bi.Dcopy(m, c[i:], ldc, work[i:], ldwork)
|
||||
}
|
||||
// W *= V1.
|
||||
bi.Dtrmm(blas.Right, blas.Lower, blas.NoTrans, blas.Unit, m, k,
|
||||
1, v, ldv,
|
||||
work, ldwork)
|
||||
if n > k {
|
||||
bi.Dgemm(blas.NoTrans, blas.NoTrans, m, k, n-k,
|
||||
1, c[k:], ldc, v[k*ldv:], ldv,
|
||||
1, work, ldwork)
|
||||
}
|
||||
// W *= T or T^T.
|
||||
bi.Dtrmm(blas.Right, blas.Upper, trans, blas.NonUnit, m, k,
|
||||
1, t, ldt,
|
||||
work, ldwork)
|
||||
if n > k {
|
||||
bi.Dgemm(blas.NoTrans, blas.Trans, m, n-k, k,
|
||||
-1, work, ldwork, v[k*ldv:], ldv,
|
||||
1, c[k:], ldc)
|
||||
}
|
||||
// C -= W * V^T.
|
||||
bi.Dtrmm(blas.Right, blas.Lower, blas.Trans, blas.Unit, m, k,
|
||||
1, v, ldv,
|
||||
work, ldwork)
|
||||
// C -= W.
|
||||
// TODO(btracey): This should use blas.Axpy.
|
||||
for i := 0; i < m; i++ {
|
||||
for j := 0; j < k; j++ {
|
||||
c[i*ldc+j] -= work[i*ldwork+j]
|
||||
}
|
||||
}
|
||||
return
|
||||
}
|
||||
// V = (V1)
|
||||
// = (V2) (last k rows)
|
||||
// Where V2 is unit upper triangular.
|
||||
if side == blas.Left {
|
||||
// Form H * C or
|
||||
// W = C^T V.
|
||||
|
||||
// W = C2^T.
|
||||
for j := 0; j < k; j++ {
|
||||
bi.Dcopy(n, c[(m-k+j)*ldc:], 1, work[j:], ldwork)
|
||||
}
|
||||
// W *= V2.
|
||||
bi.Dtrmm(blas.Right, blas.Upper, blas.NoTrans, blas.Unit, n, k,
|
||||
1, v[(m-k)*ldv:], ldv,
|
||||
work, ldwork)
|
||||
if m > k {
|
||||
// W += C1^T * V1.
|
||||
bi.Dgemm(blas.Trans, blas.NoTrans, n, k, m-k,
|
||||
1, c, ldc, v, ldv,
|
||||
1, work, ldwork)
|
||||
}
|
||||
// W *= T or T^T.
|
||||
bi.Dtrmm(blas.Right, blas.Lower, transt, blas.NonUnit, n, k,
|
||||
1, t, ldt,
|
||||
work, ldwork)
|
||||
// C -= V * W^T.
|
||||
if m > k {
|
||||
bi.Dgemm(blas.NoTrans, blas.Trans, m-k, n, k,
|
||||
-1, v, ldv, work, ldwork,
|
||||
1, c, ldc)
|
||||
}
|
||||
// W *= V2^T.
|
||||
bi.Dtrmm(blas.Right, blas.Upper, blas.Trans, blas.Unit, n, k,
|
||||
1, v[(m-k)*ldv:], ldv,
|
||||
work, ldwork)
|
||||
// C2 -= W^T.
|
||||
// TODO(btracey): This should use blas.Axpy.
|
||||
for i := 0; i < n; i++ {
|
||||
for j := 0; j < k; j++ {
|
||||
c[(m-k+j)*ldc+i] -= work[i*ldwork+j]
|
||||
}
|
||||
}
|
||||
return
|
||||
}
|
||||
// Form C * H or C * H^T where C = (C1 C2).
|
||||
// W = C * V.
|
||||
|
||||
// W = C2.
|
||||
for j := 0; j < k; j++ {
|
||||
bi.Dcopy(m, c[n-k+j:], ldc, work[j:], ldwork)
|
||||
}
|
||||
|
||||
// W = W * V2.
|
||||
bi.Dtrmm(blas.Right, blas.Upper, blas.NoTrans, blas.Unit, m, k,
|
||||
1, v[(n-k)*ldv:], ldv,
|
||||
work, ldwork)
|
||||
if n > k {
|
||||
bi.Dgemm(blas.NoTrans, blas.NoTrans, m, k, n-k,
|
||||
1, c, ldc, v, ldv,
|
||||
1, work, ldwork)
|
||||
}
|
||||
// W *= T or T^T.
|
||||
bi.Dtrmm(blas.Right, blas.Lower, trans, blas.NonUnit, m, k,
|
||||
1, t, ldt,
|
||||
work, ldwork)
|
||||
// C -= W * V^T.
|
||||
if n > k {
|
||||
// C1 -= W * V1^T.
|
||||
bi.Dgemm(blas.NoTrans, blas.Trans, m, n-k, k,
|
||||
-1, work, ldwork, v, ldv,
|
||||
1, c, ldc)
|
||||
}
|
||||
// W *= V2^T.
|
||||
bi.Dtrmm(blas.Right, blas.Upper, blas.Trans, blas.Unit, m, k,
|
||||
1, v[(n-k)*ldv:], ldv,
|
||||
work, ldwork)
|
||||
// C2 -= W.
|
||||
// TODO(btracey): This should use blas.Axpy.
|
||||
for i := 0; i < m; i++ {
|
||||
for j := 0; j < k; j++ {
|
||||
c[i*ldc+n-k+j] -= work[i*ldwork+j]
|
||||
}
|
||||
}
|
||||
return
|
||||
}
|
||||
// Store = Rowwise.
|
||||
if direct == lapack.Forward {
|
||||
// V = (V1 V2) where v1 is unit upper triangular.
|
||||
if side == blas.Left {
|
||||
// Form H * C or H^T * C where C = (C1; C2).
|
||||
// W = C^T * V^T.
|
||||
|
||||
// W = C1^T.
|
||||
for j := 0; j < k; j++ {
|
||||
bi.Dcopy(n, c[j*ldc:], 1, work[j:], ldwork)
|
||||
}
|
||||
// W *= V1^T.
|
||||
bi.Dtrmm(blas.Right, blas.Upper, blas.Trans, blas.Unit, n, k,
|
||||
1, v, ldv,
|
||||
work, ldwork)
|
||||
if m > k {
|
||||
bi.Dgemm(blas.Trans, blas.Trans, n, k, m-k,
|
||||
1, c[k*ldc:], ldc, v[k:], ldv,
|
||||
1, work, ldwork)
|
||||
}
|
||||
// W *= T or T^T.
|
||||
bi.Dtrmm(blas.Right, blas.Upper, transt, blas.NonUnit, n, k,
|
||||
1, t, ldt,
|
||||
work, ldwork)
|
||||
// C -= V^T * W^T.
|
||||
if m > k {
|
||||
bi.Dgemm(blas.Trans, blas.Trans, m-k, n, k,
|
||||
-1, v[k:], ldv, work, ldwork,
|
||||
1, c[k*ldc:], ldc)
|
||||
}
|
||||
// W *= V1.
|
||||
bi.Dtrmm(blas.Right, blas.Upper, blas.NoTrans, blas.Unit, n, k,
|
||||
1, v, ldv,
|
||||
work, ldwork)
|
||||
// C1 -= W^T.
|
||||
// TODO(btracey): This should use blas.Axpy.
|
||||
for i := 0; i < n; i++ {
|
||||
for j := 0; j < k; j++ {
|
||||
c[j*ldc+i] -= work[i*ldwork+j]
|
||||
}
|
||||
}
|
||||
return
|
||||
}
|
||||
// Form C * H or C * H^T where C = (C1 C2).
|
||||
// W = C * V^T.
|
||||
|
||||
// W = C1.
|
||||
for j := 0; j < k; j++ {
|
||||
bi.Dcopy(m, c[j:], ldc, work[j:], ldwork)
|
||||
}
|
||||
// W *= V1^T.
|
||||
bi.Dtrmm(blas.Right, blas.Upper, blas.Trans, blas.Unit, m, k,
|
||||
1, v, ldv,
|
||||
work, ldwork)
|
||||
if n > k {
|
||||
bi.Dgemm(blas.NoTrans, blas.Trans, m, k, n-k,
|
||||
1, c[k:], ldc, v[k:], ldv,
|
||||
1, work, ldwork)
|
||||
}
|
||||
// W *= T or T^T.
|
||||
bi.Dtrmm(blas.Right, blas.Upper, trans, blas.NonUnit, m, k,
|
||||
1, t, ldt,
|
||||
work, ldwork)
|
||||
// C -= W * V.
|
||||
if n > k {
|
||||
bi.Dgemm(blas.NoTrans, blas.NoTrans, m, n-k, k,
|
||||
-1, work, ldwork, v[k:], ldv,
|
||||
1, c[k:], ldc)
|
||||
}
|
||||
// W *= V1.
|
||||
bi.Dtrmm(blas.Right, blas.Upper, blas.NoTrans, blas.Unit, m, k,
|
||||
1, v, ldv,
|
||||
work, ldwork)
|
||||
// C1 -= W.
|
||||
// TODO(btracey): This should use blas.Axpy.
|
||||
for i := 0; i < m; i++ {
|
||||
for j := 0; j < k; j++ {
|
||||
c[i*ldc+j] -= work[i*ldwork+j]
|
||||
}
|
||||
}
|
||||
return
|
||||
}
|
||||
// V = (V1 V2) where V2 is the last k columns and is lower unit triangular.
|
||||
if side == blas.Left {
|
||||
// Form H * C or H^T C where C = (C1 ; C2).
|
||||
// W = C^T * V^T.
|
||||
|
||||
// W = C2^T.
|
||||
for j := 0; j < k; j++ {
|
||||
bi.Dcopy(n, c[(m-k+j)*ldc:], 1, work[j:], ldwork)
|
||||
}
|
||||
// W *= V2^T.
|
||||
bi.Dtrmm(blas.Right, blas.Lower, blas.Trans, blas.Unit, n, k,
|
||||
1, v[m-k:], ldv,
|
||||
work, ldwork)
|
||||
if m > k {
|
||||
bi.Dgemm(blas.Trans, blas.Trans, n, k, m-k,
|
||||
1, c, ldc, v, ldv,
|
||||
1, work, ldwork)
|
||||
}
|
||||
// W *= T or T^T.
|
||||
bi.Dtrmm(blas.Right, blas.Lower, transt, blas.NonUnit, n, k,
|
||||
1, t, ldt,
|
||||
work, ldwork)
|
||||
// C -= V^T * W^T.
|
||||
if m > k {
|
||||
bi.Dgemm(blas.Trans, blas.Trans, m-k, n, k,
|
||||
-1, v, ldv, work, ldwork,
|
||||
1, c, ldc)
|
||||
}
|
||||
// W *= V2.
|
||||
bi.Dtrmm(blas.Right, blas.Lower, blas.NoTrans, blas.Unit, n, k,
|
||||
1, v[m-k:], ldv,
|
||||
work, ldwork)
|
||||
// C2 -= W^T.
|
||||
// TODO(btracey): This should use blas.Axpy.
|
||||
for i := 0; i < n; i++ {
|
||||
for j := 0; j < k; j++ {
|
||||
c[(m-k+j)*ldc+i] -= work[i*ldwork+j]
|
||||
}
|
||||
}
|
||||
return
|
||||
}
|
||||
// Form C * H or C * H^T where C = (C1 C2).
|
||||
// W = C * V^T.
|
||||
// W = C2.
|
||||
for j := 0; j < k; j++ {
|
||||
bi.Dcopy(m, c[n-k+j:], ldc, work[j:], ldwork)
|
||||
}
|
||||
// W *= V2^T.
|
||||
bi.Dtrmm(blas.Right, blas.Lower, blas.Trans, blas.Unit, m, k,
|
||||
1, v[n-k:], ldv,
|
||||
work, ldwork)
|
||||
if n > k {
|
||||
bi.Dgemm(blas.NoTrans, blas.Trans, m, k, n-k,
|
||||
1, c, ldc, v, ldv,
|
||||
1, work, ldwork)
|
||||
}
|
||||
// W *= T or T^T.
|
||||
bi.Dtrmm(blas.Right, blas.Lower, trans, blas.NonUnit, m, k,
|
||||
1, t, ldt,
|
||||
work, ldwork)
|
||||
// C -= W * V.
|
||||
if n > k {
|
||||
bi.Dgemm(blas.NoTrans, blas.NoTrans, m, n-k, k,
|
||||
-1, work, ldwork, v, ldv,
|
||||
1, c, ldc)
|
||||
}
|
||||
// W *= V2.
|
||||
bi.Dtrmm(blas.Right, blas.Lower, blas.NoTrans, blas.Unit, m, k,
|
||||
1, v[n-k:], ldv,
|
||||
work, ldwork)
|
||||
// C1 -= W.
|
||||
// TODO(btracey): This should use blas.Axpy.
|
||||
for i := 0; i < m; i++ {
|
||||
for j := 0; j < k; j++ {
|
||||
c[i*ldc+n-k+j] -= work[i*ldwork+j]
|
||||
}
|
||||
}
|
||||
}
|
||||
+60
@@ -0,0 +1,60 @@
|
||||
// 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 native
|
||||
|
||||
import (
|
||||
"math"
|
||||
|
||||
"github.com/gonum/blas/blas64"
|
||||
)
|
||||
|
||||
// Dlarfg generates an elementary reflector for a Householder matrix. It creates
|
||||
// a real elementary reflector of order n such that
|
||||
// H * (alpha) = (beta)
|
||||
// ( x) ( 0)
|
||||
// H^T * H = I
|
||||
// H is represented in the form
|
||||
// H = 1 - tau * (1; v) * (1 v^T)
|
||||
// where tau is a real scalar.
|
||||
//
|
||||
// On entry, x contains the vector x, on exit it contains v.
|
||||
func (impl Implementation) Dlarfg(n int, alpha float64, x []float64, incX int) (beta, tau float64) {
|
||||
if n < 0 {
|
||||
panic(nLT0)
|
||||
}
|
||||
if n <= 1 {
|
||||
return alpha, 0
|
||||
}
|
||||
checkVector(n-1, x, incX)
|
||||
bi := blas64.Implementation()
|
||||
xnorm := bi.Dnrm2(n-1, x, incX)
|
||||
if xnorm == 0 {
|
||||
return alpha, 0
|
||||
}
|
||||
beta = -math.Copysign(impl.Dlapy2(alpha, xnorm), alpha)
|
||||
safmin := dlamchS / dlamchE
|
||||
knt := 0
|
||||
if math.Abs(beta) < safmin {
|
||||
// xnorm and beta may be innacurate, scale x and recompute.
|
||||
rsafmn := 1 / safmin
|
||||
for {
|
||||
knt++
|
||||
bi.Dscal(n-1, rsafmn, x, incX)
|
||||
beta *= rsafmn
|
||||
alpha *= rsafmn
|
||||
if math.Abs(beta) >= safmin {
|
||||
break
|
||||
}
|
||||
}
|
||||
xnorm = bi.Dnrm2(n-1, x, incX)
|
||||
beta = -math.Copysign(impl.Dlapy2(alpha, xnorm), alpha)
|
||||
}
|
||||
tau = (beta - alpha) / beta
|
||||
bi.Dscal(n-1, 1/(alpha-beta), x, incX)
|
||||
for j := 0; j < knt; j++ {
|
||||
beta *= safmin
|
||||
}
|
||||
return beta, tau
|
||||
}
|
||||
+148
@@ -0,0 +1,148 @@
|
||||
// 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 native
|
||||
|
||||
import (
|
||||
"github.com/gonum/blas"
|
||||
"github.com/gonum/blas/blas64"
|
||||
"github.com/gonum/lapack"
|
||||
)
|
||||
|
||||
// Dlarft forms the triangular factor t of a block reflector, storing the answer
|
||||
// in t.
|
||||
// H = 1 - V * T * V^T if store == lapack.ColumnWise
|
||||
// H = 1 - V^T * T * V if store == lapack.RowWise
|
||||
// H is defined by a product of the elementary reflectors where
|
||||
// H = H_1 * H_2 * ... * H_k if direct == lapack.Forward
|
||||
// H = H_k * H_k-1 * ... * H_1 if direct == lapack.Backward
|
||||
//
|
||||
// t is a k×k triangular matrix. t is upper triangular if direct = lapack.Forward
|
||||
// and lower triangular otherwise. This function will panic if t is not of
|
||||
// sufficient size.
|
||||
//
|
||||
// store describes the storage of the elementary reflectors in v. Please see
|
||||
// Dlarfb for a description of layout.
|
||||
//
|
||||
// tau contains the scalar factor of the elementary reflectors h.
|
||||
func (Implementation) Dlarft(direct lapack.Direct, store lapack.StoreV, n, k int,
|
||||
v []float64, ldv int, tau []float64, t []float64, ldt int) {
|
||||
if n == 0 {
|
||||
return
|
||||
}
|
||||
if n < 0 || k < 0 {
|
||||
panic(negDimension)
|
||||
}
|
||||
if direct != lapack.Forward && direct != lapack.Backward {
|
||||
panic(badDirect)
|
||||
}
|
||||
if store != lapack.RowWise && store != lapack.ColumnWise {
|
||||
panic(badStore)
|
||||
}
|
||||
if len(tau) < k {
|
||||
panic(badTau)
|
||||
}
|
||||
checkMatrix(k, k, t, ldt)
|
||||
bi := blas64.Implementation()
|
||||
// TODO(btracey): There are a number of minor obvious loop optimizations here.
|
||||
// TODO(btracey): It may be possible to rearrange some of the code so that
|
||||
// index of 1 is more common in the Dgemv.
|
||||
if direct == lapack.Forward {
|
||||
prevlastv := n - 1
|
||||
for i := 0; i < k; i++ {
|
||||
prevlastv = max(i, prevlastv)
|
||||
if tau[i] == 0 {
|
||||
for j := 0; j <= i; j++ {
|
||||
t[j*ldt+i] = 0
|
||||
}
|
||||
continue
|
||||
}
|
||||
var lastv int
|
||||
if store == lapack.ColumnWise {
|
||||
// skip trailing zeros
|
||||
for lastv = n - 1; lastv >= i+1; lastv-- {
|
||||
if v[lastv*ldv+i] != 0 {
|
||||
break
|
||||
}
|
||||
}
|
||||
for j := 0; j < i; j++ {
|
||||
t[j*ldt+i] = -tau[i] * v[i*ldv+j]
|
||||
}
|
||||
j := min(lastv, prevlastv)
|
||||
bi.Dgemv(blas.Trans, j-i, i,
|
||||
-tau[i], v[(i+1)*ldv:], ldv, v[(i+1)*ldv+i:], ldv,
|
||||
1, t[i:], ldt)
|
||||
} else {
|
||||
for lastv = n - 1; lastv >= i+1; lastv-- {
|
||||
if v[i*ldv+lastv] != 0 {
|
||||
break
|
||||
}
|
||||
}
|
||||
for j := 0; j < i; j++ {
|
||||
t[j*ldt+i] = -tau[i] * v[j*ldv+i]
|
||||
}
|
||||
j := min(lastv, prevlastv)
|
||||
bi.Dgemv(blas.NoTrans, i, j-i,
|
||||
-tau[i], v[i+1:], ldv, v[i*ldv+i+1:], 1,
|
||||
1, t[i:], ldt)
|
||||
}
|
||||
bi.Dtrmv(blas.Upper, blas.NoTrans, blas.NonUnit, i, t, ldt, t[i:], ldt)
|
||||
t[i*ldt+i] = tau[i]
|
||||
if i > 1 {
|
||||
prevlastv = max(prevlastv, lastv)
|
||||
} else {
|
||||
prevlastv = lastv
|
||||
}
|
||||
}
|
||||
return
|
||||
}
|
||||
prevlastv := 0
|
||||
for i := k - 1; i >= 0; i-- {
|
||||
if tau[i] == 0 {
|
||||
for j := i; j < k; j++ {
|
||||
t[j*ldt+i] = 0
|
||||
}
|
||||
continue
|
||||
}
|
||||
var lastv int
|
||||
if i < k-1 {
|
||||
if store == lapack.ColumnWise {
|
||||
for lastv = 0; lastv < i; lastv++ {
|
||||
if v[lastv*ldv+i] != 0 {
|
||||
break
|
||||
}
|
||||
}
|
||||
for j := i + 1; j < k; j++ {
|
||||
t[j*ldt+i] = -tau[i] * v[(n-k+i)*ldv+j]
|
||||
}
|
||||
j := max(lastv, prevlastv)
|
||||
bi.Dgemv(blas.Trans, n-k+i-j, k-i-1,
|
||||
-tau[i], v[j*ldv+i+1:], ldv, v[j*ldv+i:], ldv,
|
||||
1, t[(i+1)*ldt+i:], ldt)
|
||||
} else {
|
||||
for lastv := 0; lastv < i; lastv++ {
|
||||
if v[i*ldv+lastv] != 0 {
|
||||
break
|
||||
}
|
||||
}
|
||||
for j := i + 1; j < k; j++ {
|
||||
t[j*ldt+i] = -tau[i] * v[j*ldv+n-k+i]
|
||||
}
|
||||
j := max(lastv, prevlastv)
|
||||
bi.Dgemv(blas.NoTrans, k-i-1, n-k+i-j,
|
||||
-tau[i], v[(i+1)*ldv+j:], ldv, v[i*ldv+j:], 1,
|
||||
1, t[(i+1)*ldt+i:], ldt)
|
||||
}
|
||||
bi.Dtrmv(blas.Lower, blas.NoTrans, blas.NonUnit, k-i-1,
|
||||
t[(i+1)*ldt+i+1:], ldt,
|
||||
t[(i+1)*ldt+i:], ldt)
|
||||
if i > 0 {
|
||||
prevlastv = min(prevlastv, lastv)
|
||||
} else {
|
||||
prevlastv = lastv
|
||||
}
|
||||
}
|
||||
t[i*ldt+i] = tau[i]
|
||||
}
|
||||
}
|
||||
+72
@@ -0,0 +1,72 @@
|
||||
// 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 native
|
||||
|
||||
import (
|
||||
"math"
|
||||
|
||||
"github.com/gonum/lapack"
|
||||
)
|
||||
|
||||
// Dlascl multiplies a rectangular matrix by a scalar.
|
||||
func (impl Implementation) Dlascl(kind lapack.MatrixType, kl, ku int, cfrom, cto float64, m, n int, a []float64, lda int) {
|
||||
checkMatrix(m, n, a, lda)
|
||||
if cfrom == 0 {
|
||||
panic("dlascl: zero divisor")
|
||||
}
|
||||
if math.IsNaN(cfrom) || math.IsNaN(cto) {
|
||||
panic("dlascl: NaN scale factor")
|
||||
}
|
||||
if n == 0 || m == 0 {
|
||||
return
|
||||
}
|
||||
smlnum := dlamchS
|
||||
bignum := 1 / smlnum
|
||||
cfromc := cfrom
|
||||
ctoc := cto
|
||||
cfrom1 := cfromc * smlnum
|
||||
for {
|
||||
var done bool
|
||||
var mul, ctol float64
|
||||
if cfrom1 == cfromc {
|
||||
// cfromc is inf
|
||||
mul = ctoc / cfromc
|
||||
done = true
|
||||
ctol = ctoc
|
||||
} else {
|
||||
ctol = ctoc / bignum
|
||||
if ctol == ctoc {
|
||||
// ctoc is either 0 or inf.
|
||||
mul = ctoc
|
||||
done = true
|
||||
cfromc = 1
|
||||
} else if math.Abs(cfrom1) > math.Abs(ctoc) && ctoc != 0 {
|
||||
mul = smlnum
|
||||
done = false
|
||||
cfromc = cfrom1
|
||||
} else if math.Abs(ctol) > math.Abs(cfromc) {
|
||||
mul = bignum
|
||||
done = false
|
||||
ctoc = ctol
|
||||
} else {
|
||||
mul = ctoc / cfromc
|
||||
done = true
|
||||
}
|
||||
}
|
||||
switch kind {
|
||||
default:
|
||||
panic("lapack: not implemented")
|
||||
case lapack.General:
|
||||
for i := 0; i < m; i++ {
|
||||
for j := 0; j < n; j++ {
|
||||
a[i*lda+j] = a[i*lda+j] * mul
|
||||
}
|
||||
}
|
||||
}
|
||||
if done {
|
||||
break
|
||||
}
|
||||
}
|
||||
}
|
||||
+37
@@ -0,0 +1,37 @@
|
||||
// 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 native
|
||||
|
||||
import "github.com/gonum/blas"
|
||||
|
||||
// Dlaset sets the off-diagonal elements of a to alpha, and the diagonal elements
|
||||
// of a to beta. If uplo == blas.Upper, only the upper diagonal elements are set.
|
||||
// If uplo == blas.Lower, only the lower diagonal elements are set. If uplo is
|
||||
// otherwise, all of the elements of a are set.
|
||||
func (impl Implementation) Dlaset(uplo blas.Uplo, m, n int, alpha, beta float64, a []float64, lda int) {
|
||||
checkMatrix(m, n, a, lda)
|
||||
if uplo == blas.Upper {
|
||||
for i := 0; i < m; i++ {
|
||||
for j := i + 1; j < n; j++ {
|
||||
a[i*lda+j] = alpha
|
||||
}
|
||||
}
|
||||
} else if uplo == blas.Lower {
|
||||
for i := 0; i < m; i++ {
|
||||
for j := 0; j < i; j++ {
|
||||
a[i*lda+j] = alpha
|
||||
}
|
||||
}
|
||||
} else {
|
||||
for i := 0; i < m; i++ {
|
||||
for j := 0; j < n; j++ {
|
||||
a[i*lda+j] = alpha
|
||||
}
|
||||
}
|
||||
}
|
||||
for i := 0; i < min(m, n); i++ {
|
||||
a[i*lda+i] = beta
|
||||
}
|
||||
}
|
||||
+29
@@ -0,0 +1,29 @@
|
||||
// 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 native
|
||||
|
||||
import "math"
|
||||
|
||||
// Dlassq updates a sum of squares in scaled form. The input parameters scale and
|
||||
// sumsq represent the current scale and total sum of squares. These values are
|
||||
// updated with the information in the first n elements of the vector specified
|
||||
// by x and incX.
|
||||
func (impl Implementation) Dlassq(n int, x []float64, incx int, scale float64, sumsq float64) (scl, smsq float64) {
|
||||
if n <= 0 {
|
||||
return scale, sumsq
|
||||
}
|
||||
for ix := 0; ix <= (n-1)*incx; ix += incx {
|
||||
absxi := math.Abs(x[ix])
|
||||
if absxi > 0 || math.IsNaN(absxi) {
|
||||
if scale < absxi {
|
||||
sumsq = 1 + sumsq*(scale/absxi)*(scale/absxi)
|
||||
scale = absxi
|
||||
} else {
|
||||
sumsq += (absxi / scale) * (absxi / scale)
|
||||
}
|
||||
}
|
||||
}
|
||||
return scale, sumsq
|
||||
}
|
||||
+28
@@ -0,0 +1,28 @@
|
||||
// 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 native is a pure-go implementation of the LAPACK API. The LAPACK API defines
|
||||
// a set of algorithms for advanced matrix operations.
|
||||
//
|
||||
// The function definitions and implementations follow that of the netlib reference
|
||||
// implementation. Please see http://www.netlib.org/lapack/explore-html/ for more
|
||||
// information, and http://www.netlib.org/lapack/explore-html/d4/de1/_l_i_c_e_n_s_e_source.html
|
||||
// for more license information.
|
||||
//
|
||||
// Slice function arguments frequently represent vectors and matrices. The data
|
||||
// layout is identical to that found in https://godoc.org/github.com/gonum/blas/native.
|
||||
//
|
||||
// Most LAPACK functions are built on top the routines defined in the BLAS API,
|
||||
// and as such the computation time for many LAPACK functions is
|
||||
// dominated by BLAS calls. Here, BLAS is accessed through the
|
||||
// the blas64 package (https://godoc.org/github.com/gonum/blas/blas64). In particular,
|
||||
// this implies that an external BLAS library will be used if it is
|
||||
// registered in blas64.
|
||||
//
|
||||
// The full LAPACK capability has not been implemented at present. The full
|
||||
// API is very large, containing approximately 200 functions for double precision
|
||||
// alone. Future additions will be focused on supporting the gonum matrix
|
||||
// package (https://godoc.org/github.com/gonum/matrix/mat64), though pull requests
|
||||
// with implementations and tests for LAPACK function are encouraged.
|
||||
package native
|
||||
+86
@@ -0,0 +1,86 @@
|
||||
// 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 native
|
||||
|
||||
import "github.com/gonum/blas"
|
||||
|
||||
// Dorm2r multiplies a general matrix c by an orthogonal matrix from a QR factorization
|
||||
// determined by Dgeqrf.
|
||||
// C = Q * C if side == blas.Left and trans == blas.NoTrans
|
||||
// C = Q^T * C if side == blas.Left and trans == blas.Trans
|
||||
// C = C * Q if side == blas.Right and trans == blas.NoTrans
|
||||
// C = C * Q^T if side == blas.Right and trans == blas.Trans
|
||||
// If side == blas.Left, a is a matrix of size m×k, and if side == blas.Right
|
||||
// a is of size n×k.
|
||||
//
|
||||
// Tau contains the householder factors and is of length at least k and this function
|
||||
// will panic otherwise.
|
||||
//
|
||||
// Work is temporary storage of length at least n if side == blas.Left
|
||||
// and at least m if side == blas.Right and this function will panic otherwise.
|
||||
func (impl Implementation) Dorm2r(side blas.Side, trans blas.Transpose, m, n, k int, a []float64, lda int, tau, c []float64, ldc int, work []float64) {
|
||||
if side != blas.Left && side != blas.Right {
|
||||
panic(badSide)
|
||||
}
|
||||
if trans != blas.Trans && trans != blas.NoTrans {
|
||||
panic(badTrans)
|
||||
}
|
||||
|
||||
left := side == blas.Left
|
||||
notran := trans == blas.NoTrans
|
||||
if left {
|
||||
// Q is m x m
|
||||
checkMatrix(m, k, a, lda)
|
||||
if len(work) < n {
|
||||
panic(badWork)
|
||||
}
|
||||
} else {
|
||||
// Q is n x n
|
||||
checkMatrix(n, k, a, lda)
|
||||
if len(work) < m {
|
||||
panic(badWork)
|
||||
}
|
||||
}
|
||||
checkMatrix(m, n, c, ldc)
|
||||
if m == 0 || n == 0 || k == 0 {
|
||||
return
|
||||
}
|
||||
if len(tau) < k {
|
||||
panic(badTau)
|
||||
}
|
||||
if left {
|
||||
if notran {
|
||||
for i := k - 1; i >= 0; i-- {
|
||||
aii := a[i*lda+i]
|
||||
a[i*lda+i] = 1
|
||||
impl.Dlarf(side, m-i, n, a[i*lda+i:], lda, tau[i], c[i*ldc:], ldc, work)
|
||||
a[i*lda+i] = aii
|
||||
}
|
||||
return
|
||||
}
|
||||
for i := 0; i < k; i++ {
|
||||
aii := a[i*lda+i]
|
||||
a[i*lda+i] = 1
|
||||
impl.Dlarf(side, m-i, n, a[i*lda+i:], lda, tau[i], c[i*ldc:], ldc, work)
|
||||
a[i*lda+i] = aii
|
||||
}
|
||||
return
|
||||
}
|
||||
if notran {
|
||||
for i := 0; i < k; i++ {
|
||||
aii := a[i*lda+i]
|
||||
a[i*lda+i] = 1
|
||||
impl.Dlarf(side, m, n-i, a[i*lda+i:], lda, tau[i], c[i:], ldc, work)
|
||||
a[i*lda+i] = aii
|
||||
}
|
||||
return
|
||||
}
|
||||
for i := k - 1; i >= 0; i-- {
|
||||
aii := a[i*lda+i]
|
||||
a[i*lda+i] = 1
|
||||
impl.Dlarf(side, m, n-i, a[i*lda+i:], lda, tau[i], c[i:], ldc, work)
|
||||
a[i*lda+i] = aii
|
||||
}
|
||||
}
|
||||
+83
@@ -0,0 +1,83 @@
|
||||
// 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 native
|
||||
|
||||
import "github.com/gonum/blas"
|
||||
|
||||
// Dorml2 multiplies a general matrix c by an orthogonal matrix from an LQ factorization
|
||||
// determined by Dgelqf.
|
||||
// C = Q * C if side == blas.Left and trans == blas.NoTrans
|
||||
// C = Q^T * C if side == blas.Left and trans == blas.Trans
|
||||
// C = C * Q if side == blas.Right and trans == blas.NoTrans
|
||||
// C = C * Q^T if side == blas.Right and trans == blas.Trans
|
||||
// If side == blas.Left, a is a matrix of side k×m, and if side == blas.Right
|
||||
// a is of size k×n.
|
||||
//
|
||||
//
|
||||
// Tau contains the householder factors and is of length at least k and this function will
|
||||
// panic otherwise.
|
||||
//
|
||||
// Work is temporary storage of length at least n if side == blas.Left
|
||||
// and at least m if side == blas.Right and this function will panic otherwise.
|
||||
func (impl Implementation) Dorml2(side blas.Side, trans blas.Transpose, m, n, k int, a []float64, lda int, tau, c []float64, ldc int, work []float64) {
|
||||
if side != blas.Left && side != blas.Right {
|
||||
panic(badSide)
|
||||
}
|
||||
if trans != blas.Trans && trans != blas.NoTrans {
|
||||
panic(badTrans)
|
||||
}
|
||||
|
||||
left := side == blas.Left
|
||||
notran := trans == blas.NoTrans
|
||||
if left {
|
||||
checkMatrix(k, m, a, lda)
|
||||
if len(work) < n {
|
||||
panic(badWork)
|
||||
}
|
||||
} else {
|
||||
checkMatrix(k, n, a, lda)
|
||||
if len(work) < m {
|
||||
panic(badWork)
|
||||
}
|
||||
}
|
||||
checkMatrix(m, n, c, ldc)
|
||||
if m == 0 || n == 0 || k == 0 {
|
||||
return
|
||||
}
|
||||
switch {
|
||||
case left && notran:
|
||||
for i := 0; i < k; i++ {
|
||||
aii := a[i*lda+i]
|
||||
a[i*lda+i] = 1
|
||||
impl.Dlarf(side, m-i, n, a[i*lda+i:], 1, tau[i], c[i*ldc:], ldc, work)
|
||||
a[i*lda+i] = aii
|
||||
}
|
||||
return
|
||||
case left && !notran:
|
||||
for i := k - 1; i >= 0; i-- {
|
||||
aii := a[i*lda+i]
|
||||
a[i*lda+i] = 1
|
||||
impl.Dlarf(side, m-i, n, a[i*lda+i:], 1, tau[i], c[i*ldc:], ldc, work)
|
||||
a[i*lda+i] = aii
|
||||
}
|
||||
return
|
||||
case !left && notran:
|
||||
for i := k - 1; i >= 0; i-- {
|
||||
aii := a[i*lda+i]
|
||||
a[i*lda+i] = 1
|
||||
impl.Dlarf(side, m, n-i, a[i*lda+i:], 1, tau[i], c[i:], ldc, work)
|
||||
a[i*lda+i] = aii
|
||||
}
|
||||
return
|
||||
case !left && !notran:
|
||||
for i := 0; i < k; i++ {
|
||||
aii := a[i*lda+i]
|
||||
a[i*lda+i] = 1
|
||||
impl.Dlarf(side, m, n-i, a[i*lda+i:], 1, tau[i], c[i:], ldc, work)
|
||||
a[i*lda+i] = aii
|
||||
}
|
||||
return
|
||||
}
|
||||
}
|
||||
+155
@@ -0,0 +1,155 @@
|
||||
// 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 native
|
||||
|
||||
import (
|
||||
"github.com/gonum/blas"
|
||||
"github.com/gonum/lapack"
|
||||
)
|
||||
|
||||
// Dormlq multiplies the matrix c by the othogonal matrix q defined by the
|
||||
// slices a and tau. A and tau are as returned from Dgelqf.
|
||||
// C = Q * C if side == blas.Left and trans == blas.NoTrans
|
||||
// C = Q^T * C if side == blas.Left and trans == blas.Trans
|
||||
// C = C * Q if side == blas.Right and trans == blas.NoTrans
|
||||
// C = C * Q^T if side == blas.Right and trans == blas.Trans
|
||||
// If side == blas.Left, a is a matrix of side k×m, and if side == blas.Right
|
||||
// a is of size k×n. This uses a blocked algorithm.
|
||||
//
|
||||
// Work is temporary storage, and lwork specifies the usable memory length.
|
||||
// At minimum, lwork >= m if side == blas.Left and lwork >= n if side == blas.Right,
|
||||
// and this function will panic otherwise.
|
||||
// Dormlq uses a block algorithm, but the block size is limited
|
||||
// by the temporary space available. If lwork == -1, instead of performing Dormlq,
|
||||
// the optimal work length will be stored into work[0].
|
||||
//
|
||||
// Tau contains the householder scales and must have length at least k, and
|
||||
// this function will panic otherwise.
|
||||
func (impl Implementation) Dormlq(side blas.Side, trans blas.Transpose, m, n, k int, a []float64, lda int, tau, c []float64, ldc int, work []float64, lwork int) {
|
||||
if side != blas.Left && side != blas.Right {
|
||||
panic(badSide)
|
||||
}
|
||||
if trans != blas.Trans && trans != blas.NoTrans {
|
||||
panic(badTrans)
|
||||
}
|
||||
left := side == blas.Left
|
||||
notran := trans == blas.NoTrans
|
||||
if left {
|
||||
checkMatrix(k, m, a, lda)
|
||||
} else {
|
||||
checkMatrix(k, n, a, lda)
|
||||
}
|
||||
checkMatrix(m, n, c, ldc)
|
||||
if len(tau) < k {
|
||||
panic(badTau)
|
||||
}
|
||||
|
||||
const nbmax = 64
|
||||
nw := n
|
||||
if !left {
|
||||
nw = m
|
||||
}
|
||||
opts := string(side) + string(trans)
|
||||
nb := min(nbmax, impl.Ilaenv(1, "DORMLQ", opts, m, n, k, -1))
|
||||
lworkopt := max(1, nw) * nb
|
||||
if lwork == -1 {
|
||||
work[0] = float64(lworkopt)
|
||||
return
|
||||
}
|
||||
if left {
|
||||
if lwork < n {
|
||||
panic(badWork)
|
||||
}
|
||||
} else {
|
||||
if lwork < m {
|
||||
panic(badWork)
|
||||
}
|
||||
}
|
||||
|
||||
if m == 0 || n == 0 || k == 0 {
|
||||
return
|
||||
}
|
||||
nbmin := 2
|
||||
|
||||
ldwork := nb
|
||||
if nb > 1 && nb < k {
|
||||
iws := nw * nb
|
||||
if lwork < iws {
|
||||
nb = lwork / nw
|
||||
nbmin = max(2, impl.Ilaenv(2, "DORMLQ", opts, m, n, k, -1))
|
||||
}
|
||||
}
|
||||
if nb < nbmin || nb >= k {
|
||||
// Call unblocked code
|
||||
impl.Dorml2(side, trans, m, n, k, a, lda, tau, c, ldc, work)
|
||||
return
|
||||
}
|
||||
ldt := nb
|
||||
t := make([]float64, nb*ldt)
|
||||
|
||||
transt := blas.NoTrans
|
||||
if notran {
|
||||
transt = blas.Trans
|
||||
}
|
||||
|
||||
switch {
|
||||
case left && notran:
|
||||
for i := 0; i < k; i += nb {
|
||||
ib := min(nb, k-i)
|
||||
impl.Dlarft(lapack.Forward, lapack.RowWise, m-i, ib,
|
||||
a[i*lda+i:], lda,
|
||||
tau[i:],
|
||||
t, ldt)
|
||||
impl.Dlarfb(side, transt, lapack.Forward, lapack.RowWise, m-i, n, ib,
|
||||
a[i*lda+i:], lda,
|
||||
t, ldt,
|
||||
c[i*ldc:], ldc,
|
||||
work, ldwork)
|
||||
}
|
||||
return
|
||||
case left && !notran:
|
||||
for i := ((k - 1) / nb) * nb; i >= 0; i -= nb {
|
||||
ib := min(nb, k-i)
|
||||
impl.Dlarft(lapack.Forward, lapack.RowWise, m-i, ib,
|
||||
a[i*lda+i:], lda,
|
||||
tau[i:],
|
||||
t, ldt)
|
||||
impl.Dlarfb(side, transt, lapack.Forward, lapack.RowWise, m-i, n, ib,
|
||||
a[i*lda+i:], lda,
|
||||
t, ldt,
|
||||
c[i*ldc:], ldc,
|
||||
work, ldwork)
|
||||
}
|
||||
return
|
||||
case !left && notran:
|
||||
for i := ((k - 1) / nb) * nb; i >= 0; i -= nb {
|
||||
ib := min(nb, k-i)
|
||||
impl.Dlarft(lapack.Forward, lapack.RowWise, n-i, ib,
|
||||
a[i*lda+i:], lda,
|
||||
tau[i:],
|
||||
t, ldt)
|
||||
impl.Dlarfb(side, transt, lapack.Forward, lapack.RowWise, m, n-i, ib,
|
||||
a[i*lda+i:], lda,
|
||||
t, ldt,
|
||||
c[i:], ldc,
|
||||
work, ldwork)
|
||||
}
|
||||
return
|
||||
case !left && !notran:
|
||||
for i := 0; i < k; i += nb {
|
||||
ib := min(nb, k-i)
|
||||
impl.Dlarft(lapack.Forward, lapack.RowWise, n-i, ib,
|
||||
a[i*lda+i:], lda,
|
||||
tau[i:],
|
||||
t, ldt)
|
||||
impl.Dlarfb(side, transt, lapack.Forward, lapack.RowWise, m, n-i, ib,
|
||||
a[i*lda+i:], lda,
|
||||
t, ldt,
|
||||
c[i:], ldc,
|
||||
work, ldwork)
|
||||
}
|
||||
return
|
||||
}
|
||||
}
|
||||
+139
@@ -0,0 +1,139 @@
|
||||
// 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 native
|
||||
|
||||
import (
|
||||
"github.com/gonum/blas"
|
||||
"github.com/gonum/lapack"
|
||||
)
|
||||
|
||||
// Dormqr multiplies the matrix c by the othogonal matrix q defined by the
|
||||
// slices a and tau. A and tau are as returned from Dgeqrf.
|
||||
// C = Q * C if side == blas.Left and trans == blas.NoTrans
|
||||
// C = Q^T * C if side == blas.Left and trans == blas.Trans
|
||||
// C = C * Q if side == blas.Right and trans == blas.NoTrans
|
||||
// C = C * Q^T if side == blas.Right and trans == blas.Trans
|
||||
// If side == blas.Left, a is a matrix of side k×m, and if side == blas.Right
|
||||
// a is of size k×n. This uses a blocked algorithm.
|
||||
//
|
||||
// Work is temporary storage, and lwork specifies the usable memory length.
|
||||
// At minimum, lwork >= m if side == blas.Left and lwork >= n if side == blas.Right,
|
||||
// and this function will panic otherwise.
|
||||
// Dormqr uses a block algorithm, but the block size is limited
|
||||
// by the temporary space available. If lwork == -1, instead of performing Dormqr,
|
||||
// the optimal work length will be stored into work[0].
|
||||
//
|
||||
// Tau contains the householder scales and must have length at least k, and
|
||||
// this function will panic otherwise.
|
||||
func (impl Implementation) Dormqr(side blas.Side, trans blas.Transpose, m, n, k int, a []float64, lda int, tau, c []float64, ldc int, work []float64, lwork int) {
|
||||
left := side == blas.Left
|
||||
notran := trans == blas.NoTrans
|
||||
if left {
|
||||
checkMatrix(m, k, a, lda)
|
||||
} else {
|
||||
checkMatrix(n, k, a, lda)
|
||||
}
|
||||
checkMatrix(m, n, c, ldc)
|
||||
|
||||
const nbmax = 64
|
||||
nw := n
|
||||
if side == blas.Right {
|
||||
nw = m
|
||||
}
|
||||
opts := string(side) + string(trans)
|
||||
nb := min(nbmax, impl.Ilaenv(1, "DORMQR", opts, m, n, k, -1))
|
||||
lworkopt := max(1, nw) * nb
|
||||
if lwork == -1 {
|
||||
work[0] = float64(lworkopt)
|
||||
return
|
||||
}
|
||||
if left {
|
||||
if lwork < n {
|
||||
panic(badWork)
|
||||
}
|
||||
} else {
|
||||
if lwork < m {
|
||||
panic(badWork)
|
||||
}
|
||||
}
|
||||
if m == 0 || n == 0 || k == 0 {
|
||||
return
|
||||
}
|
||||
nbmin := 2
|
||||
|
||||
ldwork := nb
|
||||
if nb > 1 && nb < k {
|
||||
iws := nw * nb
|
||||
if lwork < iws {
|
||||
nb = lwork / nw
|
||||
nbmin = max(2, impl.Ilaenv(2, "DORMQR", opts, m, n, k, -1))
|
||||
}
|
||||
}
|
||||
if nb < nbmin || nb >= k {
|
||||
// Call unblocked code
|
||||
impl.Dorm2r(side, trans, m, n, k, a, lda, tau, c, ldc, work)
|
||||
return
|
||||
}
|
||||
ldt := nb
|
||||
t := make([]float64, nb*ldt)
|
||||
switch {
|
||||
case left && notran:
|
||||
for i := ((k - 1) / nb) * nb; i >= 0; i -= nb {
|
||||
ib := min(nb, k-i)
|
||||
impl.Dlarft(lapack.Forward, lapack.ColumnWise, m-i, ib,
|
||||
a[i*lda+i:], lda,
|
||||
tau[i:],
|
||||
t, ldt)
|
||||
impl.Dlarfb(side, trans, lapack.Forward, lapack.ColumnWise, m-i, n, ib,
|
||||
a[i*lda+i:], lda,
|
||||
t, ldt,
|
||||
c[i*ldc:], ldc,
|
||||
work, ldwork)
|
||||
}
|
||||
return
|
||||
case left && !notran:
|
||||
for i := 0; i < k; i += nb {
|
||||
ib := min(nb, k-i)
|
||||
impl.Dlarft(lapack.Forward, lapack.ColumnWise, m-i, ib,
|
||||
a[i*lda+i:], lda,
|
||||
tau[i:],
|
||||
t, ldt)
|
||||
impl.Dlarfb(side, trans, lapack.Forward, lapack.ColumnWise, m-i, n, ib,
|
||||
a[i*lda+i:], lda,
|
||||
t, ldt,
|
||||
c[i*ldc:], ldc,
|
||||
work, ldwork)
|
||||
}
|
||||
return
|
||||
case !left && notran:
|
||||
for i := 0; i < k; i += nb {
|
||||
ib := min(nb, k-i)
|
||||
impl.Dlarft(lapack.Forward, lapack.ColumnWise, n-i, ib,
|
||||
a[i*lda+i:], lda,
|
||||
tau[i:],
|
||||
t, ldt)
|
||||
impl.Dlarfb(side, trans, lapack.Forward, lapack.ColumnWise, m, n-i, ib,
|
||||
a[i*lda+i:], lda,
|
||||
t, ldt,
|
||||
c[i:], ldc,
|
||||
work, ldwork)
|
||||
}
|
||||
return
|
||||
case !left && !notran:
|
||||
for i := ((k - 1) / nb) * nb; i >= 0; i -= nb {
|
||||
ib := min(nb, k-i)
|
||||
impl.Dlarft(lapack.Forward, lapack.ColumnWise, n-i, ib,
|
||||
a[i*lda+i:], lda,
|
||||
tau[i:],
|
||||
t, ldt)
|
||||
impl.Dlarfb(side, trans, lapack.Forward, lapack.ColumnWise, m, n-i, ib,
|
||||
a[i*lda+i:], lda,
|
||||
t, ldt,
|
||||
c[i:], ldc,
|
||||
work, ldwork)
|
||||
}
|
||||
return
|
||||
}
|
||||
}
|
||||
+73
@@ -0,0 +1,73 @@
|
||||
// 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 native
|
||||
|
||||
import (
|
||||
"math"
|
||||
|
||||
"github.com/gonum/blas"
|
||||
"github.com/gonum/blas/blas64"
|
||||
)
|
||||
|
||||
// Dpotf2 computes the cholesky decomposition of the symmetric positive definite
|
||||
// matrix a. If ul == blas.Upper, then a is stored as an upper-triangular matrix,
|
||||
// and a = U^T U is stored in place into a. If ul == blas.Lower, then a = L L^T
|
||||
// is computed and stored in-place into a. If a is not positive definite, false
|
||||
// is returned. This is the unblocked version of the algorithm.
|
||||
func (Implementation) Dpotf2(ul blas.Uplo, n int, a []float64, lda int) (ok bool) {
|
||||
if ul != blas.Upper && ul != blas.Lower {
|
||||
panic(badUplo)
|
||||
}
|
||||
if n < 0 {
|
||||
panic(nLT0)
|
||||
}
|
||||
if lda < n {
|
||||
panic(badLdA)
|
||||
}
|
||||
if n == 0 {
|
||||
return true
|
||||
}
|
||||
bi := blas64.Implementation()
|
||||
if ul == blas.Upper {
|
||||
for j := 0; j < n; j++ {
|
||||
ajj := a[j*lda+j]
|
||||
if j != 0 {
|
||||
ajj -= bi.Ddot(j, a[j:], lda, a[j:], lda)
|
||||
}
|
||||
if ajj <= 0 || math.IsNaN(ajj) {
|
||||
a[j*lda+j] = ajj
|
||||
return false
|
||||
}
|
||||
ajj = math.Sqrt(ajj)
|
||||
a[j*lda+j] = ajj
|
||||
if j < n-1 {
|
||||
bi.Dgemv(blas.Trans, j, n-j-1,
|
||||
-1, a[j+1:], lda, a[j:], lda,
|
||||
1, a[j*lda+j+1:], 1)
|
||||
bi.Dscal(n-j-1, 1/ajj, a[j*lda+j+1:], 1)
|
||||
}
|
||||
}
|
||||
return true
|
||||
}
|
||||
for j := 0; j < n; j++ {
|
||||
ajj := a[j*lda+j]
|
||||
if j != 0 {
|
||||
ajj -= bi.Ddot(j, a[j*lda:], 1, a[j*lda:], 1)
|
||||
}
|
||||
if ajj <= 0 || math.IsNaN(ajj) {
|
||||
a[j*lda+j] = ajj
|
||||
return false
|
||||
}
|
||||
ajj = math.Sqrt(ajj)
|
||||
a[j*lda+j] = ajj
|
||||
if j < n-1 {
|
||||
bi.Dgemv(blas.NoTrans, n-j-1, j,
|
||||
-1, a[(j+1)*lda:], lda, a[j*lda:], 1,
|
||||
1, a[(j+1)*lda+j:], lda)
|
||||
bi.Dscal(n-j-1, 1/ajj, a[(j+1)*lda+j:], lda)
|
||||
}
|
||||
}
|
||||
return true
|
||||
}
|
||||
+75
@@ -0,0 +1,75 @@
|
||||
// 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 native
|
||||
|
||||
import (
|
||||
"github.com/gonum/blas"
|
||||
"github.com/gonum/blas/blas64"
|
||||
)
|
||||
|
||||
// Dpotrf computes the cholesky decomposition of the symmetric positive definite
|
||||
// matrix a. If ul == blas.Upper, then a is stored as an upper-triangular matrix,
|
||||
// and a = U U^T is stored in place into a. If ul == blas.Lower, then a = L L^T
|
||||
// is computed and stored in-place into a. If a is not positive definite, false
|
||||
// is returned. This is the blocked version of the algorithm.
|
||||
func (impl Implementation) Dpotrf(ul blas.Uplo, n int, a []float64, lda int) (ok bool) {
|
||||
bi := blas64.Implementation()
|
||||
if ul != blas.Upper && ul != blas.Lower {
|
||||
panic(badUplo)
|
||||
}
|
||||
if n < 0 {
|
||||
panic(nLT0)
|
||||
}
|
||||
if lda < n {
|
||||
panic(badLdA)
|
||||
}
|
||||
if n == 0 {
|
||||
return true
|
||||
}
|
||||
nb := impl.Ilaenv(1, "DPOTRF", string(ul), n, -1, -1, -1)
|
||||
if n <= nb {
|
||||
return impl.Dpotf2(ul, n, a, lda)
|
||||
}
|
||||
if ul == blas.Upper {
|
||||
for j := 0; j < n; j += nb {
|
||||
jb := min(nb, n-j)
|
||||
bi.Dsyrk(blas.Upper, blas.Trans, jb, j,
|
||||
-1, a[j:], lda,
|
||||
1, a[j*lda+j:], lda)
|
||||
ok = impl.Dpotf2(blas.Upper, jb, a[j*lda+j:], lda)
|
||||
if !ok {
|
||||
return ok
|
||||
}
|
||||
if j+jb < n {
|
||||
bi.Dgemm(blas.Trans, blas.NoTrans, jb, n-j-jb, j,
|
||||
-1, a[j:], lda, a[j+jb:], lda,
|
||||
1, a[j*lda+j+jb:], lda)
|
||||
bi.Dtrsm(blas.Left, blas.Upper, blas.Trans, blas.NonUnit, jb, n-j-jb,
|
||||
1, a[j*lda+j:], lda,
|
||||
a[j*lda+j+jb:], lda)
|
||||
}
|
||||
}
|
||||
return true
|
||||
}
|
||||
for j := 0; j < n; j += nb {
|
||||
jb := min(nb, n-j)
|
||||
bi.Dsyrk(blas.Lower, blas.NoTrans, jb, j,
|
||||
-1, a[j*lda:], lda,
|
||||
1, a[j*lda+j:], lda)
|
||||
ok := impl.Dpotf2(blas.Lower, jb, a[j*lda+j:], lda)
|
||||
if !ok {
|
||||
return ok
|
||||
}
|
||||
if j+jb < n {
|
||||
bi.Dgemm(blas.NoTrans, blas.Trans, n-j-jb, jb, j,
|
||||
-1, a[(j+jb)*lda:], lda, a[j*lda:], lda,
|
||||
1, a[(j+jb)*lda+j:], lda)
|
||||
bi.Dtrsm(blas.Right, blas.Lower, blas.Trans, blas.NonUnit, n-j-jb, jb,
|
||||
1, a[j*lda+j:], lda,
|
||||
a[(j+jb)*lda+j:], lda)
|
||||
}
|
||||
}
|
||||
return true
|
||||
}
|
||||
+31
@@ -0,0 +1,31 @@
|
||||
// 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 native
|
||||
|
||||
import (
|
||||
"github.com/gonum/blas"
|
||||
"github.com/gonum/blas/blas64"
|
||||
)
|
||||
|
||||
// Dtrtrs solves a triangular system of the form a * x = b or a^T * x = b. Dtrtrs
|
||||
// checks for singularity in a. If a is singular, false is returned and no solve
|
||||
// is performed. True is returned otherwise.
|
||||
func (impl Implementation) Dtrtrs(uplo blas.Uplo, trans blas.Transpose, diag blas.Diag, n, nrhs int, a []float64, lda int, b []float64, ldb int) (ok bool) {
|
||||
nounit := diag == blas.NonUnit
|
||||
if n == 0 {
|
||||
return false
|
||||
}
|
||||
// Check for singularity.
|
||||
if nounit {
|
||||
for i := 0; i < n; i++ {
|
||||
if a[i*lda+i] == 0 {
|
||||
return false
|
||||
}
|
||||
}
|
||||
}
|
||||
bi := blas64.Implementation()
|
||||
bi.Dtrsm(blas.Left, uplo, trans, diag, n, nrhs, 1, a, lda, b, ldb)
|
||||
return true
|
||||
}
|
||||
+92
@@ -0,0 +1,92 @@
|
||||
// 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 native
|
||||
|
||||
import (
|
||||
"math"
|
||||
|
||||
"github.com/gonum/lapack"
|
||||
)
|
||||
|
||||
// Implementation is the native Go implementation of LAPACK routines. It
|
||||
// is built on top of calls to the return of blas64.Implementation(), so while
|
||||
// this code is in pure Go, the underlying BLAS implementation may not be.
|
||||
type Implementation struct{}
|
||||
|
||||
var _ lapack.Float64 = Implementation{}
|
||||
|
||||
const (
|
||||
badDirect = "lapack: bad direct"
|
||||
badLdA = "lapack: index of a out of range"
|
||||
badSide = "lapack: bad side"
|
||||
badStore = "lapack: bad store"
|
||||
badTau = "lapack: tau has insufficient length"
|
||||
badTrans = "lapack: bad trans"
|
||||
badUplo = "lapack: illegal triangle"
|
||||
badWork = "lapack: insufficient working memory"
|
||||
badWorkStride = "lapack: insufficient working array stride"
|
||||
negDimension = "lapack: negative matrix dimension"
|
||||
nLT0 = "lapack: n < 0"
|
||||
shortWork = "lapack: working array shorter than declared"
|
||||
)
|
||||
|
||||
// checkMatrix verifies the parameters of a matrix input.
|
||||
func checkMatrix(m, n int, a []float64, lda int) {
|
||||
if m < 0 {
|
||||
panic("lapack: has negative number of rows")
|
||||
}
|
||||
if m < 0 {
|
||||
panic("lapack: has negative number of columns")
|
||||
}
|
||||
if lda < n {
|
||||
panic("lapack: stride less than number of columns")
|
||||
}
|
||||
if len(a) < (m-1)*lda+n {
|
||||
panic("lapack: insufficient matrix slice length")
|
||||
}
|
||||
}
|
||||
|
||||
func checkVector(n int, v []float64, inc int) {
|
||||
if n < 0 {
|
||||
panic("lapack: negative matrix length")
|
||||
}
|
||||
if (inc > 0 && (n-1)*inc >= len(v)) || (inc < 0 && (1-n)*inc >= len(v)) {
|
||||
panic("lapack: insufficient vector slice length")
|
||||
}
|
||||
}
|
||||
|
||||
func min(a, b int) int {
|
||||
if a < b {
|
||||
return a
|
||||
}
|
||||
return b
|
||||
}
|
||||
|
||||
func max(a, b int) int {
|
||||
if a > b {
|
||||
return a
|
||||
}
|
||||
return b
|
||||
}
|
||||
|
||||
// dlamch is a function in fortran, but since go forces IEEE-754 these are all
|
||||
// fixed values. Probably a way to get them as constants.
|
||||
// TODO(btracey): Is there a better way to find the smallest number such that 1+E > 1
|
||||
|
||||
var dlamchE, dlamchS, dlamchP float64
|
||||
|
||||
func init() {
|
||||
onePlusEps := math.Nextafter(1, math.Inf(1))
|
||||
eps := (math.Nextafter(1, math.Inf(1)) - 1) * 0.5
|
||||
dlamchE = eps
|
||||
sfmin := math.SmallestNonzeroFloat64
|
||||
small := 1 / math.MaxFloat64
|
||||
if small >= sfmin {
|
||||
sfmin = small * onePlusEps
|
||||
}
|
||||
dlamchS = sfmin
|
||||
radix := 2.0
|
||||
dlamchP = radix * eps
|
||||
}
|
||||
+31
@@ -0,0 +1,31 @@
|
||||
// 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 native
|
||||
|
||||
// Iladlc scans a matrix for its last non-zero column. Returns -1 if the matrix
|
||||
// is all zeros.
|
||||
func (Implementation) Iladlc(m, n int, a []float64, lda int) int {
|
||||
if n == 0 || m == 0 {
|
||||
return n - 1
|
||||
}
|
||||
checkMatrix(m, n, a, lda)
|
||||
|
||||
// Test common case where corner is non-zero.
|
||||
if a[n-1] != 0 || a[(m-1)*lda+(n-1)] != 0 {
|
||||
return n - 1
|
||||
}
|
||||
|
||||
// Scan each row tracking the highest column seen.
|
||||
highest := -1
|
||||
for i := 0; i < m; i++ {
|
||||
for j := n - 1; j >= 0; j-- {
|
||||
if a[i*lda+j] != 0 {
|
||||
highest = max(highest, j)
|
||||
break
|
||||
}
|
||||
}
|
||||
}
|
||||
return highest
|
||||
}
|
||||
+28
@@ -0,0 +1,28 @@
|
||||
// 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 native
|
||||
|
||||
// Iladlr scans a matrix for its last non-zero row. Returns -1 if the matrix
|
||||
// is all zeros.
|
||||
func (Implementation) Iladlr(m, n int, a []float64, lda int) int {
|
||||
if m == 0 {
|
||||
return m - 1
|
||||
}
|
||||
|
||||
checkMatrix(m, n, a, lda)
|
||||
|
||||
// Check the common case where the corner is non-zero
|
||||
if a[(m-1)*lda] != 0 || a[(m-1)*lda+n-1] != 0 {
|
||||
return m - 1
|
||||
}
|
||||
for i := m - 1; i >= 0; i-- {
|
||||
for j := 0; j < n; j++ {
|
||||
if a[i*lda+j] != 0 {
|
||||
return i
|
||||
}
|
||||
}
|
||||
}
|
||||
return -1
|
||||
}
|
||||
+375
@@ -0,0 +1,375 @@
|
||||
// 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 native
|
||||
|
||||
// Ilaenv returns algorithm tuning parameters for the algorithm given by the
|
||||
// input string. ispec specifies the parameter to return.
|
||||
// 1: The optimal block size
|
||||
// 2: The minimum block size for which the algorithm should be used.
|
||||
// 3: The crossover point below which an unblocked routine should be used.
|
||||
// 4: The number of shifts.
|
||||
// 5: The minumum column dimension for blocking to be used.
|
||||
// 6: The crossover point for SVD (to use QR factorization or not).
|
||||
// 7: The number of processors.
|
||||
// 8: The crossover point for multishift in QR and QZ methods for nonsymmetric eigenvalue problems.
|
||||
// 9: Maximum size of the subproblems in divide-and-conquer algorithms.
|
||||
// 10: ieee NaN arithmetic can be trusted not to trap.
|
||||
// 11: infinity arithmetic can be trusted not to trap.
|
||||
func (Implementation) Ilaenv(ispec int, s string, opts string, n1, n2, n3, n4 int) int {
|
||||
// TODO(btracey): Replace this with a constant lookup? A list of constants?
|
||||
// TODO: What is the difference between 2 and 3?
|
||||
sname := s[0] == 'S' || s[0] == 'D'
|
||||
cname := s[0] == 'C' || s[0] == 'Z'
|
||||
if !sname && !cname {
|
||||
panic("lapack: bad name")
|
||||
}
|
||||
c2 := s[1:3]
|
||||
c3 := s[3:6]
|
||||
c4 := c3[1:3]
|
||||
|
||||
switch ispec {
|
||||
default:
|
||||
panic("lapack: bad ispec")
|
||||
case 1:
|
||||
switch c2 {
|
||||
default:
|
||||
panic("lapack: bad function name")
|
||||
case "GE":
|
||||
switch c3 {
|
||||
default:
|
||||
panic("lapack: bad function name")
|
||||
case "TRF":
|
||||
if sname {
|
||||
return 64
|
||||
}
|
||||
return 64
|
||||
case "QRF", "RQF", "LQF", "QLF":
|
||||
if sname {
|
||||
return 32
|
||||
}
|
||||
return 32
|
||||
case "HRD":
|
||||
if sname {
|
||||
return 32
|
||||
}
|
||||
return 32
|
||||
case "BRD":
|
||||
if sname {
|
||||
return 32
|
||||
}
|
||||
return 32
|
||||
case "TRI":
|
||||
if sname {
|
||||
return 64
|
||||
}
|
||||
return 64
|
||||
}
|
||||
case "PO":
|
||||
switch c3 {
|
||||
default:
|
||||
panic("lapack: bad function name")
|
||||
case "TRF":
|
||||
if sname {
|
||||
return 64
|
||||
}
|
||||
return 64
|
||||
}
|
||||
case "SY":
|
||||
switch c3 {
|
||||
default:
|
||||
panic("lapack: bad function name")
|
||||
case "TRF":
|
||||
if sname {
|
||||
return 64
|
||||
}
|
||||
return 64
|
||||
case "TRD":
|
||||
return 32
|
||||
case "GST":
|
||||
return 64
|
||||
}
|
||||
case "HE":
|
||||
switch c3 {
|
||||
default:
|
||||
panic("lapack: bad function name")
|
||||
case "TRF":
|
||||
return 64
|
||||
case "TRD":
|
||||
return 32
|
||||
case "GST":
|
||||
return 64
|
||||
}
|
||||
case "OR":
|
||||
switch c3[0] {
|
||||
default:
|
||||
panic("lapack: bad function name")
|
||||
case 'G':
|
||||
switch c3[1:] {
|
||||
default:
|
||||
panic("lapack: bad function name")
|
||||
case "QR", "RQ", "LQ", "QL", "HR", "TR", "BR":
|
||||
return 32
|
||||
}
|
||||
case 'M':
|
||||
switch c3[1:] {
|
||||
default:
|
||||
panic("lapack: bad function name")
|
||||
case "QR", "RQ", "LQ", "QL", "HR", "TR", "BR":
|
||||
return 32
|
||||
}
|
||||
}
|
||||
case "UN":
|
||||
switch c3[0] {
|
||||
default:
|
||||
panic("lapack: bad function name")
|
||||
case 'G':
|
||||
switch c3[1:] {
|
||||
default:
|
||||
panic("lapack: bad function name")
|
||||
case "QR", "RQ", "LQ", "QL", "HR", "TR", "BR":
|
||||
return 32
|
||||
}
|
||||
case 'M':
|
||||
switch c3[1:] {
|
||||
default:
|
||||
panic("lapack: bad function name")
|
||||
case "QR", "RQ", "LQ", "QL", "HR", "TR", "BR":
|
||||
return 32
|
||||
}
|
||||
}
|
||||
case "GB":
|
||||
switch c3 {
|
||||
default:
|
||||
panic("lapack: bad function name")
|
||||
case "TRF":
|
||||
if sname {
|
||||
if n4 <= 64 {
|
||||
return 1
|
||||
}
|
||||
return 32
|
||||
}
|
||||
if n4 <= 64 {
|
||||
return 1
|
||||
}
|
||||
return 32
|
||||
}
|
||||
case "PB":
|
||||
switch c3 {
|
||||
default:
|
||||
panic("lapack: bad function name")
|
||||
case "TRF":
|
||||
if sname {
|
||||
if n4 <= 64 {
|
||||
return 1
|
||||
}
|
||||
return 32
|
||||
}
|
||||
if n4 <= 64 {
|
||||
return 1
|
||||
}
|
||||
return 32
|
||||
}
|
||||
case "TR":
|
||||
switch c3 {
|
||||
default:
|
||||
panic("lapack: bad function name")
|
||||
case "TRI":
|
||||
if sname {
|
||||
return 64
|
||||
}
|
||||
return 64
|
||||
}
|
||||
case "LA":
|
||||
switch c3 {
|
||||
default:
|
||||
panic("lapack: bad function name")
|
||||
case "UUM":
|
||||
if sname {
|
||||
return 64
|
||||
}
|
||||
return 64
|
||||
}
|
||||
case "ST":
|
||||
if sname && c3 == "EBZ" {
|
||||
return 1
|
||||
}
|
||||
panic("lapack: bad function name")
|
||||
}
|
||||
case 2:
|
||||
switch c2 {
|
||||
default:
|
||||
panic("lapack: bad function name")
|
||||
case "GE":
|
||||
switch c3 {
|
||||
default:
|
||||
panic("lapack: bad function name")
|
||||
case "QRF", "RQF", "LQF", "QLF":
|
||||
if sname {
|
||||
return 2
|
||||
}
|
||||
return 2
|
||||
case "HRD":
|
||||
if sname {
|
||||
return 2
|
||||
}
|
||||
return 2
|
||||
case "BRD":
|
||||
if sname {
|
||||
return 2
|
||||
}
|
||||
return 2
|
||||
case "TRI":
|
||||
if sname {
|
||||
return 2
|
||||
}
|
||||
return 2
|
||||
}
|
||||
case "SY":
|
||||
switch c3 {
|
||||
default:
|
||||
panic("lapack: bad function name")
|
||||
case "TRF":
|
||||
if sname {
|
||||
return 8
|
||||
}
|
||||
return 8
|
||||
case "TRD":
|
||||
if sname {
|
||||
return 2
|
||||
}
|
||||
panic("lapack: bad function name")
|
||||
}
|
||||
case "HE":
|
||||
if c3 == "TRD" {
|
||||
return 2
|
||||
}
|
||||
panic("lapack: bad function name")
|
||||
case "OR":
|
||||
if !sname {
|
||||
panic("lapack: bad function name")
|
||||
}
|
||||
switch c3[0] {
|
||||
default:
|
||||
panic("lapack: bad function name")
|
||||
case 'G':
|
||||
switch c4 {
|
||||
default:
|
||||
panic("lapack: bad function name")
|
||||
case "QR", "RQ", "LQ", "QL", "HR", "TR", "BR":
|
||||
return 2
|
||||
}
|
||||
case 'M':
|
||||
switch c4 {
|
||||
default:
|
||||
panic("lapack: bad function name")
|
||||
case "QR", "RQ", "LQ", "QL", "HR", "TR", "BR":
|
||||
return 2
|
||||
}
|
||||
}
|
||||
case "UN":
|
||||
switch c3[0] {
|
||||
default:
|
||||
panic("lapack: bad function name")
|
||||
case 'G':
|
||||
switch c4 {
|
||||
default:
|
||||
panic("lapack: bad function name")
|
||||
case "QR", "RQ", "LQ", "QL", "HR", "TR", "BR":
|
||||
return 2
|
||||
}
|
||||
case 'M':
|
||||
switch c4 {
|
||||
default:
|
||||
panic("lapack: bad function name")
|
||||
case "QR", "RQ", "LQ", "QL", "HR", "TR", "BR":
|
||||
return 2
|
||||
}
|
||||
}
|
||||
}
|
||||
case 3:
|
||||
switch c2 {
|
||||
default:
|
||||
panic("lapack: bad function name")
|
||||
case "GE":
|
||||
switch c3 {
|
||||
default:
|
||||
panic("lapack: bad function name")
|
||||
case "QRF", "RQF", "LQF", "QLF":
|
||||
if sname {
|
||||
return 128
|
||||
}
|
||||
return 128
|
||||
case "HRD":
|
||||
if sname {
|
||||
return 128
|
||||
}
|
||||
return 128
|
||||
case "BRD":
|
||||
if sname {
|
||||
return 128
|
||||
}
|
||||
return 128
|
||||
}
|
||||
case "SY":
|
||||
if sname && c3 == "TRD" {
|
||||
return 32
|
||||
}
|
||||
panic("lapack: bad function name")
|
||||
case "HE":
|
||||
if c3 == "TRD" {
|
||||
return 32
|
||||
}
|
||||
panic("lapack: bad function name")
|
||||
case "OR":
|
||||
switch c3[0] {
|
||||
default:
|
||||
panic("lapack: bad function name")
|
||||
case 'G':
|
||||
switch c4 {
|
||||
default:
|
||||
panic("lapack: bad function name")
|
||||
case "QR", "RQ", "LQ", "QL", "HR", "TR", "BR":
|
||||
return 128
|
||||
}
|
||||
}
|
||||
case "UN":
|
||||
switch c3[0] {
|
||||
default:
|
||||
panic("lapack: bad function name")
|
||||
case 'G':
|
||||
switch c4 {
|
||||
default:
|
||||
panic("lapack: bad function name")
|
||||
case "QR", "RQ", "LQ", "QL", "HR", "TR", "BR":
|
||||
return 128
|
||||
}
|
||||
}
|
||||
}
|
||||
case 4:
|
||||
// Used by xHSEQR
|
||||
return 6
|
||||
case 5:
|
||||
// Not used
|
||||
return 2
|
||||
case 6:
|
||||
// Used by xGELSS and xGESVD
|
||||
return min(n1, n2) * 1e6
|
||||
case 7:
|
||||
// Not used
|
||||
return 1
|
||||
case 8:
|
||||
// Used by xHSEQR
|
||||
return 50
|
||||
case 9:
|
||||
// used by xGELSD and xGESDD
|
||||
return 25
|
||||
case 10:
|
||||
// Go guarantees ieee
|
||||
return 1
|
||||
case 11:
|
||||
// Go guarantees ieee
|
||||
return 1
|
||||
}
|
||||
}
|
||||
Reference in New Issue
Block a user