gonum.org/v1/gonum@v0.14.0/lapack/gonum/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 gonum
     6  
     7  import (
     8  	"gonum.org/v1/gonum/blas"
     9  	"gonum.org/v1/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  //
    15  //	C = Q * C   if side == blas.Left and trans == blas.NoTrans
    16  //	C = Qᵀ * C  if side == blas.Left and trans == blas.Trans
    17  //	C = C * Q   if side == blas.Right and trans == blas.NoTrans
    18  //	C = C * Qᵀ  if side == blas.Right and trans == blas.Trans
    19  //
    20  // If side == blas.Left, A is a matrix of side k×m, and if side == blas.Right
    21  // A is of size k×n. This uses a blocked algorithm.
    22  //
    23  // work is temporary storage, and lwork specifies the usable memory length.
    24  // At minimum, lwork >= m if side == blas.Left and lwork >= n if side == blas.Right,
    25  // and this function will panic otherwise.
    26  // Dormlq uses a block algorithm, but the block size is limited
    27  // by the temporary space available. If lwork == -1, instead of performing Dormlq,
    28  // the optimal work length will be stored into work[0].
    29  //
    30  // tau contains the Householder scales and must have length at least k, and
    31  // this function will panic otherwise.
    32  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) {
    33  	left := side == blas.Left
    34  	nw := m
    35  	if left {
    36  		nw = n
    37  	}
    38  	switch {
    39  	case !left && side != blas.Right:
    40  		panic(badSide)
    41  	case trans != blas.Trans && trans != blas.NoTrans:
    42  		panic(badTrans)
    43  	case m < 0:
    44  		panic(mLT0)
    45  	case n < 0:
    46  		panic(nLT0)
    47  	case k < 0:
    48  		panic(kLT0)
    49  	case left && k > m:
    50  		panic(kGTM)
    51  	case !left && k > n:
    52  		panic(kGTN)
    53  	case left && lda < max(1, m):
    54  		panic(badLdA)
    55  	case !left && lda < max(1, n):
    56  		panic(badLdA)
    57  	case lwork < max(1, nw) && lwork != -1:
    58  		panic(badLWork)
    59  	case len(work) < max(1, lwork):
    60  		panic(shortWork)
    61  	}
    62  
    63  	// Quick return if possible.
    64  	if m == 0 || n == 0 || k == 0 {
    65  		work[0] = 1
    66  		return
    67  	}
    68  
    69  	const (
    70  		nbmax = 64
    71  		ldt   = nbmax
    72  		tsize = nbmax * ldt
    73  	)
    74  	opts := string(side) + string(trans)
    75  	nb := min(nbmax, impl.Ilaenv(1, "DORMLQ", opts, m, n, k, -1))
    76  	lworkopt := max(1, nw)*nb + tsize
    77  	if lwork == -1 {
    78  		work[0] = float64(lworkopt)
    79  		return
    80  	}
    81  
    82  	switch {
    83  	case left && len(a) < (k-1)*lda+m:
    84  		panic(shortA)
    85  	case !left && len(a) < (k-1)*lda+n:
    86  		panic(shortA)
    87  	case len(tau) < k:
    88  		panic(shortTau)
    89  	case len(c) < (m-1)*ldc+n:
    90  		panic(shortC)
    91  	}
    92  
    93  	nbmin := 2
    94  	if 1 < nb && nb < k {
    95  		iws := nw*nb + tsize
    96  		if lwork < iws {
    97  			nb = (lwork - tsize) / nw
    98  			nbmin = max(2, impl.Ilaenv(2, "DORMLQ", opts, m, n, k, -1))
    99  		}
   100  	}
   101  	if nb < nbmin || k <= nb {
   102  		// Call unblocked code.
   103  		impl.Dorml2(side, trans, m, n, k, a, lda, tau, c, ldc, work)
   104  		work[0] = float64(lworkopt)
   105  		return
   106  	}
   107  
   108  	t := work[:tsize]
   109  	wrk := work[tsize:]
   110  	ldwrk := nb
   111  
   112  	notrans := trans == blas.NoTrans
   113  	transt := blas.NoTrans
   114  	if notrans {
   115  		transt = blas.Trans
   116  	}
   117  
   118  	switch {
   119  	case left && notrans:
   120  		for i := 0; i < k; i += nb {
   121  			ib := min(nb, k-i)
   122  			impl.Dlarft(lapack.Forward, lapack.RowWise, m-i, ib,
   123  				a[i*lda+i:], lda,
   124  				tau[i:],
   125  				t, ldt)
   126  			impl.Dlarfb(side, transt, lapack.Forward, lapack.RowWise, m-i, n, ib,
   127  				a[i*lda+i:], lda,
   128  				t, ldt,
   129  				c[i*ldc:], ldc,
   130  				wrk, ldwrk)
   131  		}
   132  
   133  	case left && !notrans:
   134  		for i := ((k - 1) / nb) * nb; i >= 0; i -= nb {
   135  			ib := min(nb, k-i)
   136  			impl.Dlarft(lapack.Forward, lapack.RowWise, m-i, ib,
   137  				a[i*lda+i:], lda,
   138  				tau[i:],
   139  				t, ldt)
   140  			impl.Dlarfb(side, transt, lapack.Forward, lapack.RowWise, m-i, n, ib,
   141  				a[i*lda+i:], lda,
   142  				t, ldt,
   143  				c[i*ldc:], ldc,
   144  				wrk, ldwrk)
   145  		}
   146  
   147  	case !left && notrans:
   148  		for i := ((k - 1) / nb) * nb; i >= 0; i -= nb {
   149  			ib := min(nb, k-i)
   150  			impl.Dlarft(lapack.Forward, lapack.RowWise, n-i, ib,
   151  				a[i*lda+i:], lda,
   152  				tau[i:],
   153  				t, ldt)
   154  			impl.Dlarfb(side, transt, lapack.Forward, lapack.RowWise, m, n-i, ib,
   155  				a[i*lda+i:], lda,
   156  				t, ldt,
   157  				c[i:], ldc,
   158  				wrk, ldwrk)
   159  		}
   160  
   161  	case !left && !notrans:
   162  		for i := 0; i < k; i += nb {
   163  			ib := min(nb, k-i)
   164  			impl.Dlarft(lapack.Forward, lapack.RowWise, n-i, ib,
   165  				a[i*lda+i:], lda,
   166  				tau[i:],
   167  				t, ldt)
   168  			impl.Dlarfb(side, transt, lapack.Forward, lapack.RowWise, m, n-i, ib,
   169  				a[i*lda+i:], lda,
   170  				t, ldt,
   171  				c[i:], ldc,
   172  				wrk, ldwrk)
   173  		}
   174  	}
   175  	work[0] = float64(lworkopt)
   176  }