gonum.org/v1/gonum@v0.14.0/lapack/gonum/dlauum.go (about) 1 // Copyright ©2018 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 gonum 6 7 import ( 8 "gonum.org/v1/gonum/blas" 9 "gonum.org/v1/gonum/blas/blas64" 10 ) 11 12 // Dlauum computes the product 13 // 14 // U * Uᵀ if uplo is blas.Upper 15 // Lᵀ * L if uplo is blas.Lower 16 // 17 // where U or L is stored in the upper or lower triangular part of A. 18 // Only the upper or lower triangle of the result is stored, overwriting 19 // the corresponding factor in A. 20 func (impl Implementation) Dlauum(uplo blas.Uplo, n int, a []float64, lda int) { 21 switch { 22 case uplo != blas.Upper && uplo != blas.Lower: 23 panic(badUplo) 24 case n < 0: 25 panic(nLT0) 26 case lda < max(1, n): 27 panic(badLdA) 28 } 29 30 // Quick return if possible. 31 if n == 0 { 32 return 33 } 34 35 if len(a) < (n-1)*lda+n { 36 panic(shortA) 37 } 38 39 // Determine the block size. 40 opts := "U" 41 if uplo == blas.Lower { 42 opts = "L" 43 } 44 nb := impl.Ilaenv(1, "DLAUUM", opts, n, -1, -1, -1) 45 46 if nb <= 1 || n <= nb { 47 // Use unblocked code. 48 impl.Dlauu2(uplo, n, a, lda) 49 return 50 } 51 52 // Use blocked code. 53 bi := blas64.Implementation() 54 if uplo == blas.Upper { 55 // Compute the product U*Uᵀ. 56 for i := 0; i < n; i += nb { 57 ib := min(nb, n-i) 58 bi.Dtrmm(blas.Right, blas.Upper, blas.Trans, blas.NonUnit, 59 i, ib, 1, a[i*lda+i:], lda, a[i:], lda) 60 impl.Dlauu2(blas.Upper, ib, a[i*lda+i:], lda) 61 if n-i-ib > 0 { 62 bi.Dgemm(blas.NoTrans, blas.Trans, i, ib, n-i-ib, 63 1, a[i+ib:], lda, a[i*lda+i+ib:], lda, 1, a[i:], lda) 64 bi.Dsyrk(blas.Upper, blas.NoTrans, ib, n-i-ib, 65 1, a[i*lda+i+ib:], lda, 1, a[i*lda+i:], lda) 66 } 67 } 68 } else { 69 // Compute the product Lᵀ*L. 70 for i := 0; i < n; i += nb { 71 ib := min(nb, n-i) 72 bi.Dtrmm(blas.Left, blas.Lower, blas.Trans, blas.NonUnit, 73 ib, i, 1, a[i*lda+i:], lda, a[i*lda:], lda) 74 impl.Dlauu2(blas.Lower, ib, a[i*lda+i:], lda) 75 if n-i-ib > 0 { 76 bi.Dgemm(blas.Trans, blas.NoTrans, ib, i, n-i-ib, 77 1, a[(i+ib)*lda+i:], lda, a[(i+ib)*lda:], lda, 1, a[i*lda:], lda) 78 bi.Dsyrk(blas.Lower, blas.Trans, ib, n-i-ib, 79 1, a[(i+ib)*lda+i:], lda, 1, a[i*lda+i:], lda) 80 } 81 } 82 } 83 }