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