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