1// Copyright ©2013 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 mat 6 7import ( 8 "math" 9 10 "gonum.org/v1/gonum/blas" 11 "gonum.org/v1/gonum/blas/blas64" 12 "gonum.org/v1/gonum/lapack/lapack64" 13) 14 15const ( 16 badTriangle = "mat: invalid triangle" 17 badCholesky = "mat: invalid Cholesky factorization" 18) 19 20var ( 21 _ Matrix = (*Cholesky)(nil) 22 _ Symmetric = (*Cholesky)(nil) 23) 24 25// Cholesky is a symmetric positive definite matrix represented by its 26// Cholesky decomposition. 27// 28// The decomposition can be constructed using the Factorize method. The 29// factorization itself can be extracted using the UTo or LTo methods, and the 30// original symmetric matrix can be recovered with ToSym. 31// 32// Note that this matrix representation is useful for certain operations, in 33// particular finding solutions to linear equations. It is very inefficient 34// at other operations, in particular At is slow. 35// 36// Cholesky methods may only be called on a value that has been successfully 37// initialized by a call to Factorize that has returned true. Calls to methods 38// of an unsuccessful Cholesky factorization will panic. 39type Cholesky struct { 40 // The chol pointer must never be retained as a pointer outside the Cholesky 41 // struct, either by returning chol outside the struct or by setting it to 42 // a pointer coming from outside. The same prohibition applies to the data 43 // slice within chol. 44 chol *TriDense 45 cond float64 46} 47 48// updateCond updates the condition number of the Cholesky decomposition. If 49// norm > 0, then that norm is used as the norm of the original matrix A, otherwise 50// the norm is estimated from the decomposition. 51func (c *Cholesky) updateCond(norm float64) { 52 n := c.chol.mat.N 53 work := getFloats(3*n, false) 54 defer putFloats(work) 55 if norm < 0 { 56 // This is an approximation. By the definition of a norm, 57 // |AB| <= |A| |B|. 58 // Since A = Uᵀ*U, we get for the condition number κ that 59 // κ(A) := |A| |A^-1| = |Uᵀ*U| |A^-1| <= |Uᵀ| |U| |A^-1|, 60 // so this will overestimate the condition number somewhat. 61 // The norm of the original factorized matrix cannot be stored 62 // because of update possibilities. 63 unorm := lapack64.Lantr(CondNorm, c.chol.mat, work) 64 lnorm := lapack64.Lantr(CondNormTrans, c.chol.mat, work) 65 norm = unorm * lnorm 66 } 67 sym := c.chol.asSymBlas() 68 iwork := getInts(n, false) 69 v := lapack64.Pocon(sym, norm, work, iwork) 70 putInts(iwork) 71 c.cond = 1 / v 72} 73 74// Dims returns the dimensions of the matrix. 75func (ch *Cholesky) Dims() (r, c int) { 76 if !ch.valid() { 77 panic(badCholesky) 78 } 79 r, c = ch.chol.Dims() 80 return r, c 81} 82 83// At returns the element at row i, column j. 84func (c *Cholesky) At(i, j int) float64 { 85 if !c.valid() { 86 panic(badCholesky) 87 } 88 n := c.Symmetric() 89 if uint(i) >= uint(n) { 90 panic(ErrRowAccess) 91 } 92 if uint(j) >= uint(n) { 93 panic(ErrColAccess) 94 } 95 96 var val float64 97 for k := 0; k <= min(i, j); k++ { 98 val += c.chol.at(k, i) * c.chol.at(k, j) 99 } 100 return val 101} 102 103// T returns the the receiver, the transpose of a symmetric matrix. 104func (c *Cholesky) T() Matrix { 105 return c 106} 107 108// Symmetric implements the Symmetric interface and returns the number of rows 109// in the matrix (this is also the number of columns). 110func (c *Cholesky) Symmetric() int { 111 r, _ := c.chol.Dims() 112 return r 113} 114 115// Cond returns the condition number of the factorized matrix. 116func (c *Cholesky) Cond() float64 { 117 if !c.valid() { 118 panic(badCholesky) 119 } 120 return c.cond 121} 122 123// Factorize calculates the Cholesky decomposition of the matrix A and returns 124// whether the matrix is positive definite. If Factorize returns false, the 125// factorization must not be used. 126func (c *Cholesky) Factorize(a Symmetric) (ok bool) { 127 n := a.Symmetric() 128 if c.chol == nil { 129 c.chol = NewTriDense(n, Upper, nil) 130 } else { 131 c.chol = NewTriDense(n, Upper, use(c.chol.mat.Data, n*n)) 132 } 133 copySymIntoTriangle(c.chol, a) 134 135 sym := c.chol.asSymBlas() 136 work := getFloats(c.chol.mat.N, false) 137 norm := lapack64.Lansy(CondNorm, sym, work) 138 putFloats(work) 139 _, ok = lapack64.Potrf(sym) 140 if ok { 141 c.updateCond(norm) 142 } else { 143 c.Reset() 144 } 145 return ok 146} 147 148// Reset resets the factorization so that it can be reused as the receiver of a 149// dimensionally restricted operation. 150func (c *Cholesky) Reset() { 151 if c.chol != nil { 152 c.chol.Reset() 153 } 154 c.cond = math.Inf(1) 155} 156 157// IsEmpty returns whether the receiver is empty. Empty matrices can be the 158// receiver for size-restricted operations. The receiver can be emptied using 159// Reset. 160func (c *Cholesky) IsEmpty() bool { 161 return c.chol == nil || c.chol.IsEmpty() 162} 163 164// SetFromU sets the Cholesky decomposition from the given triangular matrix. 165// SetFromU panics if t is not upper triangular. If the receiver is empty it 166// is resized to be n×n, the size of t. If dst is non-empty, SetFromU panics 167// if c is not of size n×n. Note that t is copied into, not stored inside, the 168// receiver. 169func (c *Cholesky) SetFromU(t Triangular) { 170 n, kind := t.Triangle() 171 if kind != Upper { 172 panic("cholesky: matrix must be upper triangular") 173 } 174 if c.chol == nil { 175 c.chol = NewTriDense(n, Upper, nil) 176 } else { 177 c.chol.reuseAsNonZeroed(n, Upper) 178 } 179 c.chol.Copy(t) 180 c.updateCond(-1) 181} 182 183// Clone makes a copy of the input Cholesky into the receiver, overwriting the 184// previous value of the receiver. Clone does not place any restrictions on receiver 185// shape. Clone panics if the input Cholesky is not the result of a valid decomposition. 186func (c *Cholesky) Clone(chol *Cholesky) { 187 if !chol.valid() { 188 panic(badCholesky) 189 } 190 n := chol.Symmetric() 191 if c.chol == nil { 192 c.chol = NewTriDense(n, Upper, nil) 193 } else { 194 c.chol = NewTriDense(n, Upper, use(c.chol.mat.Data, n*n)) 195 } 196 c.chol.Copy(chol.chol) 197 c.cond = chol.cond 198} 199 200// Det returns the determinant of the matrix that has been factorized. 201func (c *Cholesky) Det() float64 { 202 if !c.valid() { 203 panic(badCholesky) 204 } 205 return math.Exp(c.LogDet()) 206} 207 208// LogDet returns the log of the determinant of the matrix that has been factorized. 209func (c *Cholesky) LogDet() float64 { 210 if !c.valid() { 211 panic(badCholesky) 212 } 213 var det float64 214 for i := 0; i < c.chol.mat.N; i++ { 215 det += 2 * math.Log(c.chol.mat.Data[i*c.chol.mat.Stride+i]) 216 } 217 return det 218} 219 220// SolveTo finds the matrix X that solves A * X = B where A is represented 221// by the Cholesky decomposition. The result is stored in-place into dst. 222func (c *Cholesky) SolveTo(dst *Dense, b Matrix) error { 223 if !c.valid() { 224 panic(badCholesky) 225 } 226 n := c.chol.mat.N 227 bm, bn := b.Dims() 228 if n != bm { 229 panic(ErrShape) 230 } 231 232 dst.reuseAsNonZeroed(bm, bn) 233 if b != dst { 234 dst.Copy(b) 235 } 236 lapack64.Potrs(c.chol.mat, dst.mat) 237 if c.cond > ConditionTolerance { 238 return Condition(c.cond) 239 } 240 return nil 241} 242 243// SolveCholTo finds the matrix X that solves A * X = B where A and B are represented 244// by their Cholesky decompositions a and b. The result is stored in-place into 245// dst. 246func (a *Cholesky) SolveCholTo(dst *Dense, b *Cholesky) error { 247 if !a.valid() || !b.valid() { 248 panic(badCholesky) 249 } 250 bn := b.chol.mat.N 251 if a.chol.mat.N != bn { 252 panic(ErrShape) 253 } 254 255 dst.reuseAsZeroed(bn, bn) 256 dst.Copy(b.chol.T()) 257 blas64.Trsm(blas.Left, blas.Trans, 1, a.chol.mat, dst.mat) 258 blas64.Trsm(blas.Left, blas.NoTrans, 1, a.chol.mat, dst.mat) 259 blas64.Trmm(blas.Right, blas.NoTrans, 1, b.chol.mat, dst.mat) 260 if a.cond > ConditionTolerance { 261 return Condition(a.cond) 262 } 263 return nil 264} 265 266// SolveVecTo finds the vector X that solves A * x = b where A is represented 267// by the Cholesky decomposition. The result is stored in-place into 268// dst. 269func (c *Cholesky) SolveVecTo(dst *VecDense, b Vector) error { 270 if !c.valid() { 271 panic(badCholesky) 272 } 273 n := c.chol.mat.N 274 if br, bc := b.Dims(); br != n || bc != 1 { 275 panic(ErrShape) 276 } 277 switch rv := b.(type) { 278 default: 279 dst.reuseAsNonZeroed(n) 280 return c.SolveTo(dst.asDense(), b) 281 case RawVectorer: 282 bmat := rv.RawVector() 283 if dst != b { 284 dst.checkOverlap(bmat) 285 } 286 dst.reuseAsNonZeroed(n) 287 if dst != b { 288 dst.CopyVec(b) 289 } 290 lapack64.Potrs(c.chol.mat, dst.asGeneral()) 291 if c.cond > ConditionTolerance { 292 return Condition(c.cond) 293 } 294 return nil 295 } 296} 297 298// RawU returns the Triangular matrix used to store the Cholesky decomposition of 299// the original matrix A. The returned matrix should not be modified. If it is 300// modified, the decomposition is invalid and should not be used. 301func (c *Cholesky) RawU() Triangular { 302 return c.chol 303} 304 305// UTo stores into dst the n×n upper triangular matrix U from a Cholesky 306// decomposition 307// A = Uᵀ * U. 308// If dst is empty, it is resized to be an n×n upper triangular matrix. When dst 309// is non-empty, UTo panics if dst is not n×n or not Upper. UTo will also panic 310// if the receiver does not contain a successful factorization. 311func (c *Cholesky) UTo(dst *TriDense) { 312 if !c.valid() { 313 panic(badCholesky) 314 } 315 n := c.chol.mat.N 316 if dst.IsEmpty() { 317 dst.ReuseAsTri(n, Upper) 318 } else { 319 n2, kind := dst.Triangle() 320 if n != n2 { 321 panic(ErrShape) 322 } 323 if kind != Upper { 324 panic(ErrTriangle) 325 } 326 } 327 dst.Copy(c.chol) 328} 329 330// LTo stores into dst the n×n lower triangular matrix L from a Cholesky 331// decomposition 332// A = L * Lᵀ. 333// If dst is empty, it is resized to be an n×n lower triangular matrix. When dst 334// is non-empty, LTo panics if dst is not n×n or not Lower. LTo will also panic 335// if the receiver does not contain a successful factorization. 336func (c *Cholesky) LTo(dst *TriDense) { 337 if !c.valid() { 338 panic(badCholesky) 339 } 340 n := c.chol.mat.N 341 if dst.IsEmpty() { 342 dst.ReuseAsTri(n, Lower) 343 } else { 344 n2, kind := dst.Triangle() 345 if n != n2 { 346 panic(ErrShape) 347 } 348 if kind != Lower { 349 panic(ErrTriangle) 350 } 351 } 352 dst.Copy(c.chol.TTri()) 353} 354 355// ToSym reconstructs the original positive definite matrix from its 356// Cholesky decomposition, storing the result into dst. If dst is 357// empty it is resized to be n×n. If dst is non-empty, ToSym panics 358// if dst is not of size n×n. ToSym will also panic if the receiver 359// does not contain a successful factorization. 360func (c *Cholesky) ToSym(dst *SymDense) { 361 if !c.valid() { 362 panic(badCholesky) 363 } 364 n := c.chol.mat.N 365 if dst.IsEmpty() { 366 dst.ReuseAsSym(n) 367 } else { 368 n2 := dst.Symmetric() 369 if n != n2 { 370 panic(ErrShape) 371 } 372 } 373 // Create a TriDense representing the Cholesky factor U with dst's 374 // backing slice. 375 // Operations on u are reflected in s. 376 u := &TriDense{ 377 mat: blas64.Triangular{ 378 Uplo: blas.Upper, 379 Diag: blas.NonUnit, 380 N: n, 381 Data: dst.mat.Data, 382 Stride: dst.mat.Stride, 383 }, 384 cap: n, 385 } 386 u.Copy(c.chol) 387 // Compute the product Uᵀ*U using the algorithm from LAPACK/TESTING/LIN/dpot01.f 388 a := u.mat.Data 389 lda := u.mat.Stride 390 bi := blas64.Implementation() 391 for k := n - 1; k >= 0; k-- { 392 a[k*lda+k] = bi.Ddot(k+1, a[k:], lda, a[k:], lda) 393 if k > 0 { 394 bi.Dtrmv(blas.Upper, blas.Trans, blas.NonUnit, k, a, lda, a[k:], lda) 395 } 396 } 397} 398 399// InverseTo computes the inverse of the matrix represented by its Cholesky 400// factorization and stores the result into s. If the factorized 401// matrix is ill-conditioned, a Condition error will be returned. 402// Note that matrix inversion is numerically unstable, and should generally be 403// avoided where possible, for example by using the Solve routines. 404func (c *Cholesky) InverseTo(s *SymDense) error { 405 if !c.valid() { 406 panic(badCholesky) 407 } 408 s.reuseAsNonZeroed(c.chol.mat.N) 409 // Create a TriDense representing the Cholesky factor U with the backing 410 // slice from s. 411 // Operations on u are reflected in s. 412 u := &TriDense{ 413 mat: blas64.Triangular{ 414 Uplo: blas.Upper, 415 Diag: blas.NonUnit, 416 N: s.mat.N, 417 Data: s.mat.Data, 418 Stride: s.mat.Stride, 419 }, 420 cap: s.mat.N, 421 } 422 u.Copy(c.chol) 423 424 _, ok := lapack64.Potri(u.mat) 425 if !ok { 426 return Condition(math.Inf(1)) 427 } 428 if c.cond > ConditionTolerance { 429 return Condition(c.cond) 430 } 431 return nil 432} 433 434// Scale multiplies the original matrix A by a positive constant using 435// its Cholesky decomposition, storing the result in-place into the receiver. 436// That is, if the original Cholesky factorization is 437// Uᵀ * U = A 438// the updated factorization is 439// U'ᵀ * U' = f A = A' 440// Scale panics if the constant is non-positive, or if the receiver is non-empty 441// and is of a different size from the input. 442func (c *Cholesky) Scale(f float64, orig *Cholesky) { 443 if !orig.valid() { 444 panic(badCholesky) 445 } 446 if f <= 0 { 447 panic("cholesky: scaling by a non-positive constant") 448 } 449 n := orig.Symmetric() 450 if c.chol == nil { 451 c.chol = NewTriDense(n, Upper, nil) 452 } else if c.chol.mat.N != n { 453 panic(ErrShape) 454 } 455 c.chol.ScaleTri(math.Sqrt(f), orig.chol) 456 c.cond = orig.cond // Scaling by a positive constant does not change the condition number. 457} 458 459// ExtendVecSym computes the Cholesky decomposition of the original matrix A, 460// whose Cholesky decomposition is in a, extended by a the n×1 vector v according to 461// [A w] 462// [w' k] 463// where k = v[n-1] and w = v[:n-1]. The result is stored into the receiver. 464// In order for the updated matrix to be positive definite, it must be the case 465// that k > w' A^-1 w. If this condition does not hold then ExtendVecSym will 466// return false and the receiver will not be updated. 467// 468// ExtendVecSym will panic if v.Len() != a.Symmetric()+1 or if a does not contain 469// a valid decomposition. 470func (c *Cholesky) ExtendVecSym(a *Cholesky, v Vector) (ok bool) { 471 n := a.Symmetric() 472 473 if v.Len() != n+1 { 474 panic(badSliceLength) 475 } 476 if !a.valid() { 477 panic(badCholesky) 478 } 479 480 // The algorithm is commented here, but see also 481 // https://math.stackexchange.com/questions/955874/cholesky-factor-when-adding-a-row-and-column-to-already-factorized-matrix 482 // We have A and want to compute the Cholesky of 483 // [A w] 484 // [w' k] 485 // We want 486 // [U c] 487 // [0 d] 488 // to be the updated Cholesky, and so it must be that 489 // [A w] = [U' 0] [U c] 490 // [w' k] [c' d] [0 d] 491 // Thus, we need 492 // 1) A = U'U (true by the original decomposition being valid), 493 // 2) U' * c = w => c = U'^-1 w 494 // 3) c'*c + d'*d = k => d = sqrt(k-c'*c) 495 496 // First, compute c = U'^-1 a 497 // TODO(btracey): Replace this with CopyVec when issue 167 is fixed. 498 w := NewVecDense(n, nil) 499 for i := 0; i < n; i++ { 500 w.SetVec(i, v.At(i, 0)) 501 } 502 k := v.At(n, 0) 503 504 var t VecDense 505 t.SolveVec(a.chol.T(), w) 506 507 dot := Dot(&t, &t) 508 if dot >= k { 509 return false 510 } 511 d := math.Sqrt(k - dot) 512 513 newU := NewTriDense(n+1, Upper, nil) 514 newU.Copy(a.chol) 515 for i := 0; i < n; i++ { 516 newU.SetTri(i, n, t.At(i, 0)) 517 } 518 newU.SetTri(n, n, d) 519 c.chol = newU 520 c.updateCond(-1) 521 return true 522} 523 524// SymRankOne performs a rank-1 update of the original matrix A and refactorizes 525// its Cholesky factorization, storing the result into the receiver. That is, if 526// in the original Cholesky factorization 527// Uᵀ * U = A, 528// in the updated factorization 529// U'ᵀ * U' = A + alpha * x * xᵀ = A'. 530// 531// Note that when alpha is negative, the updating problem may be ill-conditioned 532// and the results may be inaccurate, or the updated matrix A' may not be 533// positive definite and not have a Cholesky factorization. SymRankOne returns 534// whether the updated matrix A' is positive definite. 535// 536// SymRankOne updates a Cholesky factorization in O(n²) time. The Cholesky 537// factorization computation from scratch is O(n³). 538func (c *Cholesky) SymRankOne(orig *Cholesky, alpha float64, x Vector) (ok bool) { 539 if !orig.valid() { 540 panic(badCholesky) 541 } 542 n := orig.Symmetric() 543 if r, c := x.Dims(); r != n || c != 1 { 544 panic(ErrShape) 545 } 546 if orig != c { 547 if c.chol == nil { 548 c.chol = NewTriDense(n, Upper, nil) 549 } else if c.chol.mat.N != n { 550 panic(ErrShape) 551 } 552 c.chol.Copy(orig.chol) 553 } 554 555 if alpha == 0 { 556 return true 557 } 558 559 // Algorithms for updating and downdating the Cholesky factorization are 560 // described, for example, in 561 // - J. J. Dongarra, J. R. Bunch, C. B. Moler, G. W. Stewart: LINPACK 562 // Users' Guide. SIAM (1979), pages 10.10--10.14 563 // or 564 // - P. E. Gill, G. H. Golub, W. Murray, and M. A. Saunders: Methods for 565 // modifying matrix factorizations. Mathematics of Computation 28(126) 566 // (1974), Method C3 on page 521 567 // 568 // The implementation is based on LINPACK code 569 // http://www.netlib.org/linpack/dchud.f 570 // http://www.netlib.org/linpack/dchdd.f 571 // and 572 // https://icl.cs.utk.edu/lapack-forum/viewtopic.php?f=2&t=2646 573 // 574 // According to http://icl.cs.utk.edu/lapack-forum/archives/lapack/msg00301.html 575 // LINPACK is released under BSD license. 576 // 577 // See also: 578 // - M. A. Saunders: Large-scale Linear Programming Using the Cholesky 579 // Factorization. Technical Report Stanford University (1972) 580 // http://i.stanford.edu/pub/cstr/reports/cs/tr/72/252/CS-TR-72-252.pdf 581 // - Matthias Seeger: Low rank updates for the Cholesky decomposition. 582 // EPFL Technical Report 161468 (2004) 583 // http://infoscience.epfl.ch/record/161468 584 585 work := getFloats(n, false) 586 defer putFloats(work) 587 var xmat blas64.Vector 588 if rv, ok := x.(RawVectorer); ok { 589 xmat = rv.RawVector() 590 } else { 591 var tmp *VecDense 592 tmp.CopyVec(x) 593 xmat = tmp.RawVector() 594 } 595 blas64.Copy(xmat, blas64.Vector{N: n, Data: work, Inc: 1}) 596 597 if alpha > 0 { 598 // Compute rank-1 update. 599 if alpha != 1 { 600 blas64.Scal(math.Sqrt(alpha), blas64.Vector{N: n, Data: work, Inc: 1}) 601 } 602 umat := c.chol.mat 603 stride := umat.Stride 604 for i := 0; i < n; i++ { 605 // Compute parameters of the Givens matrix that zeroes 606 // the i-th element of x. 607 c, s, r, _ := blas64.Rotg(umat.Data[i*stride+i], work[i]) 608 if r < 0 { 609 // Multiply by -1 to have positive diagonal 610 // elemnts. 611 r *= -1 612 c *= -1 613 s *= -1 614 } 615 umat.Data[i*stride+i] = r 616 if i < n-1 { 617 // Multiply the extended factorization matrix by 618 // the Givens matrix from the left. Only 619 // the i-th row and x are modified. 620 blas64.Rot( 621 blas64.Vector{N: n - i - 1, Data: umat.Data[i*stride+i+1 : i*stride+n], Inc: 1}, 622 blas64.Vector{N: n - i - 1, Data: work[i+1 : n], Inc: 1}, 623 c, s) 624 } 625 } 626 c.updateCond(-1) 627 return true 628 } 629 630 // Compute rank-1 downdate. 631 alpha = math.Sqrt(-alpha) 632 if alpha != 1 { 633 blas64.Scal(alpha, blas64.Vector{N: n, Data: work, Inc: 1}) 634 } 635 // Solve Uᵀ * p = x storing the result into work. 636 ok = lapack64.Trtrs(blas.Trans, c.chol.RawTriangular(), blas64.General{ 637 Rows: n, 638 Cols: 1, 639 Stride: 1, 640 Data: work, 641 }) 642 if !ok { 643 // The original matrix is singular. Should not happen, because 644 // the factorization is valid. 645 panic(badCholesky) 646 } 647 norm := blas64.Nrm2(blas64.Vector{N: n, Data: work, Inc: 1}) 648 if norm >= 1 { 649 // The updated matrix is not positive definite. 650 return false 651 } 652 norm = math.Sqrt((1 + norm) * (1 - norm)) 653 cos := getFloats(n, false) 654 defer putFloats(cos) 655 sin := getFloats(n, false) 656 defer putFloats(sin) 657 for i := n - 1; i >= 0; i-- { 658 // Compute parameters of Givens matrices that zero elements of p 659 // backwards. 660 cos[i], sin[i], norm, _ = blas64.Rotg(norm, work[i]) 661 if norm < 0 { 662 norm *= -1 663 cos[i] *= -1 664 sin[i] *= -1 665 } 666 } 667 umat := c.chol.mat 668 stride := umat.Stride 669 for i := n - 1; i >= 0; i-- { 670 work[i] = 0 671 // Apply Givens matrices to U. 672 // TODO(vladimir-ch): Use workspace to avoid modifying the 673 // receiver in case an invalid factorization is created. 674 blas64.Rot( 675 blas64.Vector{N: n - i, Data: work[i:n], Inc: 1}, 676 blas64.Vector{N: n - i, Data: umat.Data[i*stride+i : i*stride+n], Inc: 1}, 677 cos[i], sin[i]) 678 if umat.Data[i*stride+i] == 0 { 679 // The matrix is singular (may rarely happen due to 680 // floating-point effects?). 681 ok = false 682 } else if umat.Data[i*stride+i] < 0 { 683 // Diagonal elements should be positive. If it happens 684 // that on the i-th row the diagonal is negative, 685 // multiply U from the left by an identity matrix that 686 // has -1 on the i-th row. 687 blas64.Scal(-1, blas64.Vector{N: n - i, Data: umat.Data[i*stride+i : i*stride+n], Inc: 1}) 688 } 689 } 690 if ok { 691 c.updateCond(-1) 692 } else { 693 c.Reset() 694 } 695 return ok 696} 697 698func (c *Cholesky) valid() bool { 699 return c.chol != nil && !c.chol.IsEmpty() 700} 701