gonum.org/v1/gonum@v0.14.0/lapack/gonum/dlarf.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  	"gonum.org/v1/gonum/blas"
     9  	"gonum.org/v1/gonum/blas/blas64"
    10  )
    11  
    12  // Dlarf applies an elementary reflector H to an m×n matrix C:
    13  //
    14  //	C = H * C  if side == blas.Left
    15  //	C = C * H  if side == blas.Right
    16  //
    17  // H is represented in the form
    18  //
    19  //	H = I - tau * v * vᵀ
    20  //
    21  // where tau is a scalar and v is a vector.
    22  //
    23  // work must have length at least m if side == blas.Left and
    24  // at least n if side == blas.Right.
    25  //
    26  // Dlarf is an internal routine. It is exported for testing purposes.
    27  func (impl Implementation) Dlarf(side blas.Side, m, n int, v []float64, incv int, tau float64, c []float64, ldc int, work []float64) {
    28  	switch {
    29  	case side != blas.Left && side != blas.Right:
    30  		panic(badSide)
    31  	case m < 0:
    32  		panic(mLT0)
    33  	case n < 0:
    34  		panic(nLT0)
    35  	case incv == 0:
    36  		panic(zeroIncV)
    37  	case ldc < max(1, n):
    38  		panic(badLdC)
    39  	}
    40  
    41  	if m == 0 || n == 0 {
    42  		return
    43  	}
    44  
    45  	applyleft := side == blas.Left
    46  	lenV := n
    47  	if applyleft {
    48  		lenV = m
    49  	}
    50  
    51  	switch {
    52  	case len(v) < 1+(lenV-1)*abs(incv):
    53  		panic(shortV)
    54  	case len(c) < (m-1)*ldc+n:
    55  		panic(shortC)
    56  	case (applyleft && len(work) < n) || (!applyleft && len(work) < m):
    57  		panic(shortWork)
    58  	}
    59  
    60  	lastv := -1 // last non-zero element of v
    61  	lastc := -1 // last non-zero row/column of C
    62  	if tau != 0 {
    63  		if applyleft {
    64  			lastv = m - 1
    65  		} else {
    66  			lastv = n - 1
    67  		}
    68  		var i int
    69  		if incv > 0 {
    70  			i = lastv * incv
    71  		}
    72  		// Look for the last non-zero row in v.
    73  		for lastv >= 0 && v[i] == 0 {
    74  			lastv--
    75  			i -= incv
    76  		}
    77  		if applyleft {
    78  			// Scan for the last non-zero column in C[0:lastv, :]
    79  			lastc = impl.Iladlc(lastv+1, n, c, ldc)
    80  		} else {
    81  			// Scan for the last non-zero row in C[:, 0:lastv]
    82  			lastc = impl.Iladlr(m, lastv+1, c, ldc)
    83  		}
    84  	}
    85  	if lastv == -1 || lastc == -1 {
    86  		return
    87  	}
    88  	bi := blas64.Implementation()
    89  	if applyleft {
    90  		// Form H * C
    91  		// w[0:lastc+1] = c[1:lastv+1, 1:lastc+1]ᵀ * v[1:lastv+1,1]
    92  		bi.Dgemv(blas.Trans, lastv+1, lastc+1, 1, c, ldc, v, incv, 0, work, 1)
    93  		// c[0: lastv, 0: lastc] = c[...] - w[0:lastv, 1] * v[1:lastc, 1]ᵀ
    94  		bi.Dger(lastv+1, lastc+1, -tau, v, incv, work, 1, c, ldc)
    95  	} else {
    96  		// Form C * H
    97  		// w[0:lastc+1,1] := c[0:lastc+1,0:lastv+1] * v[0:lastv+1,1]
    98  		bi.Dgemv(blas.NoTrans, lastc+1, lastv+1, 1, c, ldc, v, incv, 0, work, 1)
    99  		// c[0:lastc+1,0:lastv+1] = c[...] - w[0:lastc+1,0] * v[0:lastv+1,0]ᵀ
   100  		bi.Dger(lastc+1, lastv+1, -tau, work, 1, v, incv, c, ldc)
   101  	}
   102  }