gonum.org/v1/gonum@v0.15.1-0.20240517103525-f853624cb1bb/lapack/gonum/dgetc2.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  // Dgetc2 computes an LU factorization with complete pivoting of the n×n matrix
    14  // A. The factorization has the form
    15  //
    16  //	A = P * L * U * Q,
    17  //
    18  // where P and Q are permutation matrices, L is lower triangular with unit
    19  // diagonal elements and U is upper triangular.
    20  //
    21  // On entry, a contains the matrix A to be factored. On return, a is overwritten
    22  // with the factors L and U. The unit diagonal elements of L are not stored.
    23  //
    24  // On return, ipiv and jpiv contain the pivot indices: row i has been
    25  // interchanged with row ipiv[i] and column j has been interchanged with column
    26  // jpiv[j]. ipiv and jpiv must have length n, otherwise Dgetc2 will panic.
    27  //
    28  // If k is non-negative, then U[k,k] is likely to produce overflow when solving
    29  // for x in A*x=b and U has been perturbed to avoid the overflow.
    30  //
    31  // Dgetc2 is an internal routine. It is exported for testing purposes.
    32  func (impl Implementation) Dgetc2(n int, a []float64, lda int, ipiv, jpiv []int) (k int) {
    33  	switch {
    34  	case n < 0:
    35  		panic(nLT0)
    36  	case lda < max(1, n):
    37  		panic(badLdA)
    38  	}
    39  
    40  	// Negative k indicates U was not perturbed.
    41  	k = -1
    42  
    43  	// Quick return if possible.
    44  	if n == 0 {
    45  		return k
    46  	}
    47  
    48  	switch {
    49  	case len(a) < (n-1)*lda+n:
    50  		panic(shortA)
    51  	case len(ipiv) != n:
    52  		panic(badLenIpiv)
    53  	case len(jpiv) != n:
    54  		panic(badLenJpvt)
    55  	}
    56  
    57  	const (
    58  		eps    = dlamchP
    59  		smlnum = dlamchS / eps
    60  	)
    61  
    62  	if n == 1 {
    63  		ipiv[0], jpiv[0] = 0, 0
    64  		if math.Abs(a[0]) < smlnum {
    65  			k = 0
    66  			a[0] = smlnum
    67  		}
    68  		return k
    69  	}
    70  
    71  	// Factorize A using complete pivoting.
    72  	// Set pivots less than smin to smin.
    73  	var smin float64
    74  	var ipv, jpv int
    75  	bi := blas64.Implementation()
    76  	for i := 0; i < n-1; i++ {
    77  		var xmax float64
    78  		for ip := i; ip < n; ip++ {
    79  			for jp := i; jp < n; jp++ {
    80  				if math.Abs(a[ip*lda+jp]) >= xmax {
    81  					xmax = math.Abs(a[ip*lda+jp])
    82  					ipv = ip
    83  					jpv = jp
    84  				}
    85  			}
    86  		}
    87  		if i == 0 {
    88  			smin = math.Max(eps*xmax, smlnum)
    89  		}
    90  
    91  		// Swap rows.
    92  		if ipv != i {
    93  			bi.Dswap(n, a[ipv*lda:], 1, a[i*lda:], 1)
    94  		}
    95  		ipiv[i] = ipv
    96  
    97  		// Swap columns.
    98  		if jpv != i {
    99  			bi.Dswap(n, a[jpv:], lda, a[i:], lda)
   100  		}
   101  		jpiv[i] = jpv
   102  
   103  		// Check for singularity.
   104  		if math.Abs(a[i*lda+i]) < smin {
   105  			k = i
   106  			a[i*lda+i] = smin
   107  		}
   108  
   109  		for j := i + 1; j < n; j++ {
   110  			a[j*lda+i] /= a[i*lda+i]
   111  		}
   112  		bi.Dger(n-i-1, n-i-1, -1, a[(i+1)*lda+i:], lda, a[i*lda+i+1:], 1, a[(i+1)*lda+i+1:], lda)
   113  	}
   114  
   115  	if math.Abs(a[(n-1)*lda+n-1]) < smin {
   116  		k = n - 1
   117  		a[(n-1)*lda+(n-1)] = smin
   118  	}
   119  
   120  	// Set last pivots to last index.
   121  	ipiv[n-1] = n - 1
   122  	jpiv[n-1] = n - 1
   123  
   124  	return k
   125  }