github.com/gonum/lapack@v0.0.0-20181123203213-e4cdc5a0bff9/native/dormlq.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/lapack"
    10  )
    11  
    12  // Dormlq multiplies the matrix C by the orthogonal matrix Q defined by the
    13  // slices a and tau. A and tau are as returned from Dgelqf.
    14  //  C = Q * C    if side == blas.Left and trans == blas.NoTrans
    15  //  C = Q^T * C  if side == blas.Left and trans == blas.Trans
    16  //  C = C * Q    if side == blas.Right and trans == blas.NoTrans
    17  //  C = C * Q^T  if side == blas.Right and trans == blas.Trans
    18  // If side == blas.Left, A is a matrix of side k×m, and if side == blas.Right
    19  // A is of size k×n. This uses a blocked algorithm.
    20  //
    21  // work is temporary storage, and lwork specifies the usable memory length.
    22  // At minimum, lwork >= m if side == blas.Left and lwork >= n if side == blas.Right,
    23  // and this function will panic otherwise.
    24  // Dormlq uses a block algorithm, but the block size is limited
    25  // by the temporary space available. If lwork == -1, instead of performing Dormlq,
    26  // the optimal work length will be stored into work[0].
    27  //
    28  // tau contains the Householder scales and must have length at least k, and
    29  // this function will panic otherwise.
    30  func (impl Implementation) Dormlq(side blas.Side, trans blas.Transpose, m, n, k int, a []float64, lda int, tau, c []float64, ldc int, work []float64, lwork int) {
    31  	if side != blas.Left && side != blas.Right {
    32  		panic(badSide)
    33  	}
    34  	if trans != blas.Trans && trans != blas.NoTrans {
    35  		panic(badTrans)
    36  	}
    37  	left := side == blas.Left
    38  	if left {
    39  		checkMatrix(k, m, a, lda)
    40  	} else {
    41  		checkMatrix(k, n, a, lda)
    42  	}
    43  	checkMatrix(m, n, c, ldc)
    44  	if len(tau) < k {
    45  		panic(badTau)
    46  	}
    47  	if len(work) < lwork {
    48  		panic(shortWork)
    49  	}
    50  	nw := m
    51  	if left {
    52  		nw = n
    53  	}
    54  	if lwork < max(1, nw) && lwork != -1 {
    55  		panic(badWork)
    56  	}
    57  
    58  	if m == 0 || n == 0 || k == 0 {
    59  		work[0] = 1
    60  		return
    61  	}
    62  
    63  	const (
    64  		nbmax = 64
    65  		ldt   = nbmax
    66  		tsize = nbmax * ldt
    67  	)
    68  	opts := string(side) + string(trans)
    69  	nb := min(nbmax, impl.Ilaenv(1, "DORMLQ", opts, m, n, k, -1))
    70  	lworkopt := max(1, nw)*nb + tsize
    71  	if lwork == -1 {
    72  		work[0] = float64(lworkopt)
    73  		return
    74  	}
    75  
    76  	nbmin := 2
    77  	if 1 < nb && nb < k {
    78  		iws := nw*nb + tsize
    79  		if lwork < iws {
    80  			nb = (lwork - tsize) / nw
    81  			nbmin = max(2, impl.Ilaenv(2, "DORMLQ", opts, m, n, k, -1))
    82  		}
    83  	}
    84  	if nb < nbmin || k <= nb {
    85  		// Call unblocked code.
    86  		impl.Dorml2(side, trans, m, n, k, a, lda, tau, c, ldc, work)
    87  		work[0] = float64(lworkopt)
    88  		return
    89  	}
    90  
    91  	t := work[:tsize]
    92  	wrk := work[tsize:]
    93  	ldwrk := nb
    94  
    95  	notran := trans == blas.NoTrans
    96  	transt := blas.NoTrans
    97  	if notran {
    98  		transt = blas.Trans
    99  	}
   100  
   101  	switch {
   102  	case left && notran:
   103  		for i := 0; i < k; i += nb {
   104  			ib := min(nb, k-i)
   105  			impl.Dlarft(lapack.Forward, lapack.RowWise, m-i, ib,
   106  				a[i*lda+i:], lda,
   107  				tau[i:],
   108  				t, ldt)
   109  			impl.Dlarfb(side, transt, lapack.Forward, lapack.RowWise, m-i, n, ib,
   110  				a[i*lda+i:], lda,
   111  				t, ldt,
   112  				c[i*ldc:], ldc,
   113  				wrk, ldwrk)
   114  		}
   115  
   116  	case left && !notran:
   117  		for i := ((k - 1) / nb) * nb; i >= 0; i -= nb {
   118  			ib := min(nb, k-i)
   119  			impl.Dlarft(lapack.Forward, lapack.RowWise, m-i, ib,
   120  				a[i*lda+i:], lda,
   121  				tau[i:],
   122  				t, ldt)
   123  			impl.Dlarfb(side, transt, lapack.Forward, lapack.RowWise, m-i, n, ib,
   124  				a[i*lda+i:], lda,
   125  				t, ldt,
   126  				c[i*ldc:], ldc,
   127  				wrk, ldwrk)
   128  		}
   129  
   130  	case !left && notran:
   131  		for i := ((k - 1) / nb) * nb; i >= 0; i -= nb {
   132  			ib := min(nb, k-i)
   133  			impl.Dlarft(lapack.Forward, lapack.RowWise, n-i, ib,
   134  				a[i*lda+i:], lda,
   135  				tau[i:],
   136  				t, ldt)
   137  			impl.Dlarfb(side, transt, lapack.Forward, lapack.RowWise, m, n-i, ib,
   138  				a[i*lda+i:], lda,
   139  				t, ldt,
   140  				c[i:], ldc,
   141  				wrk, ldwrk)
   142  		}
   143  
   144  	case !left && !notran:
   145  		for i := 0; i < k; i += nb {
   146  			ib := min(nb, k-i)
   147  			impl.Dlarft(lapack.Forward, lapack.RowWise, n-i, ib,
   148  				a[i*lda+i:], lda,
   149  				tau[i:],
   150  				t, ldt)
   151  			impl.Dlarfb(side, transt, lapack.Forward, lapack.RowWise, m, n-i, ib,
   152  				a[i*lda+i:], lda,
   153  				t, ldt,
   154  				c[i:], ldc,
   155  				wrk, ldwrk)
   156  		}
   157  	}
   158  	work[0] = float64(lworkopt)
   159  }