github.com/gonum/lapack@v0.0.0-20181123203213-e4cdc5a0bff9/native/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 native
     6  
     7  import (
     8  	"math"
     9  
    10  	"github.com/gonum/blas"
    11  	"github.com/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  	checkMatrix(m, n, a, lda)
    38  	if len(jpvt) != n {
    39  		panic(badIpiv)
    40  	}
    41  	mn := min(m-offset, n)
    42  	if len(tau) < mn {
    43  		panic(badTau)
    44  	}
    45  	if len(vn1) < n {
    46  		panic(badVn1)
    47  	}
    48  	if len(vn2) < n {
    49  		panic(badVn2)
    50  	}
    51  	if len(work) < n {
    52  		panic(badWork)
    53  	}
    54  
    55  	tol3z := math.Sqrt(dlamchE)
    56  
    57  	bi := blas64.Implementation()
    58  
    59  	// Compute factorization.
    60  	for i := 0; i < mn; i++ {
    61  		offpi := offset + i
    62  
    63  		// Determine ith pivot column and swap if necessary.
    64  		p := i + bi.Idamax(n-i, vn1[i:], 1)
    65  		if p != i {
    66  			bi.Dswap(m, a[p:], lda, a[i:], lda)
    67  			jpvt[p], jpvt[i] = jpvt[i], jpvt[p]
    68  			vn1[p] = vn1[i]
    69  			vn2[p] = vn2[i]
    70  		}
    71  
    72  		// Generate elementary reflector H_i.
    73  		if offpi < m-1 {
    74  			a[offpi*lda+i], tau[i] = impl.Dlarfg(m-offpi, a[offpi*lda+i], a[(offpi+1)*lda+i:], lda)
    75  		} else {
    76  			tau[i] = 0
    77  		}
    78  
    79  		if i < n-1 {
    80  			// Apply H_i^T to A[offset+i:m, i:n] from the left.
    81  			aii := a[offpi*lda+i]
    82  			a[offpi*lda+i] = 1
    83  			impl.Dlarf(blas.Left, m-offpi, n-i-1, a[offpi*lda+i:], lda, tau[i], a[offpi*lda+i+1:], lda, work)
    84  			a[offpi*lda+i] = aii
    85  		}
    86  
    87  		// Update partial column norms.
    88  		for j := i + 1; j < n; j++ {
    89  			if vn1[j] == 0 {
    90  				continue
    91  			}
    92  
    93  			// The following marked lines follow from the
    94  			// analysis in Lapack Working Note 176.
    95  			r := math.Abs(a[offpi*lda+j]) / vn1[j] // *
    96  			temp := math.Max(0, 1-r*r)             // *
    97  			r = vn1[j] / vn2[j]                    // *
    98  			temp2 := temp * r * r                  // *
    99  			if temp2 < tol3z {
   100  				var v float64
   101  				if offpi < m-1 {
   102  					v = bi.Dnrm2(m-offpi-1, a[(offpi+1)*lda+j:], lda)
   103  				}
   104  				vn1[j] = v
   105  				vn2[j] = v
   106  			} else {
   107  				vn1[j] *= math.Sqrt(temp) // *
   108  			}
   109  		}
   110  	}
   111  }