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

     1  // Copyright ©2015 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  	"github.com/jingcheng-WU/gonum/blas"
     9  	"github.com/jingcheng-WU/gonum/blas/blas64"
    10  	"github.com/jingcheng-WU/gonum/lapack"
    11  )
    12  
    13  // Dlarft forms the triangular factor T of a block reflector H, storing the answer
    14  // in t.
    15  //  H = I - V * T * Vᵀ  if store == lapack.ColumnWise
    16  //  H = I - Vᵀ * T * V  if store == lapack.RowWise
    17  // H is defined by a product of the elementary reflectors where
    18  //  H = H_0 * H_1 * ... * H_{k-1}  if direct == lapack.Forward
    19  //  H = H_{k-1} * ... * H_1 * H_0  if direct == lapack.Backward
    20  //
    21  // t is a k×k triangular matrix. t is upper triangular if direct = lapack.Forward
    22  // and lower triangular otherwise. This function will panic if t is not of
    23  // sufficient size.
    24  //
    25  // store describes the storage of the elementary reflectors in v. See
    26  // Dlarfb for a description of layout.
    27  //
    28  // tau contains the scalar factors of the elementary reflectors H_i.
    29  //
    30  // Dlarft is an internal routine. It is exported for testing purposes.
    31  func (Implementation) Dlarft(direct lapack.Direct, store lapack.StoreV, n, k int, v []float64, ldv int, tau []float64, t []float64, ldt int) {
    32  	mv, nv := n, k
    33  	if store == lapack.RowWise {
    34  		mv, nv = k, n
    35  	}
    36  	switch {
    37  	case direct != lapack.Forward && direct != lapack.Backward:
    38  		panic(badDirect)
    39  	case store != lapack.RowWise && store != lapack.ColumnWise:
    40  		panic(badStoreV)
    41  	case n < 0:
    42  		panic(nLT0)
    43  	case k < 1:
    44  		panic(kLT1)
    45  	case ldv < max(1, nv):
    46  		panic(badLdV)
    47  	case len(tau) < k:
    48  		panic(shortTau)
    49  	case ldt < max(1, k):
    50  		panic(shortT)
    51  	}
    52  
    53  	if n == 0 {
    54  		return
    55  	}
    56  
    57  	switch {
    58  	case len(v) < (mv-1)*ldv+nv:
    59  		panic(shortV)
    60  	case len(t) < (k-1)*ldt+k:
    61  		panic(shortT)
    62  	}
    63  
    64  	bi := blas64.Implementation()
    65  
    66  	// TODO(btracey): There are a number of minor obvious loop optimizations here.
    67  	// TODO(btracey): It may be possible to rearrange some of the code so that
    68  	// index of 1 is more common in the Dgemv.
    69  	if direct == lapack.Forward {
    70  		prevlastv := n - 1
    71  		for i := 0; i < k; i++ {
    72  			prevlastv = max(i, prevlastv)
    73  			if tau[i] == 0 {
    74  				for j := 0; j <= i; j++ {
    75  					t[j*ldt+i] = 0
    76  				}
    77  				continue
    78  			}
    79  			var lastv int
    80  			if store == lapack.ColumnWise {
    81  				// skip trailing zeros
    82  				for lastv = n - 1; lastv >= i+1; lastv-- {
    83  					if v[lastv*ldv+i] != 0 {
    84  						break
    85  					}
    86  				}
    87  				for j := 0; j < i; j++ {
    88  					t[j*ldt+i] = -tau[i] * v[i*ldv+j]
    89  				}
    90  				j := min(lastv, prevlastv)
    91  				bi.Dgemv(blas.Trans, j-i, i,
    92  					-tau[i], v[(i+1)*ldv:], ldv, v[(i+1)*ldv+i:], ldv,
    93  					1, t[i:], ldt)
    94  			} else {
    95  				for lastv = n - 1; lastv >= i+1; lastv-- {
    96  					if v[i*ldv+lastv] != 0 {
    97  						break
    98  					}
    99  				}
   100  				for j := 0; j < i; j++ {
   101  					t[j*ldt+i] = -tau[i] * v[j*ldv+i]
   102  				}
   103  				j := min(lastv, prevlastv)
   104  				bi.Dgemv(blas.NoTrans, i, j-i,
   105  					-tau[i], v[i+1:], ldv, v[i*ldv+i+1:], 1,
   106  					1, t[i:], ldt)
   107  			}
   108  			bi.Dtrmv(blas.Upper, blas.NoTrans, blas.NonUnit, i, t, ldt, t[i:], ldt)
   109  			t[i*ldt+i] = tau[i]
   110  			if i > 1 {
   111  				prevlastv = max(prevlastv, lastv)
   112  			} else {
   113  				prevlastv = lastv
   114  			}
   115  		}
   116  		return
   117  	}
   118  	prevlastv := 0
   119  	for i := k - 1; i >= 0; i-- {
   120  		if tau[i] == 0 {
   121  			for j := i; j < k; j++ {
   122  				t[j*ldt+i] = 0
   123  			}
   124  			continue
   125  		}
   126  		var lastv int
   127  		if i < k-1 {
   128  			if store == lapack.ColumnWise {
   129  				for lastv = 0; lastv < i; lastv++ {
   130  					if v[lastv*ldv+i] != 0 {
   131  						break
   132  					}
   133  				}
   134  				for j := i + 1; j < k; j++ {
   135  					t[j*ldt+i] = -tau[i] * v[(n-k+i)*ldv+j]
   136  				}
   137  				j := max(lastv, prevlastv)
   138  				bi.Dgemv(blas.Trans, n-k+i-j, k-i-1,
   139  					-tau[i], v[j*ldv+i+1:], ldv, v[j*ldv+i:], ldv,
   140  					1, t[(i+1)*ldt+i:], ldt)
   141  			} else {
   142  				for lastv = 0; lastv < i; lastv++ {
   143  					if v[i*ldv+lastv] != 0 {
   144  						break
   145  					}
   146  				}
   147  				for j := i + 1; j < k; j++ {
   148  					t[j*ldt+i] = -tau[i] * v[j*ldv+n-k+i]
   149  				}
   150  				j := max(lastv, prevlastv)
   151  				bi.Dgemv(blas.NoTrans, k-i-1, n-k+i-j,
   152  					-tau[i], v[(i+1)*ldv+j:], ldv, v[i*ldv+j:], 1,
   153  					1, t[(i+1)*ldt+i:], ldt)
   154  			}
   155  			bi.Dtrmv(blas.Lower, blas.NoTrans, blas.NonUnit, k-i-1,
   156  				t[(i+1)*ldt+i+1:], ldt,
   157  				t[(i+1)*ldt+i:], ldt)
   158  			if i > 0 {
   159  				prevlastv = min(prevlastv, lastv)
   160  			} else {
   161  				prevlastv = lastv
   162  			}
   163  		}
   164  		t[i*ldt+i] = tau[i]
   165  	}
   166  }