+++ /dev/null
-// 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 testlapack
-
-import (
- "fmt"
- "math"
- "math/cmplx"
- "testing"
-
- "golang.org/x/exp/rand"
-
- "gonum.org/v1/gonum/blas"
- "gonum.org/v1/gonum/blas/blas64"
- "gonum.org/v1/gonum/floats"
- "gonum.org/v1/gonum/lapack"
-)
-
-const (
- // dlamchE is the machine epsilon. For IEEE this is 2^{-53}.
- dlamchE = 1.0 / (1 << 53)
- dlamchB = 2
- dlamchP = dlamchB * dlamchE
- // dlamchS is the smallest normal number. For IEEE this is 2^{-1022}.
- dlamchS = 1.0 / (1 << 256) / (1 << 256) / (1 << 256) / (1 << 254)
-)
-
-func max(a, b int) int {
- if a > b {
- return a
- }
- return b
-}
-
-func min(a, b int) int {
- if a < b {
- return a
- }
- return b
-}
-
-// worklen describes how much workspace a test should use.
-type worklen int
-
-const (
- minimumWork worklen = iota
- mediumWork
- optimumWork
-)
-
-// nanSlice allocates a new slice of length n filled with NaN.
-func nanSlice(n int) []float64 {
- s := make([]float64, n)
- for i := range s {
- s[i] = math.NaN()
- }
- return s
-}
-
-// randomSlice allocates a new slice of length n filled with random values.
-func randomSlice(n int, rnd *rand.Rand) []float64 {
- s := make([]float64, n)
- for i := range s {
- s[i] = rnd.NormFloat64()
- }
- return s
-}
-
-// nanGeneral allocates a new r×c general matrix filled with NaN values.
-func nanGeneral(r, c, stride int) blas64.General {
- if r < 0 || c < 0 {
- panic("bad matrix size")
- }
- if r == 0 || c == 0 {
- return blas64.General{Stride: max(1, stride)}
- }
- if stride < c {
- panic("bad stride")
- }
- return blas64.General{
- Rows: r,
- Cols: c,
- Stride: stride,
- Data: nanSlice((r-1)*stride + c),
- }
-}
-
-// randomGeneral allocates a new r×c general matrix filled with random
-// numbers. Out-of-range elements are filled with NaN values.
-func randomGeneral(r, c, stride int, rnd *rand.Rand) blas64.General {
- ans := nanGeneral(r, c, stride)
- for i := 0; i < r; i++ {
- for j := 0; j < c; j++ {
- ans.Data[i*ans.Stride+j] = rnd.NormFloat64()
- }
- }
- return ans
-}
-
-// randomHessenberg allocates a new n×n Hessenberg matrix filled with zeros
-// under the first subdiagonal and with random numbers elsewhere. Out-of-range
-// elements are filled with NaN values.
-func randomHessenberg(n, stride int, rnd *rand.Rand) blas64.General {
- ans := nanGeneral(n, n, stride)
- for i := 0; i < n; i++ {
- for j := 0; j < i-1; j++ {
- ans.Data[i*ans.Stride+j] = 0
- }
- for j := max(0, i-1); j < n; j++ {
- ans.Data[i*ans.Stride+j] = rnd.NormFloat64()
- }
- }
- return ans
-}
-
-// randomSchurCanonical returns a random, general matrix in Schur canonical
-// form, that is, block upper triangular with 1×1 and 2×2 diagonal blocks where
-// each 2×2 diagonal block has its diagonal elements equal and its off-diagonal
-// elements of opposite sign.
-func randomSchurCanonical(n, stride int, rnd *rand.Rand) blas64.General {
- t := randomGeneral(n, n, stride, rnd)
- // Zero out the lower triangle.
- for i := 0; i < t.Rows; i++ {
- for j := 0; j < i; j++ {
- t.Data[i*t.Stride+j] = 0
- }
- }
- // Randomly create 2×2 diagonal blocks.
- for i := 0; i < t.Rows; {
- if i == t.Rows-1 || rnd.Float64() < 0.5 {
- // 1×1 block.
- i++
- continue
- }
- // 2×2 block.
- // Diagonal elements equal.
- t.Data[(i+1)*t.Stride+i+1] = t.Data[i*t.Stride+i]
- // Off-diagonal elements of opposite sign.
- c := rnd.NormFloat64()
- if math.Signbit(c) == math.Signbit(t.Data[i*t.Stride+i+1]) {
- c *= -1
- }
- t.Data[(i+1)*t.Stride+i] = c
- i += 2
- }
- return t
-}
-
-// blockedUpperTriGeneral returns a normal random, general matrix in the form
-//
-// c-k-l k l
-// A = k [ 0 A12 A13 ] if r-k-l >= 0;
-// l [ 0 0 A23 ]
-// r-k-l [ 0 0 0 ]
-//
-// c-k-l k l
-// A = k [ 0 A12 A13 ] if r-k-l < 0;
-// r-k [ 0 0 A23 ]
-//
-// where the k×k matrix A12 and l×l matrix is non-singular
-// upper triangular. A23 is l×l upper triangular if r-k-l >= 0,
-// otherwise A23 is (r-k)×l upper trapezoidal.
-func blockedUpperTriGeneral(r, c, k, l, stride int, kblock bool, rnd *rand.Rand) blas64.General {
- t := l
- if kblock {
- t += k
- }
- ans := zeros(r, c, stride)
- for i := 0; i < min(r, t); i++ {
- var v float64
- for v == 0 {
- v = rnd.NormFloat64()
- }
- ans.Data[i*ans.Stride+i+(c-t)] = v
- }
- for i := 0; i < min(r, t); i++ {
- for j := i + (c - t) + 1; j < c; j++ {
- ans.Data[i*ans.Stride+j] = rnd.NormFloat64()
- }
- }
- return ans
-}
-
-// nanTriangular allocates a new r×c triangular matrix filled with NaN values.
-func nanTriangular(uplo blas.Uplo, n, stride int) blas64.Triangular {
- if n < 0 {
- panic("bad matrix size")
- }
- if n == 0 {
- return blas64.Triangular{
- Stride: max(1, stride),
- Uplo: uplo,
- Diag: blas.NonUnit,
- }
- }
- if stride < n {
- panic("bad stride")
- }
- return blas64.Triangular{
- N: n,
- Stride: stride,
- Data: nanSlice((n-1)*stride + n),
- Uplo: uplo,
- Diag: blas.NonUnit,
- }
-}
-
-// generalOutsideAllNaN returns whether all out-of-range elements have NaN
-// values.
-func generalOutsideAllNaN(a blas64.General) bool {
- // Check after last column.
- for i := 0; i < a.Rows-1; i++ {
- for _, v := range a.Data[i*a.Stride+a.Cols : i*a.Stride+a.Stride] {
- if !math.IsNaN(v) {
- return false
- }
- }
- }
- // Check after last element.
- last := (a.Rows-1)*a.Stride + a.Cols
- if a.Rows == 0 || a.Cols == 0 {
- last = 0
- }
- for _, v := range a.Data[last:] {
- if !math.IsNaN(v) {
- return false
- }
- }
- return true
-}
-
-// triangularOutsideAllNaN returns whether all out-of-triangle elements have NaN
-// values.
-func triangularOutsideAllNaN(a blas64.Triangular) bool {
- if a.Uplo == blas.Upper {
- // Check below diagonal.
- for i := 0; i < a.N; i++ {
- for _, v := range a.Data[i*a.Stride : i*a.Stride+i] {
- if !math.IsNaN(v) {
- return false
- }
- }
- }
- // Check after last column.
- for i := 0; i < a.N-1; i++ {
- for _, v := range a.Data[i*a.Stride+a.N : i*a.Stride+a.Stride] {
- if !math.IsNaN(v) {
- return false
- }
- }
- }
- } else {
- // Check above diagonal.
- for i := 0; i < a.N-1; i++ {
- for _, v := range a.Data[i*a.Stride+i+1 : i*a.Stride+a.Stride] {
- if !math.IsNaN(v) {
- return false
- }
- }
- }
- }
- // Check after last element.
- for _, v := range a.Data[max(0, a.N-1)*a.Stride+a.N:] {
- if !math.IsNaN(v) {
- return false
- }
- }
- return true
-}
-
-// transposeGeneral returns a new general matrix that is the transpose of the
-// input. Nothing is done with data outside the {rows, cols} limit of the general.
-func transposeGeneral(a blas64.General) blas64.General {
- ans := blas64.General{
- Rows: a.Cols,
- Cols: a.Rows,
- Stride: a.Rows,
- Data: make([]float64, a.Cols*a.Rows),
- }
- for i := 0; i < a.Rows; i++ {
- for j := 0; j < a.Cols; j++ {
- ans.Data[j*ans.Stride+i] = a.Data[i*a.Stride+j]
- }
- }
- return ans
-}
-
-// columnNorms returns the column norms of a.
-func columnNorms(m, n int, a []float64, lda int) []float64 {
- bi := blas64.Implementation()
- norms := make([]float64, n)
- for j := 0; j < n; j++ {
- norms[j] = bi.Dnrm2(m, a[j:], lda)
- }
- return norms
-}
-
-// extractVMat collects the single reflectors from a into a matrix.
-func extractVMat(m, n int, a []float64, lda int, direct lapack.Direct, store lapack.StoreV) blas64.General {
- k := min(m, n)
- switch {
- default:
- panic("not implemented")
- case direct == lapack.Forward && store == lapack.ColumnWise:
- v := blas64.General{
- Rows: m,
- Cols: k,
- Stride: k,
- Data: make([]float64, m*k),
- }
- for i := 0; i < k; i++ {
- for j := 0; j < i; j++ {
- v.Data[j*v.Stride+i] = 0
- }
- v.Data[i*v.Stride+i] = 1
- for j := i + 1; j < m; j++ {
- v.Data[j*v.Stride+i] = a[j*lda+i]
- }
- }
- return v
- case direct == lapack.Forward && store == lapack.RowWise:
- v := blas64.General{
- Rows: k,
- Cols: n,
- Stride: n,
- Data: make([]float64, k*n),
- }
- for i := 0; i < k; i++ {
- for j := 0; j < i; j++ {
- v.Data[i*v.Stride+j] = 0
- }
- v.Data[i*v.Stride+i] = 1
- for j := i + 1; j < n; j++ {
- v.Data[i*v.Stride+j] = a[i*lda+j]
- }
- }
- return v
- }
-}
-
-// constructBidiagonal constructs a bidiagonal matrix with the given diagonal
-// and off-diagonal elements.
-func constructBidiagonal(uplo blas.Uplo, n int, d, e []float64) blas64.General {
- bMat := blas64.General{
- Rows: n,
- Cols: n,
- Stride: n,
- Data: make([]float64, n*n),
- }
-
- for i := 0; i < n-1; i++ {
- bMat.Data[i*bMat.Stride+i] = d[i]
- if uplo == blas.Upper {
- bMat.Data[i*bMat.Stride+i+1] = e[i]
- } else {
- bMat.Data[(i+1)*bMat.Stride+i] = e[i]
- }
- }
- bMat.Data[(n-1)*bMat.Stride+n-1] = d[n-1]
- return bMat
-}
-
-// constructVMat transforms the v matrix based on the storage.
-func constructVMat(vMat blas64.General, store lapack.StoreV, direct lapack.Direct) blas64.General {
- m := vMat.Rows
- k := vMat.Cols
- switch {
- default:
- panic("not implemented")
- case store == lapack.ColumnWise && direct == lapack.Forward:
- ldv := k
- v := make([]float64, m*k)
- for i := 0; i < m; i++ {
- for j := 0; j < k; j++ {
- if j > i {
- v[i*ldv+j] = 0
- } else if j == i {
- v[i*ldv+i] = 1
- } else {
- v[i*ldv+j] = vMat.Data[i*vMat.Stride+j]
- }
- }
- }
- return blas64.General{
- Rows: m,
- Cols: k,
- Stride: k,
- Data: v,
- }
- case store == lapack.RowWise && direct == lapack.Forward:
- ldv := m
- v := make([]float64, m*k)
- for i := 0; i < m; i++ {
- for j := 0; j < k; j++ {
- if j > i {
- v[j*ldv+i] = 0
- } else if j == i {
- v[j*ldv+i] = 1
- } else {
- v[j*ldv+i] = vMat.Data[i*vMat.Stride+j]
- }
- }
- }
- return blas64.General{
- Rows: k,
- Cols: m,
- Stride: m,
- Data: v,
- }
- case store == lapack.ColumnWise && direct == lapack.Backward:
- rowsv := m
- ldv := k
- v := make([]float64, m*k)
- for i := 0; i < m; i++ {
- for j := 0; j < k; j++ {
- vrow := rowsv - i - 1
- vcol := k - j - 1
- if j > i {
- v[vrow*ldv+vcol] = 0
- } else if j == i {
- v[vrow*ldv+vcol] = 1
- } else {
- v[vrow*ldv+vcol] = vMat.Data[i*vMat.Stride+j]
- }
- }
- }
- return blas64.General{
- Rows: rowsv,
- Cols: ldv,
- Stride: ldv,
- Data: v,
- }
- case store == lapack.RowWise && direct == lapack.Backward:
- rowsv := k
- ldv := m
- v := make([]float64, m*k)
- for i := 0; i < m; i++ {
- for j := 0; j < k; j++ {
- vcol := ldv - i - 1
- vrow := k - j - 1
- if j > i {
- v[vrow*ldv+vcol] = 0
- } else if j == i {
- v[vrow*ldv+vcol] = 1
- } else {
- v[vrow*ldv+vcol] = vMat.Data[i*vMat.Stride+j]
- }
- }
- }
- return blas64.General{
- Rows: rowsv,
- Cols: ldv,
- Stride: ldv,
- Data: v,
- }
- }
-}
-
-func constructH(tau []float64, v blas64.General, store lapack.StoreV, direct lapack.Direct) blas64.General {
- m := v.Rows
- k := v.Cols
- if store == lapack.RowWise {
- m, k = k, m
- }
- h := blas64.General{
- Rows: m,
- Cols: m,
- Stride: m,
- Data: make([]float64, m*m),
- }
- for i := 0; i < m; i++ {
- h.Data[i*m+i] = 1
- }
- for i := 0; i < k; i++ {
- vecData := make([]float64, m)
- if store == lapack.ColumnWise {
- for j := 0; j < m; j++ {
- vecData[j] = v.Data[j*v.Cols+i]
- }
- } else {
- for j := 0; j < m; j++ {
- vecData[j] = v.Data[i*v.Cols+j]
- }
- }
- vec := blas64.Vector{
- Inc: 1,
- Data: vecData,
- }
-
- hi := blas64.General{
- Rows: m,
- Cols: m,
- Stride: m,
- Data: make([]float64, m*m),
- }
- for i := 0; i < m; i++ {
- hi.Data[i*m+i] = 1
- }
- // hi = I - tau * v * v^T
- blas64.Ger(-tau[i], vec, vec, hi)
-
- hcopy := blas64.General{
- Rows: m,
- Cols: m,
- Stride: m,
- Data: make([]float64, m*m),
- }
- copy(hcopy.Data, h.Data)
- if direct == lapack.Forward {
- // H = H * H_I in forward mode
- blas64.Gemm(blas.NoTrans, blas.NoTrans, 1, hcopy, hi, 0, h)
- } else {
- // H = H_I * H in backward mode
- blas64.Gemm(blas.NoTrans, blas.NoTrans, 1, hi, hcopy, 0, h)
- }
- }
- return h
-}
-
-// constructQ constructs the Q matrix from the result of dgeqrf and dgeqr2.
-func constructQ(kind string, m, n int, a []float64, lda int, tau []float64) blas64.General {
- k := min(m, n)
- return constructQK(kind, m, n, k, a, lda, tau)
-}
-
-// constructQK constructs the Q matrix from the result of dgeqrf and dgeqr2 using
-// the first k reflectors.
-func constructQK(kind string, m, n, k int, a []float64, lda int, tau []float64) blas64.General {
- var sz int
- switch kind {
- case "QR":
- sz = m
- case "LQ", "RQ":
- sz = n
- }
-
- q := blas64.General{
- Rows: sz,
- Cols: sz,
- Stride: sz,
- Data: make([]float64, sz*sz),
- }
- for i := 0; i < sz; i++ {
- q.Data[i*sz+i] = 1
- }
- qCopy := blas64.General{
- Rows: q.Rows,
- Cols: q.Cols,
- Stride: q.Stride,
- Data: make([]float64, len(q.Data)),
- }
- for i := 0; i < k; i++ {
- h := blas64.General{
- Rows: sz,
- Cols: sz,
- Stride: sz,
- Data: make([]float64, sz*sz),
- }
- for j := 0; j < sz; j++ {
- h.Data[j*sz+j] = 1
- }
- vVec := blas64.Vector{
- Inc: 1,
- Data: make([]float64, sz),
- }
- switch kind {
- case "QR":
- vVec.Data[i] = 1
- for j := i + 1; j < sz; j++ {
- vVec.Data[j] = a[lda*j+i]
- }
- case "LQ":
- vVec.Data[i] = 1
- for j := i + 1; j < sz; j++ {
- vVec.Data[j] = a[i*lda+j]
- }
- case "RQ":
- for j := 0; j < n-k+i; j++ {
- vVec.Data[j] = a[(m-k+i)*lda+j]
- }
- vVec.Data[n-k+i] = 1
- }
- blas64.Ger(-tau[i], vVec, vVec, h)
- copy(qCopy.Data, q.Data)
- // Multiply q by the new h.
- switch kind {
- case "QR", "RQ":
- blas64.Gemm(blas.NoTrans, blas.NoTrans, 1, qCopy, h, 0, q)
- case "LQ":
- blas64.Gemm(blas.NoTrans, blas.NoTrans, 1, h, qCopy, 0, q)
- }
- }
- return q
-}
-
-// checkBidiagonal checks the bidiagonal decomposition from dlabrd and dgebd2.
-// The input to this function is the answer returned from the routines, stored
-// in a, d, e, tauP, and tauQ. The data of original A matrix (before
-// decomposition) is input in aCopy.
-//
-// checkBidiagonal constructs the V and U matrices, and from them constructs Q
-// and P. Using these constructions, it checks that Q^T * A * P and checks that
-// the result is bidiagonal.
-func checkBidiagonal(t *testing.T, m, n, nb int, a []float64, lda int, d, e, tauP, tauQ, aCopy []float64) {
- // Check the answer.
- // Construct V and U.
- qMat := constructQPBidiagonal(lapack.ApplyQ, m, n, nb, a, lda, tauQ)
- pMat := constructQPBidiagonal(lapack.ApplyP, m, n, nb, a, lda, tauP)
-
- // Compute Q^T * A * P.
- aMat := blas64.General{
- Rows: m,
- Cols: n,
- Stride: lda,
- Data: make([]float64, len(aCopy)),
- }
- copy(aMat.Data, aCopy)
-
- tmp1 := blas64.General{
- Rows: m,
- Cols: n,
- Stride: n,
- Data: make([]float64, m*n),
- }
- blas64.Gemm(blas.Trans, blas.NoTrans, 1, qMat, aMat, 0, tmp1)
- tmp2 := blas64.General{
- Rows: m,
- Cols: n,
- Stride: n,
- Data: make([]float64, m*n),
- }
- blas64.Gemm(blas.NoTrans, blas.NoTrans, 1, tmp1, pMat, 0, tmp2)
-
- // Check that the first nb rows and cols of tm2 are upper bidiagonal
- // if m >= n, and lower bidiagonal otherwise.
- correctDiag := true
- matchD := true
- matchE := true
- for i := 0; i < m; i++ {
- for j := 0; j < n; j++ {
- if i >= nb && j >= nb {
- continue
- }
- v := tmp2.Data[i*tmp2.Stride+j]
- if i == j {
- if math.Abs(d[i]-v) > 1e-12 {
- matchD = false
- }
- continue
- }
- if m >= n && i == j-1 {
- if math.Abs(e[j-1]-v) > 1e-12 {
- matchE = false
- }
- continue
- }
- if m < n && i-1 == j {
- if math.Abs(e[i-1]-v) > 1e-12 {
- matchE = false
- }
- continue
- }
- if math.Abs(v) > 1e-12 {
- correctDiag = false
- }
- }
- }
- if !correctDiag {
- t.Errorf("Updated A not bi-diagonal")
- }
- if !matchD {
- fmt.Println("d = ", d)
- t.Errorf("D Mismatch")
- }
- if !matchE {
- t.Errorf("E mismatch")
- }
-}
-
-// constructQPBidiagonal constructs Q or P from the Bidiagonal decomposition
-// computed by dlabrd and bgebd2.
-func constructQPBidiagonal(vect lapack.DecompUpdate, m, n, nb int, a []float64, lda int, tau []float64) blas64.General {
- sz := n
- if vect == lapack.ApplyQ {
- sz = m
- }
-
- var ldv int
- var v blas64.General
- if vect == lapack.ApplyQ {
- ldv = nb
- v = blas64.General{
- Rows: m,
- Cols: nb,
- Stride: ldv,
- Data: make([]float64, m*ldv),
- }
- } else {
- ldv = n
- v = blas64.General{
- Rows: nb,
- Cols: n,
- Stride: ldv,
- Data: make([]float64, m*ldv),
- }
- }
-
- if vect == lapack.ApplyQ {
- if m >= n {
- for i := 0; i < m; i++ {
- for j := 0; j <= min(nb-1, i); j++ {
- if i == j {
- v.Data[i*ldv+j] = 1
- continue
- }
- v.Data[i*ldv+j] = a[i*lda+j]
- }
- }
- } else {
- for i := 1; i < m; i++ {
- for j := 0; j <= min(nb-1, i-1); j++ {
- if i-1 == j {
- v.Data[i*ldv+j] = 1
- continue
- }
- v.Data[i*ldv+j] = a[i*lda+j]
- }
- }
- }
- } else {
- if m < n {
- for i := 0; i < nb; i++ {
- for j := i; j < n; j++ {
- if i == j {
- v.Data[i*ldv+j] = 1
- continue
- }
- v.Data[i*ldv+j] = a[i*lda+j]
- }
- }
- } else {
- for i := 0; i < nb; i++ {
- for j := i + 1; j < n; j++ {
- if j-1 == i {
- v.Data[i*ldv+j] = 1
- continue
- }
- v.Data[i*ldv+j] = a[i*lda+j]
- }
- }
- }
- }
-
- // The variable name is a computation of Q, but the algorithm is mostly the
- // same for computing P (just with different data).
- qMat := blas64.General{
- Rows: sz,
- Cols: sz,
- Stride: sz,
- Data: make([]float64, sz*sz),
- }
- hMat := blas64.General{
- Rows: sz,
- Cols: sz,
- Stride: sz,
- Data: make([]float64, sz*sz),
- }
- // set Q to I
- for i := 0; i < sz; i++ {
- qMat.Data[i*qMat.Stride+i] = 1
- }
- for i := 0; i < nb; i++ {
- qCopy := blas64.General{Rows: qMat.Rows, Cols: qMat.Cols, Stride: qMat.Stride, Data: make([]float64, len(qMat.Data))}
- copy(qCopy.Data, qMat.Data)
-
- // Set g and h to I
- for i := 0; i < sz; i++ {
- for j := 0; j < sz; j++ {
- if i == j {
- hMat.Data[i*sz+j] = 1
- } else {
- hMat.Data[i*sz+j] = 0
- }
- }
- }
- var vi blas64.Vector
- // H -= tauQ[i] * v[i] * v[i]^t
- if vect == lapack.ApplyQ {
- vi = blas64.Vector{
- Inc: v.Stride,
- Data: v.Data[i:],
- }
- } else {
- vi = blas64.Vector{
- Inc: 1,
- Data: v.Data[i*v.Stride:],
- }
- }
- blas64.Ger(-tau[i], vi, vi, hMat)
- // Q = Q * G[1]
- blas64.Gemm(blas.NoTrans, blas.NoTrans, 1, qCopy, hMat, 0, qMat)
- }
- return qMat
-}
-
-// printRowise prints the matrix with one row per line. This is useful for debugging.
-// If beyond is true, it prints beyond the final column to lda. If false, only
-// the columns are printed.
-func printRowise(a []float64, m, n, lda int, beyond bool) {
- for i := 0; i < m; i++ {
- end := n
- if beyond {
- end = lda
- }
- fmt.Println(a[i*lda : i*lda+end])
- }
-}
-
-// isOrthonormal checks that a general matrix is orthonormal.
-func isOrthonormal(q blas64.General) bool {
- n := q.Rows
- for i := 0; i < n; i++ {
- for j := i; j < n; j++ {
- dot := blas64.Dot(n,
- blas64.Vector{Inc: 1, Data: q.Data[i*q.Stride:]},
- blas64.Vector{Inc: 1, Data: q.Data[j*q.Stride:]},
- )
- if math.IsNaN(dot) {
- return false
- }
- if i == j {
- if math.Abs(dot-1) > 1e-10 {
- return false
- }
- } else {
- if math.Abs(dot) > 1e-10 {
- return false
- }
- }
- }
- }
- return true
-}
-
-// copyMatrix copies an m×n matrix src of stride n into an m×n matrix dst of stride ld.
-func copyMatrix(m, n int, dst []float64, ld int, src []float64) {
- for i := 0; i < m; i++ {
- copy(dst[i*ld:i*ld+n], src[i*n:i*n+n])
- }
-}
-
-func copyGeneral(dst, src blas64.General) {
- r := min(dst.Rows, src.Rows)
- c := min(dst.Cols, src.Cols)
- for i := 0; i < r; i++ {
- copy(dst.Data[i*dst.Stride:i*dst.Stride+c], src.Data[i*src.Stride:i*src.Stride+c])
- }
-}
-
-// cloneGeneral allocates and returns an exact copy of the given general matrix.
-func cloneGeneral(a blas64.General) blas64.General {
- c := a
- c.Data = make([]float64, len(a.Data))
- copy(c.Data, a.Data)
- return c
-}
-
-// equalApprox returns whether the matrices A and B of order n are approximately
-// equal within given tolerance.
-func equalApprox(m, n int, a []float64, lda int, b []float64, tol float64) bool {
- for i := 0; i < m; i++ {
- for j := 0; j < n; j++ {
- diff := a[i*lda+j] - b[i*n+j]
- if math.IsNaN(diff) || math.Abs(diff) > tol {
- return false
- }
- }
- }
- return true
-}
-
-// equalApproxGeneral returns whether the general matrices a and b are
-// approximately equal within given tolerance.
-func equalApproxGeneral(a, b blas64.General, tol float64) bool {
- if a.Rows != b.Rows || a.Cols != b.Cols {
- panic("bad input")
- }
- for i := 0; i < a.Rows; i++ {
- for j := 0; j < a.Cols; j++ {
- diff := a.Data[i*a.Stride+j] - b.Data[i*b.Stride+j]
- if math.IsNaN(diff) || math.Abs(diff) > tol {
- return false
- }
- }
- }
- return true
-}
-
-// equalApproxTriangular returns whether the triangular matrices A and B of
-// order n are approximately equal within given tolerance.
-func equalApproxTriangular(upper bool, n int, a []float64, lda int, b []float64, tol float64) bool {
- if upper {
- for i := 0; i < n; i++ {
- for j := i; j < n; j++ {
- diff := a[i*lda+j] - b[i*n+j]
- if math.IsNaN(diff) || math.Abs(diff) > tol {
- return false
- }
- }
- }
- return true
- }
- for i := 0; i < n; i++ {
- for j := 0; j <= i; j++ {
- diff := a[i*lda+j] - b[i*n+j]
- if math.IsNaN(diff) || math.Abs(diff) > tol {
- return false
- }
- }
- }
- return true
-}
-
-func equalApproxSymmetric(a, b blas64.Symmetric, tol float64) bool {
- if a.Uplo != b.Uplo {
- return false
- }
- if a.N != b.N {
- return false
- }
- if a.Uplo == blas.Upper {
- for i := 0; i < a.N; i++ {
- for j := i; j < a.N; j++ {
- if !floats.EqualWithinAbsOrRel(a.Data[i*a.Stride+j], b.Data[i*b.Stride+j], tol, tol) {
- return false
- }
- }
- }
- return true
- }
- for i := 0; i < a.N; i++ {
- for j := 0; j <= i; j++ {
- if !floats.EqualWithinAbsOrRel(a.Data[i*a.Stride+j], b.Data[i*b.Stride+j], tol, tol) {
- return false
- }
- }
- }
- return true
-}
-
-// randSymBand creates a random symmetric banded matrix, and returns both the
-// random matrix and the equivalent Symmetric matrix for testing. rnder
-// specifies the random number
-func randSymBand(ul blas.Uplo, n, ldab, kb int, rnd *rand.Rand) (blas64.Symmetric, blas64.SymmetricBand) {
- // A matrix is positive definite if and only if it has a Cholesky
- // decomposition. Generate a random banded lower triangular matrix
- // to construct the random symmetric matrix.
- a := make([]float64, n*n)
- for i := 0; i < n; i++ {
- for j := max(0, i-kb); j <= i; j++ {
- a[i*n+j] = rnd.NormFloat64()
- }
- a[i*n+i] = math.Abs(a[i*n+i])
- // Add an extra amound to the diagonal in order to improve the condition number.
- a[i*n+i] += 1.5 * rnd.Float64()
- }
- agen := blas64.General{
- Rows: n,
- Cols: n,
- Stride: n,
- Data: a,
- }
-
- // Construct the SymDense from a*a^T
- c := make([]float64, n*n)
- cgen := blas64.General{
- Rows: n,
- Cols: n,
- Stride: n,
- Data: c,
- }
- blas64.Gemm(blas.NoTrans, blas.Trans, 1, agen, agen, 0, cgen)
- sym := blas64.Symmetric{
- N: n,
- Stride: n,
- Data: c,
- Uplo: ul,
- }
-
- b := symToSymBand(ul, c, n, n, kb, ldab)
- band := blas64.SymmetricBand{
- N: n,
- K: kb,
- Stride: ldab,
- Data: b,
- Uplo: ul,
- }
-
- return sym, band
-}
-
-// symToSymBand takes the data in a Symmetric matrix and returns a
-// SymmetricBanded matrix.
-func symToSymBand(ul blas.Uplo, a []float64, n, lda, kb, ldab int) []float64 {
- if ul == blas.Upper {
- band := make([]float64, (n-1)*ldab+kb+1)
- for i := 0; i < n; i++ {
- for j := i; j < min(i+kb+1, n); j++ {
- band[i*ldab+j-i] = a[i*lda+j]
- }
- }
- return band
- }
- band := make([]float64, (n-1)*ldab+kb+1)
- for i := 0; i < n; i++ {
- for j := max(0, i-kb); j <= i; j++ {
- band[i*ldab+j-i+kb] = a[i*lda+j]
- }
- }
- return band
-}
-
-// symBandToSym takes a banded symmetric matrix and returns the same data as
-// a Symmetric matrix.
-func symBandToSym(ul blas.Uplo, band []float64, n, kb, ldab int) blas64.Symmetric {
- sym := make([]float64, n*n)
- if ul == blas.Upper {
- for i := 0; i < n; i++ {
- for j := 0; j < min(kb+1+i, n)-i; j++ {
- sym[i*n+i+j] = band[i*ldab+j]
- }
- }
- } else {
- for i := 0; i < n; i++ {
- for j := kb - min(i, kb); j < kb+1; j++ {
- sym[i*n+i-kb+j] = band[i*ldab+j]
- }
- }
- }
- return blas64.Symmetric{
- N: n,
- Stride: n,
- Data: sym,
- Uplo: ul,
- }
-}
-
-// eye returns an identity matrix of given order and stride.
-func eye(n, stride int) blas64.General {
- ans := nanGeneral(n, n, stride)
- for i := 0; i < n; i++ {
- for j := 0; j < n; j++ {
- ans.Data[i*ans.Stride+j] = 0
- }
- ans.Data[i*ans.Stride+i] = 1
- }
- return ans
-}
-
-// zeros returns an m×n matrix with given stride filled with zeros.
-func zeros(m, n, stride int) blas64.General {
- a := nanGeneral(m, n, stride)
- for i := 0; i < m; i++ {
- for j := 0; j < n; j++ {
- a.Data[i*a.Stride+j] = 0
- }
- }
- return a
-}
-
-// extract2x2Block returns the elements of T at [0,0], [0,1], [1,0], and [1,1].
-func extract2x2Block(t []float64, ldt int) (a, b, c, d float64) {
- return t[0], t[1], t[ldt], t[ldt+1]
-}
-
-// isSchurCanonical returns whether the 2×2 matrix [a b; c d] is in Schur
-// canonical form.
-func isSchurCanonical(a, b, c, d float64) bool {
- return c == 0 || (a == d && math.Signbit(b) != math.Signbit(c))
-}
-
-// isSchurCanonicalGeneral returns whether T is block upper triangular with 1×1
-// and 2×2 diagonal blocks, each 2×2 block in Schur canonical form. The function
-// checks only along the diagonal and the first subdiagonal, otherwise the lower
-// triangle is not accessed.
-func isSchurCanonicalGeneral(t blas64.General) bool {
- if t.Rows != t.Cols {
- panic("invalid matrix")
- }
- for i := 0; i < t.Rows-1; {
- if t.Data[(i+1)*t.Stride+i] == 0 {
- // 1×1 block.
- i++
- continue
- }
- // 2×2 block.
- a, b, c, d := extract2x2Block(t.Data[i*t.Stride+i:], t.Stride)
- if !isSchurCanonical(a, b, c, d) {
- return false
- }
- i += 2
- }
- return true
-}
-
-// schurBlockEigenvalues returns the two eigenvalues of the 2×2 matrix [a b; c d]
-// that must be in Schur canonical form.
-func schurBlockEigenvalues(a, b, c, d float64) (ev1, ev2 complex128) {
- if !isSchurCanonical(a, b, c, d) {
- panic("block not in Schur canonical form")
- }
- if c == 0 {
- return complex(a, 0), complex(d, 0)
- }
- im := math.Sqrt(-b * c)
- return complex(a, im), complex(a, -im)
-}
-
-// schurBlockSize returns the size of the diagonal block at i-th row in the
-// upper quasi-triangular matrix t in Schur canonical form, and whether i points
-// to the first row of the block. For zero-sized matrices the function returns 0
-// and true.
-func schurBlockSize(t blas64.General, i int) (size int, first bool) {
- if t.Rows != t.Cols {
- panic("matrix not square")
- }
- if t.Rows == 0 {
- return 0, true
- }
- if i < 0 || t.Rows <= i {
- panic("index out of range")
- }
-
- first = true
- if i > 0 && t.Data[i*t.Stride+i-1] != 0 {
- // There is a non-zero element to the left, therefore i must
- // point to the second row in a 2×2 diagonal block.
- first = false
- i--
- }
- size = 1
- if i+1 < t.Rows && t.Data[(i+1)*t.Stride+i] != 0 {
- // There is a non-zero element below, this must be a 2×2
- // diagonal block.
- size = 2
- }
- return size, first
-}
-
-// containsComplex returns whether z is approximately equal to one of the complex
-// numbers in v. If z is found, its index in v will be also returned.
-func containsComplex(v []complex128, z complex128, tol float64) (found bool, index int) {
- for i := range v {
- if cmplx.Abs(v[i]-z) < tol {
- return true, i
- }
- }
- return false, -1
-}
-
-// isAllNaN returns whether x contains only NaN values.
-func isAllNaN(x []float64) bool {
- for _, v := range x {
- if !math.IsNaN(v) {
- return false
- }
- }
- return true
-}
-
-// isUpperHessenberg returns whether h contains only zeros below the
-// subdiagonal.
-func isUpperHessenberg(h blas64.General) bool {
- if h.Rows != h.Cols {
- panic("matrix not square")
- }
- n := h.Rows
- for i := 0; i < n; i++ {
- for j := 0; j < n; j++ {
- if i > j+1 && h.Data[i*h.Stride+j] != 0 {
- return false
- }
- }
- }
- return true
-}
-
-// isUpperTriangular returns whether a contains only zeros below the diagonal.
-func isUpperTriangular(a blas64.General) bool {
- n := a.Rows
- for i := 1; i < n; i++ {
- for j := 0; j < i; j++ {
- if a.Data[i*a.Stride+j] != 0 {
- return false
- }
- }
- }
- return true
-}
-
-// unbalancedSparseGeneral returns an m×n dense matrix with a random sparse
-// structure consisting of nz nonzero elements. The matrix will be unbalanced by
-// multiplying each element randomly by its row or column index.
-func unbalancedSparseGeneral(m, n, stride int, nonzeros int, rnd *rand.Rand) blas64.General {
- a := zeros(m, n, stride)
- for k := 0; k < nonzeros; k++ {
- i := rnd.Intn(n)
- j := rnd.Intn(n)
- if rnd.Float64() < 0.5 {
- a.Data[i*stride+j] = float64(i+1) * rnd.NormFloat64()
- } else {
- a.Data[i*stride+j] = float64(j+1) * rnd.NormFloat64()
- }
- }
- return a
-}
-
-// columnOf returns a copy of the j-th column of a.
-func columnOf(a blas64.General, j int) []float64 {
- if j < 0 || a.Cols <= j {
- panic("bad column index")
- }
- col := make([]float64, a.Rows)
- for i := range col {
- col[i] = a.Data[i*a.Stride+j]
- }
- return col
-}
-
-// isRightEigenvectorOf returns whether the vector xRe+i*xIm, where i is the
-// imaginary unit, is the right eigenvector of A corresponding to the eigenvalue
-// lambda.
-//
-// A right eigenvector corresponding to a complex eigenvalue λ is a complex
-// non-zero vector x such that
-// A x = λ x.
-func isRightEigenvectorOf(a blas64.General, xRe, xIm []float64, lambda complex128, tol float64) bool {
- if a.Rows != a.Cols {
- panic("matrix not square")
- }
-
- if imag(lambda) != 0 && xIm == nil {
- // Complex eigenvalue of a real matrix cannot have a real
- // eigenvector.
- return false
- }
-
- n := a.Rows
-
- // Compute A real(x) and store the result into xReAns.
- xReAns := make([]float64, n)
- blas64.Gemv(blas.NoTrans, 1, a, blas64.Vector{1, xRe}, 0, blas64.Vector{1, xReAns})
-
- if imag(lambda) == 0 && xIm == nil {
- // Real eigenvalue and eigenvector.
-
- // Compute λx and store the result into lambdax.
- lambdax := make([]float64, n)
- floats.AddScaled(lambdax, real(lambda), xRe)
-
- // This is expressed as the inverse to catch the case
- // xReAns_i = Inf and lambdax_i = Inf of the same sign.
- return !(floats.Distance(xReAns, lambdax, math.Inf(1)) > tol)
- }
-
- // Complex eigenvector, and real or complex eigenvalue.
-
- // Compute A imag(x) and store the result into xImAns.
- xImAns := make([]float64, n)
- blas64.Gemv(blas.NoTrans, 1, a, blas64.Vector{1, xIm}, 0, blas64.Vector{1, xImAns})
-
- // Compute λx and store the result into lambdax.
- lambdax := make([]complex128, n)
- for i := range lambdax {
- lambdax[i] = lambda * complex(xRe[i], xIm[i])
- }
-
- for i, v := range lambdax {
- ax := complex(xReAns[i], xImAns[i])
- if cmplx.Abs(v-ax) > tol {
- return false
- }
- }
- return true
-}
-
-// isLeftEigenvectorOf returns whether the vector yRe+i*yIm, where i is the
-// imaginary unit, is the left eigenvector of A corresponding to the eigenvalue
-// lambda.
-//
-// A left eigenvector corresponding to a complex eigenvalue λ is a complex
-// non-zero vector y such that
-// y^H A = λ y^H,
-// which is equivalent for real A to
-// A^T y = conj(λ) y,
-func isLeftEigenvectorOf(a blas64.General, yRe, yIm []float64, lambda complex128, tol float64) bool {
- if a.Rows != a.Cols {
- panic("matrix not square")
- }
-
- if imag(lambda) != 0 && yIm == nil {
- // Complex eigenvalue of a real matrix cannot have a real
- // eigenvector.
- return false
- }
-
- n := a.Rows
-
- // Compute A^T real(y) and store the result into yReAns.
- yReAns := make([]float64, n)
- blas64.Gemv(blas.Trans, 1, a, blas64.Vector{1, yRe}, 0, blas64.Vector{1, yReAns})
-
- if imag(lambda) == 0 && yIm == nil {
- // Real eigenvalue and eigenvector.
-
- // Compute λy and store the result into lambday.
- lambday := make([]float64, n)
- floats.AddScaled(lambday, real(lambda), yRe)
-
- // This is expressed as the inverse to catch the case
- // yReAns_i = Inf and lambday_i = Inf of the same sign.
- return !(floats.Distance(yReAns, lambday, math.Inf(1)) > tol)
- }
-
- // Complex eigenvector, and real or complex eigenvalue.
-
- // Compute A^T imag(y) and store the result into yImAns.
- yImAns := make([]float64, n)
- blas64.Gemv(blas.Trans, 1, a, blas64.Vector{1, yIm}, 0, blas64.Vector{1, yImAns})
-
- // Compute conj(λ)y and store the result into lambday.
- lambda = cmplx.Conj(lambda)
- lambday := make([]complex128, n)
- for i := range lambday {
- lambday[i] = lambda * complex(yRe[i], yIm[i])
- }
-
- for i, v := range lambday {
- ay := complex(yReAns[i], yImAns[i])
- if cmplx.Abs(v-ay) > tol {
- return false
- }
- }
- return true
-}
-
-// rootsOfUnity returns the n complex numbers whose n-th power is equal to 1.
-func rootsOfUnity(n int) []complex128 {
- w := make([]complex128, n)
- for i := 0; i < n; i++ {
- angle := math.Pi * float64(2*i) / float64(n)
- w[i] = complex(math.Cos(angle), math.Sin(angle))
- }
- return w
-}
-
-// randomOrthogonal returns an n×n random orthogonal matrix.
-func randomOrthogonal(n int, rnd *rand.Rand) blas64.General {
- q := eye(n, n)
- x := make([]float64, n)
- v := make([]float64, n)
- for j := 0; j < n-1; j++ {
- // x represents the j-th column of a random matrix.
- for i := 0; i < j; i++ {
- x[i] = 0
- }
- for i := j; i < n; i++ {
- x[i] = rnd.NormFloat64()
- }
- // Compute v that represents the elementary reflector that
- // annihilates the subdiagonal elements of x.
- reflector(v, x, j)
- // Compute Q * H_j and store the result into Q.
- applyReflector(q, q, v)
- }
- if !isOrthonormal(q) {
- panic("Q not orthogonal")
- }
- return q
-}
-
-// reflector generates a Householder reflector v that zeros out subdiagonal
-// entries in the j-th column of a matrix.
-func reflector(v, col []float64, j int) {
- n := len(col)
- if len(v) != n {
- panic("slice length mismatch")
- }
- if j < 0 || n <= j {
- panic("invalid column index")
- }
-
- for i := range v {
- v[i] = 0
- }
- if j == n-1 {
- return
- }
- s := floats.Norm(col[j:], 2)
- if s == 0 {
- return
- }
- v[j] = col[j] + math.Copysign(s, col[j])
- copy(v[j+1:], col[j+1:])
- s = floats.Norm(v[j:], 2)
- floats.Scale(1/s, v[j:])
-}
-
-// applyReflector computes Q*H where H is a Householder matrix represented by
-// the Householder reflector v.
-func applyReflector(qh blas64.General, q blas64.General, v []float64) {
- n := len(v)
- if qh.Rows != n || qh.Cols != n {
- panic("bad size of qh")
- }
- if q.Rows != n || q.Cols != n {
- panic("bad size of q")
- }
- qv := make([]float64, n)
- blas64.Gemv(blas.NoTrans, 1, q, blas64.Vector{1, v}, 0, blas64.Vector{1, qv})
- for i := 0; i < n; i++ {
- for j := 0; j < n; j++ {
- qh.Data[i*qh.Stride+j] = q.Data[i*q.Stride+j]
- }
- }
- for i := 0; i < n; i++ {
- for j := 0; j < n; j++ {
- qh.Data[i*qh.Stride+j] -= 2 * qv[i] * v[j]
- }
- }
- var norm2 float64
- for _, vi := range v {
- norm2 += vi * vi
- }
- norm2inv := 1 / norm2
- for i := 0; i < n; i++ {
- for j := 0; j < n; j++ {
- qh.Data[i*qh.Stride+j] *= norm2inv
- }
- }
-}
-
-// constructGSVDresults returns the matrices [ 0 R ], D1 and D2 described
-// in the documentation of Dtgsja and Dggsvd3, and the result matrix in
-// the documentation for Dggsvp3.
-func constructGSVDresults(n, p, m, k, l int, a, b blas64.General, alpha, beta []float64) (zeroR, d1, d2 blas64.General) {
- // [ 0 R ]
- zeroR = zeros(k+l, n, n)
- dst := zeroR
- dst.Rows = min(m, k+l)
- dst.Cols = k + l
- dst.Data = zeroR.Data[n-k-l:]
- src := a
- src.Rows = min(m, k+l)
- src.Cols = k + l
- src.Data = a.Data[n-k-l:]
- copyGeneral(dst, src)
- if m < k+l {
- // [ 0 R ]
- dst.Rows = k + l - m
- dst.Cols = k + l - m
- dst.Data = zeroR.Data[m*zeroR.Stride+n-(k+l-m):]
- src = b
- src.Rows = k + l - m
- src.Cols = k + l - m
- src.Data = b.Data[(m-k)*b.Stride+n+m-k-l:]
- copyGeneral(dst, src)
- }
-
- // D1
- d1 = zeros(m, k+l, k+l)
- for i := 0; i < k; i++ {
- d1.Data[i*d1.Stride+i] = 1
- }
- for i := k; i < min(m, k+l); i++ {
- d1.Data[i*d1.Stride+i] = alpha[i]
- }
-
- // D2
- d2 = zeros(p, k+l, k+l)
- for i := 0; i < min(l, m-k); i++ {
- d2.Data[i*d2.Stride+i+k] = beta[k+i]
- }
- for i := m - k; i < l; i++ {
- d2.Data[i*d2.Stride+i+k] = 1
- }
-
- return zeroR, d1, d2
-}
-
-func constructGSVPresults(n, p, m, k, l int, a, b blas64.General) (zeroA, zeroB blas64.General) {
- zeroA = zeros(m, n, n)
- dst := zeroA
- dst.Rows = min(m, k+l)
- dst.Cols = k + l
- dst.Data = zeroA.Data[n-k-l:]
- src := a
- dst.Rows = min(m, k+l)
- src.Cols = k + l
- src.Data = a.Data[n-k-l:]
- copyGeneral(dst, src)
-
- zeroB = zeros(p, n, n)
- dst = zeroB
- dst.Rows = l
- dst.Cols = l
- dst.Data = zeroB.Data[n-l:]
- src = b
- dst.Rows = l
- src.Cols = l
- src.Data = b.Data[n-l:]
- copyGeneral(dst, src)
-
- return zeroA, zeroB
-}
OSDN Git Service
