vendor/gonum.org/v1/gonum/lapack/testlapack/general.go

   1 // Copyright ©2015 The Gonum Authors. All rights reserved.
   2 // Use of this source code is governed by a BSD-style
   3 // license that can be found in the LICENSE file.
   4
   5 package testlapack
   6
   7 import (
   8         "fmt"
   9         "math"
  10         "math/cmplx"
  11         "testing"
  12
  13         "golang.org/x/exp/rand"
  14
  15         "gonum.org/v1/gonum/blas"
  16         "gonum.org/v1/gonum/blas/blas64"
  17         "gonum.org/v1/gonum/floats"
  18         "gonum.org/v1/gonum/lapack"
  19 )
  20
  21 const (
  22         // dlamchE is the machine epsilon. For IEEE this is 2^{-53}.
  23         dlamchE = 1.0 / (1 << 53)
  24         dlamchB = 2
  25         dlamchP = dlamchB * dlamchE
  26         // dlamchS is the smallest normal number. For IEEE this is 2^{-1022}.
  27         dlamchS = 1.0 / (1 << 256) / (1 << 256) / (1 << 256) / (1 << 254)
  28 )
  29
  30 func max(a, b int) int {
  31         if a > b {
  32                 return a
  33         }
  34         return b
  35 }
  36
  37 func min(a, b int) int {
  38         if a < b {
  39                 return a
  40         }
  41         return b
  42 }
  43
  44 // worklen describes how much workspace a test should use.
  45 type worklen int
  46
  47 const (
  48         minimumWork worklen = iota
  49         mediumWork
  50         optimumWork
  51 )
  52
  53 // nanSlice allocates a new slice of length n filled with NaN.
  54 func nanSlice(n int) []float64 {
  55         s := make([]float64, n)
  56         for i := range s {
  57                 s[i] = math.NaN()
  58         }
  59         return s
  60 }
  61
  62 // randomSlice allocates a new slice of length n filled with random values.
  63 func randomSlice(n int, rnd *rand.Rand) []float64 {
  64         s := make([]float64, n)
  65         for i := range s {
  66                 s[i] = rnd.NormFloat64()
  67         }
  68         return s
  69 }
  70
  71 // nanGeneral allocates a new r×c general matrix filled with NaN values.
  72 func nanGeneral(r, c, stride int) blas64.General {
  73         if r < 0 || c < 0 {
  74                 panic("bad matrix size")
  75         }
  76         if r == 0 || c == 0 {
  77                 return blas64.General{Stride: max(1, stride)}
  78         }
  79         if stride < c {
  80                 panic("bad stride")
  81         }
  82         return blas64.General{
  83                 Rows:   r,
  84                 Cols:   c,
  85                 Stride: stride,
  86                 Data:   nanSlice((r-1)*stride + c),
  87         }
  88 }
  89
  90 // randomGeneral allocates a new r×c general matrix filled with random
  91 // numbers. Out-of-range elements are filled with NaN values.
  92 func randomGeneral(r, c, stride int, rnd *rand.Rand) blas64.General {
  93         ans := nanGeneral(r, c, stride)
  94         for i := 0; i < r; i++ {
  95                 for j := 0; j < c; j++ {
  96                         ans.Data[i*ans.Stride+j] = rnd.NormFloat64()
  97                 }
  98         }
  99         return ans
 100 }
 101
 102 // randomHessenberg allocates a new n×n Hessenberg matrix filled with zeros
 103 // under the first subdiagonal and with random numbers elsewhere. Out-of-range
 104 // elements are filled with NaN values.
 105 func randomHessenberg(n, stride int, rnd *rand.Rand) blas64.General {
 106         ans := nanGeneral(n, n, stride)
 107         for i := 0; i < n; i++ {
 108                 for j := 0; j < i-1; j++ {
 109                         ans.Data[i*ans.Stride+j] = 0
 110                 }
 111                 for j := max(0, i-1); j < n; j++ {
 112                         ans.Data[i*ans.Stride+j] = rnd.NormFloat64()
 113                 }
 114         }
 115         return ans
 116 }
 117
 118 // randomSchurCanonical returns a random, general matrix in Schur canonical
 119 // form, that is, block upper triangular with 1×1 and 2×2 diagonal blocks where
 120 // each 2×2 diagonal block has its diagonal elements equal and its off-diagonal
 121 // elements of opposite sign.
 122 func randomSchurCanonical(n, stride int, rnd *rand.Rand) blas64.General {
 123         t := randomGeneral(n, n, stride, rnd)
 124         // Zero out the lower triangle.
 125         for i := 0; i < t.Rows; i++ {
 126                 for j := 0; j < i; j++ {
 127                         t.Data[i*t.Stride+j] = 0
 128                 }
 129         }
 130         // Randomly create 2×2 diagonal blocks.
 131         for i := 0; i < t.Rows; {
 132                 if i == t.Rows-1 || rnd.Float64() < 0.5 {
 133                         // 1×1 block.
 134                         i++
 135                         continue
 136                 }
 137                 // 2×2 block.
 138                 // Diagonal elements equal.
 139                 t.Data[(i+1)*t.Stride+i+1] = t.Data[i*t.Stride+i]
 140                 // Off-diagonal elements of opposite sign.
 141                 c := rnd.NormFloat64()
 142                 if math.Signbit(c) == math.Signbit(t.Data[i*t.Stride+i+1]) {
 143                         c *= -1
 144                 }
 145                 t.Data[(i+1)*t.Stride+i] = c
 146                 i += 2
 147         }
 148         return t
 149 }
 150
 151 // blockedUpperTriGeneral returns a normal random, general matrix in the form
 152 //
 153 //            c-k-l  k    l
 154 //  A =    k [  0   A12  A13 ] if r-k-l >= 0;
 155 //         l [  0    0   A23 ]
 156 //     r-k-l [  0    0    0  ]
 157 //
 158 //          c-k-l  k    l
 159 //  A =  k [  0   A12  A13 ] if r-k-l < 0;
 160 //     r-k [  0    0   A23 ]
 161 //
 162 // where the k×k matrix A12 and l×l matrix is non-singular
 163 // upper triangular. A23 is l×l upper triangular if r-k-l >= 0,
 164 // otherwise A23 is (r-k)×l upper trapezoidal.
 165 func blockedUpperTriGeneral(r, c, k, l, stride int, kblock bool, rnd *rand.Rand) blas64.General {
 166         t := l
 167         if kblock {
 168                 t += k
 169         }
 170         ans := zeros(r, c, stride)
 171         for i := 0; i < min(r, t); i++ {
 172                 var v float64
 173                 for v == 0 {
 174                         v = rnd.NormFloat64()
 175                 }
 176                 ans.Data[i*ans.Stride+i+(c-t)] = v
 177         }
 178         for i := 0; i < min(r, t); i++ {
 179                 for j := i + (c - t) + 1; j < c; j++ {
 180                         ans.Data[i*ans.Stride+j] = rnd.NormFloat64()
 181                 }
 182         }
 183         return ans
 184 }
 185
 186 // nanTriangular allocates a new r×c triangular matrix filled with NaN values.
 187 func nanTriangular(uplo blas.Uplo, n, stride int) blas64.Triangular {
 188         if n < 0 {
 189                 panic("bad matrix size")
 190         }
 191         if n == 0 {
 192                 return blas64.Triangular{
 193                         Stride: max(1, stride),
 194                         Uplo:   uplo,
 195                         Diag:   blas.NonUnit,
 196                 }
 197         }
 198         if stride < n {
 199                 panic("bad stride")
 200         }
 201         return blas64.Triangular{
 202                 N:      n,
 203                 Stride: stride,
 204                 Data:   nanSlice((n-1)*stride + n),
 205                 Uplo:   uplo,
 206                 Diag:   blas.NonUnit,
 207         }
 208 }
 209
 210 // generalOutsideAllNaN returns whether all out-of-range elements have NaN
 211 // values.
 212 func generalOutsideAllNaN(a blas64.General) bool {
 213         // Check after last column.
 214         for i := 0; i < a.Rows-1; i++ {
 215                 for _, v := range a.Data[i*a.Stride+a.Cols : i*a.Stride+a.Stride] {
 216                         if !math.IsNaN(v) {
 217                                 return false
 218                         }
 219                 }
 220         }
 221         // Check after last element.
 222         last := (a.Rows-1)*a.Stride + a.Cols
 223         if a.Rows == 0 || a.Cols == 0 {
 224                 last = 0
 225         }
 226         for _, v := range a.Data[last:] {
 227                 if !math.IsNaN(v) {
 228                         return false
 229                 }
 230         }
 231         return true
 232 }
 233
 234 // triangularOutsideAllNaN returns whether all out-of-triangle elements have NaN
 235 // values.
 236 func triangularOutsideAllNaN(a blas64.Triangular) bool {
 237         if a.Uplo == blas.Upper {
 238                 // Check below diagonal.
 239                 for i := 0; i < a.N; i++ {
 240                         for _, v := range a.Data[i*a.Stride : i*a.Stride+i] {
 241                                 if !math.IsNaN(v) {
 242                                         return false
 243                                 }
 244                         }
 245                 }
 246                 // Check after last column.
 247                 for i := 0; i < a.N-1; i++ {
 248                         for _, v := range a.Data[i*a.Stride+a.N : i*a.Stride+a.Stride] {
 249                                 if !math.IsNaN(v) {
 250                                         return false
 251                                 }
 252                         }
 253                 }
 254         } else {
 255                 // Check above diagonal.
 256                 for i := 0; i < a.N-1; i++ {
 257                         for _, v := range a.Data[i*a.Stride+i+1 : i*a.Stride+a.Stride] {
 258                                 if !math.IsNaN(v) {
 259                                         return false
 260                                 }
 261                         }
 262                 }
 263         }
 264         // Check after last element.
 265         for _, v := range a.Data[max(0, a.N-1)*a.Stride+a.N:] {
 266                 if !math.IsNaN(v) {
 267                         return false
 268                 }
 269         }
 270         return true
 271 }
 272
 273 // transposeGeneral returns a new general matrix that is the transpose of the
 274 // input. Nothing is done with data outside the {rows, cols} limit of the general.
 275 func transposeGeneral(a blas64.General) blas64.General {
 276         ans := blas64.General{
 277                 Rows:   a.Cols,
 278                 Cols:   a.Rows,
 279                 Stride: a.Rows,
 280                 Data:   make([]float64, a.Cols*a.Rows),
 281         }
 282         for i := 0; i < a.Rows; i++ {
 283                 for j := 0; j < a.Cols; j++ {
 284                         ans.Data[j*ans.Stride+i] = a.Data[i*a.Stride+j]
 285                 }
 286         }
 287         return ans
 288 }
 289
 290 // columnNorms returns the column norms of a.
 291 func columnNorms(m, n int, a []float64, lda int) []float64 {
 292         bi := blas64.Implementation()
 293         norms := make([]float64, n)
 294         for j := 0; j < n; j++ {
 295                 norms[j] = bi.Dnrm2(m, a[j:], lda)
 296         }
 297         return norms
 298 }
 299
 300 // extractVMat collects the single reflectors from a into a matrix.
 301 func extractVMat(m, n int, a []float64, lda int, direct lapack.Direct, store lapack.StoreV) blas64.General {
 302         k := min(m, n)
 303         switch {
 304         default:
 305                 panic("not implemented")
 306         case direct == lapack.Forward && store == lapack.ColumnWise:
 307                 v := blas64.General{
 308                         Rows:   m,
 309                         Cols:   k,
 310                         Stride: k,
 311                         Data:   make([]float64, m*k),
 312                 }
 313                 for i := 0; i < k; i++ {
 314                         for j := 0; j < i; j++ {
 315                                 v.Data[j*v.Stride+i] = 0
 316                         }
 317                         v.Data[i*v.Stride+i] = 1
 318                         for j := i + 1; j < m; j++ {
 319                                 v.Data[j*v.Stride+i] = a[j*lda+i]
 320                         }
 321                 }
 322                 return v
 323         case direct == lapack.Forward && store == lapack.RowWise:
 324                 v := blas64.General{
 325                         Rows:   k,
 326                         Cols:   n,
 327                         Stride: n,
 328                         Data:   make([]float64, k*n),
 329                 }
 330                 for i := 0; i < k; i++ {
 331                         for j := 0; j < i; j++ {
 332                                 v.Data[i*v.Stride+j] = 0
 333                         }
 334                         v.Data[i*v.Stride+i] = 1
 335                         for j := i + 1; j < n; j++ {
 336                                 v.Data[i*v.Stride+j] = a[i*lda+j]
 337                         }
 338                 }
 339                 return v
 340         }
 341 }
 342
 343 // constructBidiagonal constructs a bidiagonal matrix with the given diagonal
 344 // and off-diagonal elements.
 345 func constructBidiagonal(uplo blas.Uplo, n int, d, e []float64) blas64.General {
 346         bMat := blas64.General{
 347                 Rows:   n,
 348                 Cols:   n,
 349                 Stride: n,
 350                 Data:   make([]float64, n*n),
 351         }
 352
 353         for i := 0; i < n-1; i++ {
 354                 bMat.Data[i*bMat.Stride+i] = d[i]
 355                 if uplo == blas.Upper {
 356                         bMat.Data[i*bMat.Stride+i+1] = e[i]
 357                 } else {
 358                         bMat.Data[(i+1)*bMat.Stride+i] = e[i]
 359                 }
 360         }
 361         bMat.Data[(n-1)*bMat.Stride+n-1] = d[n-1]
 362         return bMat
 363 }
 364
 365 // constructVMat transforms the v matrix based on the storage.
 366 func constructVMat(vMat blas64.General, store lapack.StoreV, direct lapack.Direct) blas64.General {
 367         m := vMat.Rows
 368         k := vMat.Cols
 369         switch {
 370         default:
 371                 panic("not implemented")
 372         case store == lapack.ColumnWise && direct == lapack.Forward:
 373                 ldv := k
 374                 v := make([]float64, m*k)
 375                 for i := 0; i < m; i++ {
 376                         for j := 0; j < k; j++ {
 377                                 if j > i {
 378                                         v[i*ldv+j] = 0
 379                                 } else if j == i {
 380                                         v[i*ldv+i] = 1
 381                                 } else {
 382                                         v[i*ldv+j] = vMat.Data[i*vMat.Stride+j]
 383                                 }
 384                         }
 385                 }
 386                 return blas64.General{
 387                         Rows:   m,
 388                         Cols:   k,
 389                         Stride: k,
 390                         Data:   v,
 391                 }
 392         case store == lapack.RowWise && direct == lapack.Forward:
 393                 ldv := m
 394                 v := make([]float64, m*k)
 395                 for i := 0; i < m; i++ {
 396                         for j := 0; j < k; j++ {
 397                                 if j > i {
 398                                         v[j*ldv+i] = 0
 399                                 } else if j == i {
 400                                         v[j*ldv+i] = 1
 401                                 } else {
 402                                         v[j*ldv+i] = vMat.Data[i*vMat.Stride+j]
 403                                 }
 404                         }
 405                 }
 406                 return blas64.General{
 407                         Rows:   k,
 408                         Cols:   m,
 409                         Stride: m,
 410                         Data:   v,
 411                 }
 412         case store == lapack.ColumnWise && direct == lapack.Backward:
 413                 rowsv := m
 414                 ldv := k
 415                 v := make([]float64, m*k)
 416                 for i := 0; i < m; i++ {
 417                         for j := 0; j < k; j++ {
 418                                 vrow := rowsv - i - 1
 419                                 vcol := k - j - 1
 420                                 if j > i {
 421                                         v[vrow*ldv+vcol] = 0
 422                                 } else if j == i {
 423                                         v[vrow*ldv+vcol] = 1
 424                                 } else {
 425                                         v[vrow*ldv+vcol] = vMat.Data[i*vMat.Stride+j]
 426                                 }
 427                         }
 428                 }
 429                 return blas64.General{
 430                         Rows:   rowsv,
 431                         Cols:   ldv,
 432                         Stride: ldv,
 433                         Data:   v,
 434                 }
 435         case store == lapack.RowWise && direct == lapack.Backward:
 436                 rowsv := k
 437                 ldv := m
 438                 v := make([]float64, m*k)
 439                 for i := 0; i < m; i++ {
 440                         for j := 0; j < k; j++ {
 441                                 vcol := ldv - i - 1
 442                                 vrow := k - j - 1
 443                                 if j > i {
 444                                         v[vrow*ldv+vcol] = 0
 445                                 } else if j == i {
 446                                         v[vrow*ldv+vcol] = 1
 447                                 } else {
 448                                         v[vrow*ldv+vcol] = vMat.Data[i*vMat.Stride+j]
 449                                 }
 450                         }
 451                 }
 452                 return blas64.General{
 453                         Rows:   rowsv,
 454                         Cols:   ldv,
 455                         Stride: ldv,
 456                         Data:   v,
 457                 }
 458         }
 459 }
 460
 461 func constructH(tau []float64, v blas64.General, store lapack.StoreV, direct lapack.Direct) blas64.General {
 462         m := v.Rows
 463         k := v.Cols
 464         if store == lapack.RowWise {
 465                 m, k = k, m
 466         }
 467         h := blas64.General{
 468                 Rows:   m,
 469                 Cols:   m,
 470                 Stride: m,
 471                 Data:   make([]float64, m*m),
 472         }
 473         for i := 0; i < m; i++ {
 474                 h.Data[i*m+i] = 1
 475         }
 476         for i := 0; i < k; i++ {
 477                 vecData := make([]float64, m)
 478                 if store == lapack.ColumnWise {
 479                         for j := 0; j < m; j++ {
 480                                 vecData[j] = v.Data[j*v.Cols+i]
 481                         }
 482                 } else {
 483                         for j := 0; j < m; j++ {
 484                                 vecData[j] = v.Data[i*v.Cols+j]
 485                         }
 486                 }
 487                 vec := blas64.Vector{
 488                         Inc:  1,
 489                         Data: vecData,
 490                 }
 491
 492                 hi := blas64.General{
 493                         Rows:   m,
 494                         Cols:   m,
 495                         Stride: m,
 496                         Data:   make([]float64, m*m),
 497                 }
 498                 for i := 0; i < m; i++ {
 499                         hi.Data[i*m+i] = 1
 500                 }
 501                 // hi = I - tau * v * v^T
 502                 blas64.Ger(-tau[i], vec, vec, hi)
 503
 504                 hcopy := blas64.General{
 505                         Rows:   m,
 506                         Cols:   m,
 507                         Stride: m,
 508                         Data:   make([]float64, m*m),
 509                 }
 510                 copy(hcopy.Data, h.Data)
 511                 if direct == lapack.Forward {
 512                         // H = H * H_I in forward mode
 513                         blas64.Gemm(blas.NoTrans, blas.NoTrans, 1, hcopy, hi, 0, h)
 514                 } else {
 515                         // H = H_I * H in backward mode
 516                         blas64.Gemm(blas.NoTrans, blas.NoTrans, 1, hi, hcopy, 0, h)
 517                 }
 518         }
 519         return h
 520 }
 521
 522 // constructQ constructs the Q matrix from the result of dgeqrf and dgeqr2.
 523 func constructQ(kind string, m, n int, a []float64, lda int, tau []float64) blas64.General {
 524         k := min(m, n)
 525         return constructQK(kind, m, n, k, a, lda, tau)
 526 }
 527
 528 // constructQK constructs the Q matrix from the result of dgeqrf and dgeqr2 using
 529 // the first k reflectors.
 530 func constructQK(kind string, m, n, k int, a []float64, lda int, tau []float64) blas64.General {
 531         var sz int
 532         switch kind {
 533         case "QR":
 534                 sz = m
 535         case "LQ", "RQ":
 536                 sz = n
 537         }
 538
 539         q := blas64.General{
 540                 Rows:   sz,
 541                 Cols:   sz,
 542                 Stride: sz,
 543                 Data:   make([]float64, sz*sz),
 544         }
 545         for i := 0; i < sz; i++ {
 546                 q.Data[i*sz+i] = 1
 547         }
 548         qCopy := blas64.General{
 549                 Rows:   q.Rows,
 550                 Cols:   q.Cols,
 551                 Stride: q.Stride,
 552                 Data:   make([]float64, len(q.Data)),
 553         }
 554         for i := 0; i < k; i++ {
 555                 h := blas64.General{
 556                         Rows:   sz,
 557                         Cols:   sz,
 558                         Stride: sz,
 559                         Data:   make([]float64, sz*sz),
 560                 }
 561                 for j := 0; j < sz; j++ {
 562                         h.Data[j*sz+j] = 1
 563                 }
 564                 vVec := blas64.Vector{
 565                         Inc:  1,
 566                         Data: make([]float64, sz),
 567                 }
 568                 switch kind {
 569                 case "QR":
 570                         vVec.Data[i] = 1
 571                         for j := i + 1; j < sz; j++ {
 572                                 vVec.Data[j] = a[lda*j+i]
 573                         }
 574                 case "LQ":
 575                         vVec.Data[i] = 1
 576                         for j := i + 1; j < sz; j++ {
 577                                 vVec.Data[j] = a[i*lda+j]
 578                         }
 579                 case "RQ":
 580                         for j := 0; j < n-k+i; j++ {
 581                                 vVec.Data[j] = a[(m-k+i)*lda+j]
 582                         }
 583                         vVec.Data[n-k+i] = 1
 584                 }
 585                 blas64.Ger(-tau[i], vVec, vVec, h)
 586                 copy(qCopy.Data, q.Data)
 587                 // Multiply q by the new h.
 588                 switch kind {
 589                 case "QR", "RQ":
 590                         blas64.Gemm(blas.NoTrans, blas.NoTrans, 1, qCopy, h, 0, q)
 591                 case "LQ":
 592                         blas64.Gemm(blas.NoTrans, blas.NoTrans, 1, h, qCopy, 0, q)
 593                 }
 594         }
 595         return q
 596 }
 597
 598 // checkBidiagonal checks the bidiagonal decomposition from dlabrd and dgebd2.
 599 // The input to this function is the answer returned from the routines, stored
 600 // in a, d, e, tauP, and tauQ. The data of original A matrix (before
 601 // decomposition) is input in aCopy.
 602 //
 603 // checkBidiagonal constructs the V and U matrices, and from them constructs Q
 604 // and P. Using these constructions, it checks that Q^T * A * P and checks that
 605 // the result is bidiagonal.
 606 func checkBidiagonal(t *testing.T, m, n, nb int, a []float64, lda int, d, e, tauP, tauQ, aCopy []float64) {
 607         // Check the answer.
 608         // Construct V and U.
 609         qMat := constructQPBidiagonal(lapack.ApplyQ, m, n, nb, a, lda, tauQ)
 610         pMat := constructQPBidiagonal(lapack.ApplyP, m, n, nb, a, lda, tauP)
 611
 612         // Compute Q^T * A * P.
 613         aMat := blas64.General{
 614                 Rows:   m,
 615                 Cols:   n,
 616                 Stride: lda,
 617                 Data:   make([]float64, len(aCopy)),
 618         }
 619         copy(aMat.Data, aCopy)
 620
 621         tmp1 := blas64.General{
 622                 Rows:   m,
 623                 Cols:   n,
 624                 Stride: n,
 625                 Data:   make([]float64, m*n),
 626         }
 627         blas64.Gemm(blas.Trans, blas.NoTrans, 1, qMat, aMat, 0, tmp1)
 628         tmp2 := blas64.General{
 629                 Rows:   m,
 630                 Cols:   n,
 631                 Stride: n,
 632                 Data:   make([]float64, m*n),
 633         }
 634         blas64.Gemm(blas.NoTrans, blas.NoTrans, 1, tmp1, pMat, 0, tmp2)
 635
 636         // Check that the first nb rows and cols of tm2 are upper bidiagonal
 637         // if m >= n, and lower bidiagonal otherwise.
 638         correctDiag := true
 639         matchD := true
 640         matchE := true
 641         for i := 0; i < m; i++ {
 642                 for j := 0; j < n; j++ {
 643                         if i >= nb && j >= nb {
 644                                 continue
 645                         }
 646                         v := tmp2.Data[i*tmp2.Stride+j]
 647                         if i == j {
 648                                 if math.Abs(d[i]-v) > 1e-12 {
 649                                         matchD = false
 650                                 }
 651                                 continue
 652                         }
 653                         if m >= n && i == j-1 {
 654                                 if math.Abs(e[j-1]-v) > 1e-12 {
 655                                         matchE = false
 656                                 }
 657                                 continue
 658                         }
 659                         if m < n && i-1 == j {
 660                                 if math.Abs(e[i-1]-v) > 1e-12 {
 661                                         matchE = false
 662                                 }
 663                                 continue
 664                         }
 665                         if math.Abs(v) > 1e-12 {
 666                                 correctDiag = false
 667                         }
 668                 }
 669         }
 670         if !correctDiag {
 671                 t.Errorf("Updated A not bi-diagonal")
 672         }
 673         if !matchD {
 674                 fmt.Println("d = ", d)
 675                 t.Errorf("D Mismatch")
 676         }
 677         if !matchE {
 678                 t.Errorf("E mismatch")
 679         }
 680 }
 681
 682 // constructQPBidiagonal constructs Q or P from the Bidiagonal decomposition
 683 // computed by dlabrd and bgebd2.
 684 func constructQPBidiagonal(vect lapack.DecompUpdate, m, n, nb int, a []float64, lda int, tau []float64) blas64.General {
 685         sz := n
 686         if vect == lapack.ApplyQ {
 687                 sz = m
 688         }
 689
 690         var ldv int
 691         var v blas64.General
 692         if vect == lapack.ApplyQ {
 693                 ldv = nb
 694                 v = blas64.General{
 695                         Rows:   m,
 696                         Cols:   nb,
 697                         Stride: ldv,
 698                         Data:   make([]float64, m*ldv),
 699                 }
 700         } else {
 701                 ldv = n
 702                 v = blas64.General{
 703                         Rows:   nb,
 704                         Cols:   n,
 705                         Stride: ldv,
 706                         Data:   make([]float64, m*ldv),
 707                 }
 708         }
 709
 710         if vect == lapack.ApplyQ {
 711                 if m >= n {
 712                         for i := 0; i < m; i++ {
 713                                 for j := 0; j <= min(nb-1, i); j++ {
 714                                         if i == j {
 715                                                 v.Data[i*ldv+j] = 1
 716                                                 continue
 717                                         }
 718                                         v.Data[i*ldv+j] = a[i*lda+j]
 719                                 }
 720                         }
 721                 } else {
 722                         for i := 1; i < m; i++ {
 723                                 for j := 0; j <= min(nb-1, i-1); j++ {
 724                                         if i-1 == j {
 725                                                 v.Data[i*ldv+j] = 1
 726                                                 continue
 727                                         }
 728                                         v.Data[i*ldv+j] = a[i*lda+j]
 729                                 }
 730                         }
 731                 }
 732         } else {
 733                 if m < n {
 734                         for i := 0; i < nb; i++ {
 735                                 for j := i; j < n; j++ {
 736                                         if i == j {
 737                                                 v.Data[i*ldv+j] = 1
 738                                                 continue
 739                                         }
 740                                         v.Data[i*ldv+j] = a[i*lda+j]
 741                                 }
 742                         }
 743                 } else {
 744                         for i := 0; i < nb; i++ {
 745                                 for j := i + 1; j < n; j++ {
 746                                         if j-1 == i {
 747                                                 v.Data[i*ldv+j] = 1
 748                                                 continue
 749                                         }
 750                                         v.Data[i*ldv+j] = a[i*lda+j]
 751                                 }
 752                         }
 753                 }
 754         }
 755
 756         // The variable name is a computation of Q, but the algorithm is mostly the
 757         // same for computing P (just with different data).
 758         qMat := blas64.General{
 759                 Rows:   sz,
 760                 Cols:   sz,
 761                 Stride: sz,
 762                 Data:   make([]float64, sz*sz),
 763         }
 764         hMat := blas64.General{
 765                 Rows:   sz,
 766                 Cols:   sz,
 767                 Stride: sz,
 768                 Data:   make([]float64, sz*sz),
 769         }
 770         // set Q to I
 771         for i := 0; i < sz; i++ {
 772                 qMat.Data[i*qMat.Stride+i] = 1
 773         }
 774         for i := 0; i < nb; i++ {
 775                 qCopy := blas64.General{Rows: qMat.Rows, Cols: qMat.Cols, Stride: qMat.Stride, Data: make([]float64, len(qMat.Data))}
 776                 copy(qCopy.Data, qMat.Data)
 777
 778                 // Set g and h to I
 779                 for i := 0; i < sz; i++ {
 780                         for j := 0; j < sz; j++ {
 781                                 if i == j {
 782                                         hMat.Data[i*sz+j] = 1
 783                                 } else {
 784                                         hMat.Data[i*sz+j] = 0
 785                                 }
 786                         }
 787                 }
 788                 var vi blas64.Vector
 789                 // H -= tauQ[i] * v[i] * v[i]^t
 790                 if vect == lapack.ApplyQ {
 791                         vi = blas64.Vector{
 792                                 Inc:  v.Stride,
 793                                 Data: v.Data[i:],
 794                         }
 795                 } else {
 796                         vi = blas64.Vector{
 797                                 Inc:  1,
 798                                 Data: v.Data[i*v.Stride:],
 799                         }
 800                 }
 801                 blas64.Ger(-tau[i], vi, vi, hMat)
 802                 // Q = Q * G[1]
 803                 blas64.Gemm(blas.NoTrans, blas.NoTrans, 1, qCopy, hMat, 0, qMat)
 804         }
 805         return qMat
 806 }
 807
 808 // printRowise prints the matrix with one row per line. This is useful for debugging.
 809 // If beyond is true, it prints beyond the final column to lda. If false, only
 810 // the columns are printed.
 811 func printRowise(a []float64, m, n, lda int, beyond bool) {
 812         for i := 0; i < m; i++ {
 813                 end := n
 814                 if beyond {
 815                         end = lda
 816                 }
 817                 fmt.Println(a[i*lda : i*lda+end])
 818         }
 819 }
 820
 821 // isOrthonormal checks that a general matrix is orthonormal.
 822 func isOrthonormal(q blas64.General) bool {
 823         n := q.Rows
 824         for i := 0; i < n; i++ {
 825                 for j := i; j < n; j++ {
 826                         dot := blas64.Dot(n,
 827                                 blas64.Vector{Inc: 1, Data: q.Data[i*q.Stride:]},
 828                                 blas64.Vector{Inc: 1, Data: q.Data[j*q.Stride:]},
 829                         )
 830                         if math.IsNaN(dot) {
 831                                 return false
 832                         }
 833                         if i == j {
 834                                 if math.Abs(dot-1) > 1e-10 {
 835                                         return false
 836                                 }
 837                         } else {
 838                                 if math.Abs(dot) > 1e-10 {
 839                                         return false
 840                                 }
 841                         }
 842                 }
 843         }
 844         return true
 845 }
 846
 847 // copyMatrix copies an m×n matrix src of stride n into an m×n matrix dst of stride ld.
 848 func copyMatrix(m, n int, dst []float64, ld int, src []float64) {
 849         for i := 0; i < m; i++ {
 850                 copy(dst[i*ld:i*ld+n], src[i*n:i*n+n])
 851         }
 852 }
 853
 854 func copyGeneral(dst, src blas64.General) {
 855         r := min(dst.Rows, src.Rows)
 856         c := min(dst.Cols, src.Cols)
 857         for i := 0; i < r; i++ {
 858                 copy(dst.Data[i*dst.Stride:i*dst.Stride+c], src.Data[i*src.Stride:i*src.Stride+c])
 859         }
 860 }
 861
 862 // cloneGeneral allocates and returns an exact copy of the given general matrix.
 863 func cloneGeneral(a blas64.General) blas64.General {
 864         c := a
 865         c.Data = make([]float64, len(a.Data))
 866         copy(c.Data, a.Data)
 867         return c
 868 }
 869
 870 // equalApprox returns whether the matrices A and B of order n are approximately
 871 // equal within given tolerance.
 872 func equalApprox(m, n int, a []float64, lda int, b []float64, tol float64) bool {
 873         for i := 0; i < m; i++ {
 874                 for j := 0; j < n; j++ {
 875                         diff := a[i*lda+j] - b[i*n+j]
 876                         if math.IsNaN(diff) || math.Abs(diff) > tol {
 877                                 return false
 878                         }
 879                 }
 880         }
 881         return true
 882 }
 883
 884 // equalApproxGeneral returns whether the general matrices a and b are
 885 // approximately equal within given tolerance.
 886 func equalApproxGeneral(a, b blas64.General, tol float64) bool {
 887         if a.Rows != b.Rows || a.Cols != b.Cols {
 888                 panic("bad input")
 889         }
 890         for i := 0; i < a.Rows; i++ {
 891                 for j := 0; j < a.Cols; j++ {
 892                         diff := a.Data[i*a.Stride+j] - b.Data[i*b.Stride+j]
 893                         if math.IsNaN(diff) || math.Abs(diff) > tol {
 894                                 return false
 895                         }
 896                 }
 897         }
 898         return true
 899 }
 900
 901 // equalApproxTriangular returns whether the triangular matrices A and B of
 902 // order n are approximately equal within given tolerance.
 903 func equalApproxTriangular(upper bool, n int, a []float64, lda int, b []float64, tol float64) bool {
 904         if upper {
 905                 for i := 0; i < n; i++ {
 906                         for j := i; j < n; j++ {
 907                                 diff := a[i*lda+j] - b[i*n+j]
 908                                 if math.IsNaN(diff) || math.Abs(diff) > tol {
 909                                         return false
 910                                 }
 911                         }
 912                 }
 913                 return true
 914         }
 915         for i := 0; i < n; i++ {
 916                 for j := 0; j <= i; j++ {
 917                         diff := a[i*lda+j] - b[i*n+j]
 918                         if math.IsNaN(diff) || math.Abs(diff) > tol {
 919                                 return false
 920                         }
 921                 }
 922         }
 923         return true
 924 }
 925
 926 func equalApproxSymmetric(a, b blas64.Symmetric, tol float64) bool {
 927         if a.Uplo != b.Uplo {
 928                 return false
 929         }
 930         if a.N != b.N {
 931                 return false
 932         }
 933         if a.Uplo == blas.Upper {
 934                 for i := 0; i < a.N; i++ {
 935                         for j := i; j < a.N; j++ {
 936                                 if !floats.EqualWithinAbsOrRel(a.Data[i*a.Stride+j], b.Data[i*b.Stride+j], tol, tol) {
 937                                         return false
 938                                 }
 939                         }
 940                 }
 941                 return true
 942         }
 943         for i := 0; i < a.N; i++ {
 944                 for j := 0; j <= i; j++ {
 945                         if !floats.EqualWithinAbsOrRel(a.Data[i*a.Stride+j], b.Data[i*b.Stride+j], tol, tol) {
 946                                 return false
 947                         }
 948                 }
 949         }
 950         return true
 951 }
 952
 953 // randSymBand creates a random symmetric banded matrix, and returns both the
 954 // random matrix and the equivalent Symmetric matrix for testing. rnder
 955 // specifies the random number
 956 func randSymBand(ul blas.Uplo, n, ldab, kb int, rnd *rand.Rand) (blas64.Symmetric, blas64.SymmetricBand) {
 957         // A matrix is positive definite if and only if it has a Cholesky
 958         // decomposition. Generate a random banded lower triangular matrix
 959         // to construct the random symmetric matrix.
 960         a := make([]float64, n*n)
 961         for i := 0; i < n; i++ {
 962                 for j := max(0, i-kb); j <= i; j++ {
 963                         a[i*n+j] = rnd.NormFloat64()
 964                 }
 965                 a[i*n+i] = math.Abs(a[i*n+i])
 966                 // Add an extra amound to the diagonal in order to improve the condition number.
 967                 a[i*n+i] += 1.5 * rnd.Float64()
 968         }
 969         agen := blas64.General{
 970                 Rows:   n,
 971                 Cols:   n,
 972                 Stride: n,
 973                 Data:   a,
 974         }
 975
 976         // Construct the SymDense from a*a^T
 977         c := make([]float64, n*n)
 978         cgen := blas64.General{
 979                 Rows:   n,
 980                 Cols:   n,
 981                 Stride: n,
 982                 Data:   c,
 983         }
 984         blas64.Gemm(blas.NoTrans, blas.Trans, 1, agen, agen, 0, cgen)
 985         sym := blas64.Symmetric{
 986                 N:      n,
 987                 Stride: n,
 988                 Data:   c,
 989                 Uplo:   ul,
 990         }
 991
 992         b := symToSymBand(ul, c, n, n, kb, ldab)
 993         band := blas64.SymmetricBand{
 994                 N:      n,
 995                 K:      kb,
 996                 Stride: ldab,
 997                 Data:   b,
 998                 Uplo:   ul,
 999         }
1000
1001         return sym, band
1002 }
1003
1004 // symToSymBand takes the data in a Symmetric matrix and returns a
1005 // SymmetricBanded matrix.
1006 func symToSymBand(ul blas.Uplo, a []float64, n, lda, kb, ldab int) []float64 {
1007         if ul == blas.Upper {
1008                 band := make([]float64, (n-1)*ldab+kb+1)
1009                 for i := 0; i < n; i++ {
1010                         for j := i; j < min(i+kb+1, n); j++ {
1011                                 band[i*ldab+j-i] = a[i*lda+j]
1012                         }
1013                 }
1014                 return band
1015         }
1016         band := make([]float64, (n-1)*ldab+kb+1)
1017         for i := 0; i < n; i++ {
1018                 for j := max(0, i-kb); j <= i; j++ {
1019                         band[i*ldab+j-i+kb] = a[i*lda+j]
1020                 }
1021         }
1022         return band
1023 }
1024
1025 // symBandToSym takes a banded symmetric matrix and returns the same data as
1026 // a Symmetric matrix.
1027 func symBandToSym(ul blas.Uplo, band []float64, n, kb, ldab int) blas64.Symmetric {
1028         sym := make([]float64, n*n)
1029         if ul == blas.Upper {
1030                 for i := 0; i < n; i++ {
1031                         for j := 0; j < min(kb+1+i, n)-i; j++ {
1032                                 sym[i*n+i+j] = band[i*ldab+j]
1033                         }
1034                 }
1035         } else {
1036                 for i := 0; i < n; i++ {
1037                         for j := kb - min(i, kb); j < kb+1; j++ {
1038                                 sym[i*n+i-kb+j] = band[i*ldab+j]
1039                         }
1040                 }
1041         }
1042         return blas64.Symmetric{
1043                 N:      n,
1044                 Stride: n,
1045                 Data:   sym,
1046                 Uplo:   ul,
1047         }
1048 }
1049
1050 // eye returns an identity matrix of given order and stride.
1051 func eye(n, stride int) blas64.General {
1052         ans := nanGeneral(n, n, stride)
1053         for i := 0; i < n; i++ {
1054                 for j := 0; j < n; j++ {
1055                         ans.Data[i*ans.Stride+j] = 0
1056                 }
1057                 ans.Data[i*ans.Stride+i] = 1
1058         }
1059         return ans
1060 }
1061
1062 // zeros returns an m×n matrix with given stride filled with zeros.
1063 func zeros(m, n, stride int) blas64.General {
1064         a := nanGeneral(m, n, stride)
1065         for i := 0; i < m; i++ {
1066                 for j := 0; j < n; j++ {
1067                         a.Data[i*a.Stride+j] = 0
1068                 }
1069         }
1070         return a
1071 }
1072
1073 // extract2x2Block returns the elements of T at [0,0], [0,1], [1,0], and [1,1].
1074 func extract2x2Block(t []float64, ldt int) (a, b, c, d float64) {
1075         return t[0], t[1], t[ldt], t[ldt+1]
1076 }
1077
1078 // isSchurCanonical returns whether the 2×2 matrix [a b; c d] is in Schur
1079 // canonical form.
1080 func isSchurCanonical(a, b, c, d float64) bool {
1081         return c == 0 || (a == d && math.Signbit(b) != math.Signbit(c))
1082 }
1083
1084 // isSchurCanonicalGeneral returns whether T is block upper triangular with 1×1
1085 // and 2×2 diagonal blocks, each 2×2 block in Schur canonical form. The function
1086 // checks only along the diagonal and the first subdiagonal, otherwise the lower
1087 // triangle is not accessed.
1088 func isSchurCanonicalGeneral(t blas64.General) bool {
1089         if t.Rows != t.Cols {
1090                 panic("invalid matrix")
1091         }
1092         for i := 0; i < t.Rows-1; {
1093                 if t.Data[(i+1)*t.Stride+i] == 0 {
1094                         // 1×1 block.
1095                         i++
1096                         continue
1097                 }
1098                 // 2×2 block.
1099                 a, b, c, d := extract2x2Block(t.Data[i*t.Stride+i:], t.Stride)
1100                 if !isSchurCanonical(a, b, c, d) {
1101                         return false
1102                 }
1103                 i += 2
1104         }
1105         return true
1106 }
1107
1108 // schurBlockEigenvalues returns the two eigenvalues of the 2×2 matrix [a b; c d]
1109 // that must be in Schur canonical form.
1110 func schurBlockEigenvalues(a, b, c, d float64) (ev1, ev2 complex128) {
1111         if !isSchurCanonical(a, b, c, d) {
1112                 panic("block not in Schur canonical form")
1113         }
1114         if c == 0 {
1115                 return complex(a, 0), complex(d, 0)
1116         }
1117         im := math.Sqrt(-b * c)
1118         return complex(a, im), complex(a, -im)
1119 }
1120
1121 // schurBlockSize returns the size of the diagonal block at i-th row in the
1122 // upper quasi-triangular matrix t in Schur canonical form, and whether i points
1123 // to the first row of the block. For zero-sized matrices the function returns 0
1124 // and true.
1125 func schurBlockSize(t blas64.General, i int) (size int, first bool) {
1126         if t.Rows != t.Cols {
1127                 panic("matrix not square")
1128         }
1129         if t.Rows == 0 {
1130                 return 0, true
1131         }
1132         if i < 0 || t.Rows <= i {
1133                 panic("index out of range")
1134         }
1135
1136         first = true
1137         if i > 0 && t.Data[i*t.Stride+i-1] != 0 {
1138                 // There is a non-zero element to the left, therefore i must
1139                 // point to the second row in a 2×2 diagonal block.
1140                 first = false
1141                 i--
1142         }
1143         size = 1
1144         if i+1 < t.Rows && t.Data[(i+1)*t.Stride+i] != 0 {
1145                 // There is a non-zero element below, this must be a 2×2
1146                 // diagonal block.
1147                 size = 2
1148         }
1149         return size, first
1150 }
1151
1152 // containsComplex returns whether z is approximately equal to one of the complex
1153 // numbers in v. If z is found, its index in v will be also returned.
1154 func containsComplex(v []complex128, z complex128, tol float64) (found bool, index int) {
1155         for i := range v {
1156                 if cmplx.Abs(v[i]-z) < tol {
1157                         return true, i
1158                 }
1159         }
1160         return false, -1
1161 }
1162
1163 // isAllNaN returns whether x contains only NaN values.
1164 func isAllNaN(x []float64) bool {
1165         for _, v := range x {
1166                 if !math.IsNaN(v) {
1167                         return false
1168                 }
1169         }
1170         return true
1171 }
1172
1173 // isUpperHessenberg returns whether h contains only zeros below the
1174 // subdiagonal.
1175 func isUpperHessenberg(h blas64.General) bool {
1176         if h.Rows != h.Cols {
1177                 panic("matrix not square")
1178         }
1179         n := h.Rows
1180         for i := 0; i < n; i++ {
1181                 for j := 0; j < n; j++ {
1182                         if i > j+1 && h.Data[i*h.Stride+j] != 0 {
1183                                 return false
1184                         }
1185                 }
1186         }
1187         return true
1188 }
1189
1190 // isUpperTriangular returns whether a contains only zeros below the diagonal.
1191 func isUpperTriangular(a blas64.General) bool {
1192         n := a.Rows
1193         for i := 1; i < n; i++ {
1194                 for j := 0; j < i; j++ {
1195                         if a.Data[i*a.Stride+j] != 0 {
1196                                 return false
1197                         }
1198                 }
1199         }
1200         return true
1201 }
1202
1203 // unbalancedSparseGeneral returns an m×n dense matrix with a random sparse
1204 // structure consisting of nz nonzero elements. The matrix will be unbalanced by
1205 // multiplying each element randomly by its row or column index.
1206 func unbalancedSparseGeneral(m, n, stride int, nonzeros int, rnd *rand.Rand) blas64.General {
1207         a := zeros(m, n, stride)
1208         for k := 0; k < nonzeros; k++ {
1209                 i := rnd.Intn(n)
1210                 j := rnd.Intn(n)
1211                 if rnd.Float64() < 0.5 {
1212                         a.Data[i*stride+j] = float64(i+1) * rnd.NormFloat64()
1213                 } else {
1214                         a.Data[i*stride+j] = float64(j+1) * rnd.NormFloat64()
1215                 }
1216         }
1217         return a
1218 }
1219
1220 // columnOf returns a copy of the j-th column of a.
1221 func columnOf(a blas64.General, j int) []float64 {
1222         if j < 0 || a.Cols <= j {
1223                 panic("bad column index")
1224         }
1225         col := make([]float64, a.Rows)
1226         for i := range col {
1227                 col[i] = a.Data[i*a.Stride+j]
1228         }
1229         return col
1230 }
1231
1232 // isRightEigenvectorOf returns whether the vector xRe+i*xIm, where i is the
1233 // imaginary unit, is the right eigenvector of A corresponding to the eigenvalue
1234 // lambda.
1235 //
1236 // A right eigenvector corresponding to a complex eigenvalue λ is a complex
1237 // non-zero vector x such that
1238 //  A x = λ x.
1239 func isRightEigenvectorOf(a blas64.General, xRe, xIm []float64, lambda complex128, tol float64) bool {
1240         if a.Rows != a.Cols {
1241                 panic("matrix not square")
1242         }
1243
1244         if imag(lambda) != 0 && xIm == nil {
1245                 // Complex eigenvalue of a real matrix cannot have a real
1246                 // eigenvector.
1247                 return false
1248         }
1249
1250         n := a.Rows
1251
1252         // Compute A real(x) and store the result into xReAns.
1253         xReAns := make([]float64, n)
1254         blas64.Gemv(blas.NoTrans, 1, a, blas64.Vector{1, xRe}, 0, blas64.Vector{1, xReAns})
1255
1256         if imag(lambda) == 0 && xIm == nil {
1257                 // Real eigenvalue and eigenvector.
1258
1259                 // Compute λx and store the result into lambdax.
1260                 lambdax := make([]float64, n)
1261                 floats.AddScaled(lambdax, real(lambda), xRe)
1262
1263                 // This is expressed as the inverse to catch the case
1264                 // xReAns_i = Inf and lambdax_i = Inf of the same sign.
1265                 return !(floats.Distance(xReAns, lambdax, math.Inf(1)) > tol)
1266         }
1267
1268         // Complex eigenvector, and real or complex eigenvalue.
1269
1270         // Compute A imag(x) and store the result into xImAns.
1271         xImAns := make([]float64, n)
1272         blas64.Gemv(blas.NoTrans, 1, a, blas64.Vector{1, xIm}, 0, blas64.Vector{1, xImAns})
1273
1274         // Compute λx and store the result into lambdax.
1275         lambdax := make([]complex128, n)
1276         for i := range lambdax {
1277                 lambdax[i] = lambda * complex(xRe[i], xIm[i])
1278         }
1279
1280         for i, v := range lambdax {
1281                 ax := complex(xReAns[i], xImAns[i])
1282                 if cmplx.Abs(v-ax) > tol {
1283                         return false
1284                 }
1285         }
1286         return true
1287 }
1288
1289 // isLeftEigenvectorOf returns whether the vector yRe+i*yIm, where i is the
1290 // imaginary unit, is the left eigenvector of A corresponding to the eigenvalue
1291 // lambda.
1292 //
1293 // A left eigenvector corresponding to a complex eigenvalue λ is a complex
1294 // non-zero vector y such that
1295 //  y^H A = λ y^H,
1296 // which is equivalent for real A to
1297 //  A^T y = conj(λ) y,
1298 func isLeftEigenvectorOf(a blas64.General, yRe, yIm []float64, lambda complex128, tol float64) bool {
1299         if a.Rows != a.Cols {
1300                 panic("matrix not square")
1301         }
1302
1303         if imag(lambda) != 0 && yIm == nil {
1304                 // Complex eigenvalue of a real matrix cannot have a real
1305                 // eigenvector.
1306                 return false
1307         }
1308
1309         n := a.Rows
1310
1311         // Compute A^T real(y) and store the result into yReAns.
1312         yReAns := make([]float64, n)
1313         blas64.Gemv(blas.Trans, 1, a, blas64.Vector{1, yRe}, 0, blas64.Vector{1, yReAns})
1314
1315         if imag(lambda) == 0 && yIm == nil {
1316                 // Real eigenvalue and eigenvector.
1317
1318                 // Compute λy and store the result into lambday.
1319                 lambday := make([]float64, n)
1320                 floats.AddScaled(lambday, real(lambda), yRe)
1321
1322                 // This is expressed as the inverse to catch the case
1323                 // yReAns_i = Inf and lambday_i = Inf of the same sign.
1324                 return !(floats.Distance(yReAns, lambday, math.Inf(1)) > tol)
1325         }
1326
1327         // Complex eigenvector, and real or complex eigenvalue.
1328
1329         // Compute A^T imag(y) and store the result into yImAns.
1330         yImAns := make([]float64, n)
1331         blas64.Gemv(blas.Trans, 1, a, blas64.Vector{1, yIm}, 0, blas64.Vector{1, yImAns})
1332
1333         // Compute conj(λ)y and store the result into lambday.
1334         lambda = cmplx.Conj(lambda)
1335         lambday := make([]complex128, n)
1336         for i := range lambday {
1337                 lambday[i] = lambda * complex(yRe[i], yIm[i])
1338         }
1339
1340         for i, v := range lambday {
1341                 ay := complex(yReAns[i], yImAns[i])
1342                 if cmplx.Abs(v-ay) > tol {
1343                         return false
1344                 }
1345         }
1346         return true
1347 }
1348
1349 // rootsOfUnity returns the n complex numbers whose n-th power is equal to 1.
1350 func rootsOfUnity(n int) []complex128 {
1351         w := make([]complex128, n)
1352         for i := 0; i < n; i++ {
1353                 angle := math.Pi * float64(2*i) / float64(n)
1354                 w[i] = complex(math.Cos(angle), math.Sin(angle))
1355         }
1356         return w
1357 }
1358
1359 // randomOrthogonal returns an n×n random orthogonal matrix.
1360 func randomOrthogonal(n int, rnd *rand.Rand) blas64.General {
1361         q := eye(n, n)
1362         x := make([]float64, n)
1363         v := make([]float64, n)
1364         for j := 0; j < n-1; j++ {
1365                 // x represents the j-th column of a random matrix.
1366                 for i := 0; i < j; i++ {
1367                         x[i] = 0
1368                 }
1369                 for i := j; i < n; i++ {
1370                         x[i] = rnd.NormFloat64()
1371                 }
1372                 // Compute v that represents the elementary reflector that
1373                 // annihilates the subdiagonal elements of x.
1374                 reflector(v, x, j)
1375                 // Compute Q * H_j and store the result into Q.
1376                 applyReflector(q, q, v)
1377         }
1378         if !isOrthonormal(q) {
1379                 panic("Q not orthogonal")
1380         }
1381         return q
1382 }
1383
1384 // reflector generates a Householder reflector v that zeros out subdiagonal
1385 // entries in the j-th column of a matrix.
1386 func reflector(v, col []float64, j int) {
1387         n := len(col)
1388         if len(v) != n {
1389                 panic("slice length mismatch")
1390         }
1391         if j < 0 || n <= j {
1392                 panic("invalid column index")
1393         }
1394
1395         for i := range v {
1396                 v[i] = 0
1397         }
1398         if j == n-1 {
1399                 return
1400         }
1401         s := floats.Norm(col[j:], 2)
1402         if s == 0 {
1403                 return
1404         }
1405         v[j] = col[j] + math.Copysign(s, col[j])
1406         copy(v[j+1:], col[j+1:])
1407         s = floats.Norm(v[j:], 2)
1408         floats.Scale(1/s, v[j:])
1409 }
1410
1411 // applyReflector computes Q*H where H is a Householder matrix represented by
1412 // the Householder reflector v.
1413 func applyReflector(qh blas64.General, q blas64.General, v []float64) {
1414         n := len(v)
1415         if qh.Rows != n || qh.Cols != n {
1416                 panic("bad size of qh")
1417         }
1418         if q.Rows != n || q.Cols != n {
1419                 panic("bad size of q")
1420         }
1421         qv := make([]float64, n)
1422         blas64.Gemv(blas.NoTrans, 1, q, blas64.Vector{1, v}, 0, blas64.Vector{1, qv})
1423         for i := 0; i < n; i++ {
1424                 for j := 0; j < n; j++ {
1425                         qh.Data[i*qh.Stride+j] = q.Data[i*q.Stride+j]
1426                 }
1427         }
1428         for i := 0; i < n; i++ {
1429                 for j := 0; j < n; j++ {
1430                         qh.Data[i*qh.Stride+j] -= 2 * qv[i] * v[j]
1431                 }
1432         }
1433         var norm2 float64
1434         for _, vi := range v {
1435                 norm2 += vi * vi
1436         }
1437         norm2inv := 1 / norm2
1438         for i := 0; i < n; i++ {
1439                 for j := 0; j < n; j++ {
1440                         qh.Data[i*qh.Stride+j] *= norm2inv
1441                 }
1442         }
1443 }
1444
1445 // constructGSVDresults returns the matrices [ 0 R ], D1 and D2 described
1446 // in the documentation of Dtgsja and Dggsvd3, and the result matrix in
1447 // the documentation for Dggsvp3.
1448 func constructGSVDresults(n, p, m, k, l int, a, b blas64.General, alpha, beta []float64) (zeroR, d1, d2 blas64.General) {
1449         // [ 0 R ]
1450         zeroR = zeros(k+l, n, n)
1451         dst := zeroR
1452         dst.Rows = min(m, k+l)
1453         dst.Cols = k + l
1454         dst.Data = zeroR.Data[n-k-l:]
1455         src := a
1456         src.Rows = min(m, k+l)
1457         src.Cols = k + l
1458         src.Data = a.Data[n-k-l:]
1459         copyGeneral(dst, src)
1460         if m < k+l {
1461                 // [ 0 R ]
1462                 dst.Rows = k + l - m
1463                 dst.Cols = k + l - m
1464                 dst.Data = zeroR.Data[m*zeroR.Stride+n-(k+l-m):]
1465                 src = b
1466                 src.Rows = k + l - m
1467                 src.Cols = k + l - m
1468                 src.Data = b.Data[(m-k)*b.Stride+n+m-k-l:]
1469                 copyGeneral(dst, src)
1470         }
1471
1472         // D1
1473         d1 = zeros(m, k+l, k+l)
1474         for i := 0; i < k; i++ {
1475                 d1.Data[i*d1.Stride+i] = 1
1476         }
1477         for i := k; i < min(m, k+l); i++ {
1478                 d1.Data[i*d1.Stride+i] = alpha[i]
1479         }
1480
1481         // D2
1482         d2 = zeros(p, k+l, k+l)
1483         for i := 0; i < min(l, m-k); i++ {
1484                 d2.Data[i*d2.Stride+i+k] = beta[k+i]
1485         }
1486         for i := m - k; i < l; i++ {
1487                 d2.Data[i*d2.Stride+i+k] = 1
1488         }
1489
1490         return zeroR, d1, d2
1491 }
1492
1493 func constructGSVPresults(n, p, m, k, l int, a, b blas64.General) (zeroA, zeroB blas64.General) {
1494         zeroA = zeros(m, n, n)
1495         dst := zeroA
1496         dst.Rows = min(m, k+l)
1497         dst.Cols = k + l
1498         dst.Data = zeroA.Data[n-k-l:]
1499         src := a
1500         dst.Rows = min(m, k+l)
1501         src.Cols = k + l
1502         src.Data = a.Data[n-k-l:]
1503         copyGeneral(dst, src)
1504
1505         zeroB = zeros(p, n, n)
1506         dst = zeroB
1507         dst.Rows = l
1508         dst.Cols = l
1509         dst.Data = zeroB.Data[n-l:]
1510         src = b
1511         dst.Rows = l
1512         src.Cols = l
1513         src.Data = b.Data[n-l:]
1514         copyGeneral(dst, src)
1515
1516         return zeroA, zeroB
1517 }