github.com/gonum/lapack@v0.0.0-20181123203213-e4cdc5a0bff9/native/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 native
     6  
     7  import (
     8  	"github.com/gonum/blas"
     9  	"github.com/gonum/blas/blas64"
    10  )
    11  
    12  // Dlarf applies an elementary reflector to a general rectangular matrix c.
    13  // This computes
    14  //  c = h * c if side == Left
    15  //  c = c * h if side == right
    16  // where
    17  //  h = 1 - tau * v * v^T
    18  // and c is an m * n matrix.
    19  //
    20  // work is temporary storage of length at least m if side == Left and at least
    21  // n if side == Right. This function will panic if this length requirement is not met.
    22  //
    23  // Dlarf is an internal routine. It is exported for testing purposes.
    24  func (impl Implementation) Dlarf(side blas.Side, m, n int, v []float64, incv int, tau float64, c []float64, ldc int, work []float64) {
    25  	applyleft := side == blas.Left
    26  	if (applyleft && len(work) < n) || (!applyleft && len(work) < m) {
    27  		panic(badWork)
    28  	}
    29  	checkMatrix(m, n, c, ldc)
    30  
    31  	// v has length m if applyleft and n otherwise.
    32  	lenV := n
    33  	if applyleft {
    34  		lenV = m
    35  	}
    36  
    37  	checkVector(lenV, v, incv)
    38  
    39  	lastv := 0 // last non-zero element of v
    40  	lastc := 0 // last non-zero row/column of c
    41  	if tau != 0 {
    42  		var i int
    43  		if applyleft {
    44  			lastv = m - 1
    45  		} else {
    46  			lastv = n - 1
    47  		}
    48  		if incv > 0 {
    49  			i = lastv * incv
    50  		}
    51  
    52  		// Look for the last non-zero row in v.
    53  		for lastv >= 0 && v[i] == 0 {
    54  			lastv--
    55  			i -= incv
    56  		}
    57  		if applyleft {
    58  			// Scan for the last non-zero column in C[0:lastv, :]
    59  			lastc = impl.Iladlc(lastv+1, n, c, ldc)
    60  		} else {
    61  			// Scan for the last non-zero row in C[:, 0:lastv]
    62  			lastc = impl.Iladlr(m, lastv+1, c, ldc)
    63  		}
    64  	}
    65  	if lastv == -1 || lastc == -1 {
    66  		return
    67  	}
    68  	// Sometimes 1-indexing is nicer ...
    69  	bi := blas64.Implementation()
    70  	if applyleft {
    71  		// Form H * C
    72  		// w[0:lastc+1] = c[1:lastv+1, 1:lastc+1]^T * v[1:lastv+1,1]
    73  		bi.Dgemv(blas.Trans, lastv+1, lastc+1, 1, c, ldc, v, incv, 0, work, 1)
    74  		// c[0: lastv, 0: lastc] = c[...] - w[0:lastv, 1] * v[1:lastc, 1]^T
    75  		bi.Dger(lastv+1, lastc+1, -tau, v, incv, work, 1, c, ldc)
    76  		return
    77  	}
    78  	// Form C*H
    79  	// w[0:lastc+1,1] := c[0:lastc+1,0:lastv+1] * v[0:lastv+1,1]
    80  	bi.Dgemv(blas.NoTrans, lastc+1, lastv+1, 1, c, ldc, v, incv, 0, work, 1)
    81  	// c[0:lastc+1,0:lastv+1] = c[...] - w[0:lastc+1,0] * v[0:lastv+1,0]^T
    82  	bi.Dger(lastc+1, lastv+1, -tau, work, 1, v, incv, c, ldc)
    83  }