vendor/gonum.org/v1/gonum/lapack/gonum/dtrtri.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 gonum
   6
   7 import (
   8         "gonum.org/v1/gonum/blas"
   9         "gonum.org/v1/gonum/blas/blas64"
  10 )
  11
  12 // Dtrtri computes the inverse of a triangular matrix, storing the result in place
  13 // into a. This is the BLAS level 3 version of the algorithm which builds upon
  14 // Dtrti2 to operate on matrix blocks instead of only individual columns.
  15 //
  16 // Dtrtri will not perform the inversion if the matrix is singular, and returns
  17 // a boolean indicating whether the inversion was successful.
  18 func (impl Implementation) Dtrtri(uplo blas.Uplo, diag blas.Diag, n int, a []float64, lda int) (ok bool) {
  19         checkMatrix(n, n, a, lda)
  20         if uplo != blas.Upper && uplo != blas.Lower {
  21                 panic(badUplo)
  22         }
  23         if diag != blas.NonUnit && diag != blas.Unit {
  24                 panic(badDiag)
  25         }
  26         if n == 0 {
  27                 return false
  28         }
  29         nonUnit := diag == blas.NonUnit
  30         if nonUnit {
  31                 for i := 0; i < n; i++ {
  32                         if a[i*lda+i] == 0 {
  33                                 return false
  34                         }
  35                 }
  36         }
  37
  38         bi := blas64.Implementation()
  39
  40         nb := impl.Ilaenv(1, "DTRTRI", "UD", n, -1, -1, -1)
  41         if nb <= 1 || nb > n {
  42                 impl.Dtrti2(uplo, diag, n, a, lda)
  43                 return true
  44         }
  45         if uplo == blas.Upper {
  46                 for j := 0; j < n; j += nb {
  47                         jb := min(nb, n-j)
  48                         bi.Dtrmm(blas.Left, blas.Upper, blas.NoTrans, diag, j, jb, 1, a, lda, a[j:], lda)
  49                         bi.Dtrsm(blas.Right, blas.Upper, blas.NoTrans, diag, j, jb, -1, a[j*lda+j:], lda, a[j:], lda)
  50                         impl.Dtrti2(blas.Upper, diag, jb, a[j*lda+j:], lda)
  51                 }
  52                 return true
  53         }
  54         nn := ((n - 1) / nb) * nb
  55         for j := nn; j >= 0; j -= nb {
  56                 jb := min(nb, n-j)
  57                 if j+jb <= n-1 {
  58                         bi.Dtrmm(blas.Left, blas.Lower, blas.NoTrans, diag, n-j-jb, jb, 1, a[(j+jb)*lda+j+jb:], lda, a[(j+jb)*lda+j:], lda)
  59                         bi.Dtrsm(blas.Right, blas.Lower, blas.NoTrans, diag, n-j-jb, jb, -1, a[j*lda+j:], lda, a[(j+jb)*lda+j:], lda)
  60                 }
  61                 impl.Dtrti2(blas.Lower, diag, jb, a[j*lda+j:], lda)
  62         }
  63         return true
  64 }