gonum.org/v1/gonum@v0.14.0/lapack/gonum/dgesc2.go (about)

     1  // Copyright ©2021 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  // Dgesc2 solves a system of linear equations
    14  //
    15  //	A * x = scale * b
    16  //
    17  // with a general n×n matrix A represented by the LU factorization with complete
    18  // pivoting
    19  //
    20  //	A = P * L * U * Q
    21  //
    22  // as computed by Dgetc2.
    23  //
    24  // On entry, rhs contains the right hand side vector b. On return, it is
    25  // overwritten with the solution vector x.
    26  //
    27  // Dgesc2 returns a scale factor
    28  //
    29  //	0 <= scale <= 1
    30  //
    31  // chosen to prevent overflow in the solution.
    32  //
    33  // Dgesc2 is an internal routine. It is exported for testing purposes.
    34  func (impl Implementation) Dgesc2(n int, a []float64, lda int, rhs []float64, ipiv, jpiv []int) (scale float64) {
    35  	switch {
    36  	case n < 0:
    37  		panic(nLT0)
    38  	case lda < max(1, n):
    39  		panic(badLdA)
    40  	}
    41  
    42  	// Quick return if possible.
    43  	if n == 0 {
    44  		return 0
    45  	}
    46  
    47  	switch {
    48  	case len(a) < (n-1)*lda+n:
    49  		panic(shortA)
    50  	case len(rhs) < n:
    51  		panic(shortRHS)
    52  	case len(ipiv) != n:
    53  		panic(badLenIpiv)
    54  	case len(jpiv) != n:
    55  		panic(badLenJpiv)
    56  	}
    57  
    58  	const smlnum = dlamchS / dlamchP
    59  
    60  	// Apply permutations ipiv to rhs.
    61  	impl.Dlaswp(1, rhs, 1, 0, n-1, ipiv[:n], 1)
    62  
    63  	// Solve for L part.
    64  	for i := 0; i < n-1; i++ {
    65  		for j := i + 1; j < n; j++ {
    66  			rhs[j] -= float64(a[j*lda+i] * rhs[i])
    67  		}
    68  	}
    69  
    70  	// Check for scaling.
    71  	scale = 1.0
    72  	bi := blas64.Implementation()
    73  	i := bi.Idamax(n, rhs, 1)
    74  	if 2*smlnum*math.Abs(rhs[i]) > math.Abs(a[(n-1)*lda+(n-1)]) {
    75  		temp := 0.5 / math.Abs(rhs[i])
    76  		bi.Dscal(n, temp, rhs, 1)
    77  		scale *= temp
    78  	}
    79  
    80  	// Solve for U part.
    81  	for i := n - 1; i >= 0; i-- {
    82  		temp := 1.0 / a[i*lda+i]
    83  		rhs[i] *= temp
    84  		for j := i + 1; j < n; j++ {
    85  			rhs[i] -= float64(rhs[j] * (a[i*lda+j] * temp))
    86  		}
    87  	}
    88  
    89  	// Apply permutations jpiv to the solution (rhs).
    90  	impl.Dlaswp(1, rhs, 1, 0, n-1, jpiv[:n], -1)
    91  
    92  	return scale
    93  }