vendor/gonum.org/v1/gonum/lapack/gonum/dlaqp2.go

   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         "gonum.org/v1/gonum/blas"
  11         "gonum.org/v1/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 }