github.com/jingcheng-WU/gonum@v0.9.1-0.20210323123734-f1a2a11a8f7b/lapack/gonum/dlaqp2.go (about)

     1  // Copyright ©2017 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/jingcheng-WU/gonum/blas"
    11  	"github.com/jingcheng-WU/gonum/blas/blas64"
    12  )
    13  
    14  // Dlaqp2 computes a QR factorization with column pivoting of the block A[offset:m, 0:n]
    15  // of the m×n matrix A. The block A[0:offset, 0:n] is accordingly pivoted, but not factorized.
    16  //
    17  // On exit, the upper triangle of block A[offset:m, 0:n] is the triangular factor obtained.
    18  // The elements in block A[offset:m, 0:n] below the diagonal, together with tau, represent
    19  // the orthogonal matrix Q as a product of elementary reflectors.
    20  //
    21  // offset is number of rows of the matrix A that must be pivoted but not factorized.
    22  // offset must not be negative otherwise Dlaqp2 will panic.
    23  //
    24  // On exit, jpvt holds the permutation that was applied; the jth column of A*P was the
    25  // jpvt[j] column of A. jpvt must have length n, otherwise Dlaqp2 will panic.
    26  //
    27  // On exit tau holds the scalar factors of the elementary reflectors. It must have length
    28  // at least min(m-offset, n) otherwise Dlaqp2 will panic.
    29  //
    30  // vn1 and vn2 hold the partial and complete column norms respectively. They must have length n,
    31  // otherwise Dlaqp2 will panic.
    32  //
    33  // work must have length n, otherwise Dlaqp2 will panic.
    34  //
    35  // Dlaqp2 is an internal routine. It is exported for testing purposes.
    36  func (impl Implementation) Dlaqp2(m, n, offset int, a []float64, lda int, jpvt []int, tau, vn1, vn2, work []float64) {
    37  	switch {
    38  	case m < 0:
    39  		panic(mLT0)
    40  	case n < 0:
    41  		panic(nLT0)
    42  	case offset < 0:
    43  		panic(offsetLT0)
    44  	case offset > m:
    45  		panic(offsetGTM)
    46  	case lda < max(1, n):
    47  		panic(badLdA)
    48  	}
    49  
    50  	// Quick return if possible.
    51  	if m == 0 || n == 0 {
    52  		return
    53  	}
    54  
    55  	mn := min(m-offset, n)
    56  	switch {
    57  	case len(a) < (m-1)*lda+n:
    58  		panic(shortA)
    59  	case len(jpvt) != n:
    60  		panic(badLenJpvt)
    61  	case len(tau) < mn:
    62  		panic(shortTau)
    63  	case len(vn1) < n:
    64  		panic(shortVn1)
    65  	case len(vn2) < n:
    66  		panic(shortVn2)
    67  	case len(work) < n:
    68  		panic(shortWork)
    69  	}
    70  
    71  	tol3z := math.Sqrt(dlamchE)
    72  
    73  	bi := blas64.Implementation()
    74  
    75  	// Compute factorization.
    76  	for i := 0; i < mn; i++ {
    77  		offpi := offset + i
    78  
    79  		// Determine ith pivot column and swap if necessary.
    80  		p := i + bi.Idamax(n-i, vn1[i:], 1)
    81  		if p != i {
    82  			bi.Dswap(m, a[p:], lda, a[i:], lda)
    83  			jpvt[p], jpvt[i] = jpvt[i], jpvt[p]
    84  			vn1[p] = vn1[i]
    85  			vn2[p] = vn2[i]
    86  		}
    87  
    88  		// Generate elementary reflector H_i.
    89  		if offpi < m-1 {
    90  			a[offpi*lda+i], tau[i] = impl.Dlarfg(m-offpi, a[offpi*lda+i], a[(offpi+1)*lda+i:], lda)
    91  		} else {
    92  			tau[i] = 0
    93  		}
    94  
    95  		if i < n-1 {
    96  			// Apply H_iᵀ to A[offset+i:m, i:n] from the left.
    97  			aii := a[offpi*lda+i]
    98  			a[offpi*lda+i] = 1
    99  			impl.Dlarf(blas.Left, m-offpi, n-i-1, a[offpi*lda+i:], lda, tau[i], a[offpi*lda+i+1:], lda, work)
   100  			a[offpi*lda+i] = aii
   101  		}
   102  
   103  		// Update partial column norms.
   104  		for j := i + 1; j < n; j++ {
   105  			if vn1[j] == 0 {
   106  				continue
   107  			}
   108  
   109  			// The following marked lines follow from the
   110  			// analysis in Lapack Working Note 176.
   111  			r := math.Abs(a[offpi*lda+j]) / vn1[j] // *
   112  			temp := math.Max(0, 1-r*r)             // *
   113  			r = vn1[j] / vn2[j]                    // *
   114  			temp2 := temp * r * r                  // *
   115  			if temp2 < tol3z {
   116  				var v float64
   117  				if offpi < m-1 {
   118  					v = bi.Dnrm2(m-offpi-1, a[(offpi+1)*lda+j:], lda)
   119  				}
   120  				vn1[j] = v
   121  				vn2[j] = v
   122  			} else {
   123  				vn1[j] *= math.Sqrt(temp) // *
   124  			}
   125  		}
   126  	}
   127  }