gonum.org/v1/gonum@v0.14.0/lapack/gonum/dormhr.go (about)

     1  // Copyright ©2016 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 "gonum.org/v1/gonum/blas"
     8  
     9  // Dormhr multiplies an m×n general matrix C with an nq×nq orthogonal matrix Q
    10  //
    11  //	Q * C   if side == blas.Left  and trans == blas.NoTrans,
    12  //	Qᵀ * C  if side == blas.Left  and trans == blas.Trans,
    13  //	C * Q   if side == blas.Right and trans == blas.NoTrans,
    14  //	C * Qᵀ  if side == blas.Right and trans == blas.Trans,
    15  //
    16  // where nq == m if side == blas.Left and nq == n if side == blas.Right.
    17  //
    18  // Q is defined implicitly as the product of ihi-ilo elementary reflectors, as
    19  // returned by Dgehrd:
    20  //
    21  //	Q = H_{ilo} H_{ilo+1} ... H_{ihi-1}.
    22  //
    23  // Q is equal to the identity matrix except in the submatrix
    24  // Q[ilo+1:ihi+1,ilo+1:ihi+1].
    25  //
    26  // ilo and ihi must have the same values as in the previous call of Dgehrd. It
    27  // must hold that
    28  //
    29  //	0 <= ilo <= ihi < m   if m > 0 and side == blas.Left,
    30  //	ilo = 0 and ihi = -1  if m = 0 and side == blas.Left,
    31  //	0 <= ilo <= ihi < n   if n > 0 and side == blas.Right,
    32  //	ilo = 0 and ihi = -1  if n = 0 and side == blas.Right.
    33  //
    34  // a and lda represent an m×m matrix if side == blas.Left and an n×n matrix if
    35  // side == blas.Right. The matrix contains vectors which define the elementary
    36  // reflectors, as returned by Dgehrd.
    37  //
    38  // tau contains the scalar factors of the elementary reflectors, as returned by
    39  // Dgehrd. tau must have length m-1 if side == blas.Left and n-1 if side ==
    40  // blas.Right.
    41  //
    42  // c and ldc represent the m×n matrix C. On return, c is overwritten by the
    43  // product with Q.
    44  //
    45  // work must have length at least max(1,lwork), and lwork must be at least
    46  // max(1,n), if side == blas.Left, and max(1,m), if side == blas.Right. For
    47  // optimum performance lwork should be at least n*nb if side == blas.Left and
    48  // m*nb if side == blas.Right, where nb is the optimal block size. On return,
    49  // work[0] will contain the optimal value of lwork.
    50  //
    51  // If lwork == -1, instead of performing Dormhr, only the optimal value of lwork
    52  // will be stored in work[0].
    53  //
    54  // If any requirement on input sizes is not met, Dormhr will panic.
    55  //
    56  // Dormhr is an internal routine. It is exported for testing purposes.
    57  func (impl Implementation) Dormhr(side blas.Side, trans blas.Transpose, m, n, ilo, ihi int, a []float64, lda int, tau, c []float64, ldc int, work []float64, lwork int) {
    58  	nq := n // The order of Q.
    59  	nw := m // The minimum length of work.
    60  	if side == blas.Left {
    61  		nq = m
    62  		nw = n
    63  	}
    64  	switch {
    65  	case side != blas.Left && side != blas.Right:
    66  		panic(badSide)
    67  	case trans != blas.NoTrans && trans != blas.Trans:
    68  		panic(badTrans)
    69  	case m < 0:
    70  		panic(mLT0)
    71  	case n < 0:
    72  		panic(nLT0)
    73  	case ilo < 0 || max(1, nq) <= ilo:
    74  		panic(badIlo)
    75  	case ihi < min(ilo, nq-1) || nq <= ihi:
    76  		panic(badIhi)
    77  	case lda < max(1, nq):
    78  		panic(badLdA)
    79  	case lwork < max(1, nw) && lwork != -1:
    80  		panic(badLWork)
    81  	case len(work) < max(1, lwork):
    82  		panic(shortWork)
    83  	}
    84  
    85  	// Quick return if possible.
    86  	if m == 0 || n == 0 {
    87  		work[0] = 1
    88  		return
    89  	}
    90  
    91  	nh := ihi - ilo
    92  	var nb int
    93  	if side == blas.Left {
    94  		opts := "LN"
    95  		if trans == blas.Trans {
    96  			opts = "LT"
    97  		}
    98  		nb = impl.Ilaenv(1, "DORMQR", opts, nh, n, nh, -1)
    99  	} else {
   100  		opts := "RN"
   101  		if trans == blas.Trans {
   102  			opts = "RT"
   103  		}
   104  		nb = impl.Ilaenv(1, "DORMQR", opts, m, nh, nh, -1)
   105  	}
   106  	lwkopt := max(1, nw) * nb
   107  	if lwork == -1 {
   108  		work[0] = float64(lwkopt)
   109  		return
   110  	}
   111  
   112  	if nh == 0 {
   113  		work[0] = 1
   114  		return
   115  	}
   116  
   117  	switch {
   118  	case len(a) < (nq-1)*lda+nq:
   119  		panic(shortA)
   120  	case len(c) < (m-1)*ldc+n:
   121  		panic(shortC)
   122  	case len(tau) != nq-1:
   123  		panic(badLenTau)
   124  	}
   125  
   126  	if side == blas.Left {
   127  		impl.Dormqr(side, trans, nh, n, nh, a[(ilo+1)*lda+ilo:], lda,
   128  			tau[ilo:ihi], c[(ilo+1)*ldc:], ldc, work, lwork)
   129  	} else {
   130  		impl.Dormqr(side, trans, m, nh, nh, a[(ilo+1)*lda+ilo:], lda,
   131  			tau[ilo:ihi], c[ilo+1:], ldc, work, lwork)
   132  	}
   133  	work[0] = float64(lwkopt)
   134  }