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