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