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