gonum.org/v1/gonum@v0.15.1-0.20240517103525-f853624cb1bb/lapack/gonum/dptcon.go (about)

     1  // Copyright ©2023 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  	"math"
     9  
    10  	"gonum.org/v1/gonum/blas/blas64"
    11  )
    12  
    13  // Dptcon computes and returns the reciprocal of the condition number (in the
    14  // 1-norm) of a symmetric positive definite tridiagonal matrix A using the
    15  // factorization A = L*D*Lᵀ or A = Uᵀ*D*U computed by Dpttrf.
    16  //
    17  // The reciprocal of the condition number is computed as
    18  //
    19  //	rcond = 1 / (anorm * ‖A⁻¹‖)
    20  //
    21  // and ‖A⁻¹‖ is computed by a direct method.
    22  //
    23  // d and e contain, respectively, the n diagonal elements of the diagonal matrix
    24  // D and the (n-1) off-diagonal elements of the unit bidiagonal factor U or L
    25  // from the factorization of A, as computed by Dpttrf.
    26  //
    27  // anorm is the 1-norm of the original matrix A.
    28  //
    29  // work must have length n, otherwise Dptcon will panic.
    30  func (impl Implementation) Dptcon(n int, d, e []float64, anorm float64, work []float64) (rcond float64) {
    31  	switch {
    32  	case n < 0:
    33  		panic(nLT0)
    34  	case anorm < 0:
    35  		panic(badNorm)
    36  	}
    37  
    38  	// Quick return if possible.
    39  	if n == 0 {
    40  		return 1
    41  	}
    42  
    43  	switch {
    44  	case len(d) < n:
    45  		panic(shortD)
    46  	case len(e) < n-1:
    47  		panic(shortE)
    48  	case len(work) < n:
    49  		panic(shortWork)
    50  	}
    51  
    52  	// Quick return if possible.
    53  	switch {
    54  	case anorm == 0:
    55  		return 0
    56  	case math.IsNaN(anorm):
    57  		// Propagate NaN.
    58  		return anorm
    59  	case math.IsInf(anorm, 1):
    60  		return 0
    61  	}
    62  
    63  	// Check that d[0:n] is positive.
    64  	for _, di := range d[:n] {
    65  		if di <= 0 {
    66  			return 0
    67  		}
    68  	}
    69  
    70  	// Solve M(A) * x = e, where M(A) = (m[i,j]) is given by
    71  	//
    72  	// 	m[i,j] =  abs(A[i,j]), i == j,
    73  	// 	m[i,j] = -abs(A[i,j]), i != j,
    74  	//
    75  	// and e = [1,1,...,1]ᵀ. Note M(A) = M(L)*D*M(L)ᵀ.
    76  
    77  	// Solve M(L) * b = e.
    78  	work[0] = 1
    79  	for i := 1; i < n; i++ {
    80  		work[i] = 1 + work[i-1]*math.Abs(e[i-1])
    81  	}
    82  
    83  	// Solve D * M(L)ᵀ * x = b.
    84  	work[n-1] /= d[n-1]
    85  	for i := n - 2; i >= 0; i-- {
    86  		work[i] = work[i]/d[i] + work[i+1]*math.Abs(e[i])
    87  	}
    88  
    89  	// Compute ainvnm = max(x[i]), 0<=i<n.
    90  	bi := blas64.Implementation()
    91  	ix := bi.Idamax(n, work, 1)
    92  	ainvnm := math.Abs(work[ix])
    93  	if ainvnm == 0 {
    94  		return 0
    95  	}
    96  
    97  	// Compute the reciprocal condition number.
    98  	return 1 / ainvnm / anorm
    99  }