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
5package gonum
6
7import (
8	"gonum.org/v1/gonum/blas"
9	"gonum.org/v1/gonum/lapack"
10)
11
12// Dormbr applies a multiplicative update to the matrix C based on a
13// decomposition computed by Dgebrd.
14//
15// Dormbr overwrites the m×n matrix C with
16//  Q * C   if vect == lapack.ApplyQ, side == blas.Left, and trans == blas.NoTrans
17//  C * Q   if vect == lapack.ApplyQ, side == blas.Right, and trans == blas.NoTrans
18//  Q^T * C if vect == lapack.ApplyQ, side == blas.Left, and trans == blas.Trans
19//  C * Q^T if vect == lapack.ApplyQ, side == blas.Right, and trans == blas.Trans
20//
21//  P * C   if vect == lapack.ApplyP, side == blas.Left, and trans == blas.NoTrans
22//  C * P   if vect == lapack.ApplyP, side == blas.Right, and trans == blas.NoTrans
23//  P^T * C if vect == lapack.ApplyP, side == blas.Left, and trans == blas.Trans
24//  C * P^T if vect == lapack.ApplyP, side == blas.Right, and trans == blas.Trans
25// where P and Q are the orthogonal matrices determined by Dgebrd when reducing
26// a matrix A to bidiagonal form: A = Q * B * P^T. See Dgebrd for the
27// definitions of Q and P.
28//
29// If vect == lapack.ApplyQ, A is assumed to have been an nq×k matrix, while if
30// vect == lapack.ApplyP, A is assumed to have been a k×nq matrix. nq = m if
31// side == blas.Left, while nq = n if side == blas.Right.
32//
33// tau must have length min(nq,k), and Dormbr will panic otherwise. tau contains
34// the elementary reflectors to construct Q or P depending on the value of
35// vect.
36//
37// work must have length at least max(1,lwork), and lwork must be either -1 or
38// at least max(1,n) if side == blas.Left, and at least max(1,m) if side ==
39// blas.Right. For optimum performance lwork should be at least n*nb if side ==
40// blas.Left, and at least m*nb if side == blas.Right, where nb is the optimal
41// block size. On return, work[0] will contain the optimal value of lwork.
42//
43// If lwork == -1, the function only calculates the optimal value of lwork and
44// returns it in work[0].
45//
46// Dormbr is an internal routine. It is exported for testing purposes.
47func (impl Implementation) Dormbr(vect lapack.ApplyOrtho, side blas.Side, trans blas.Transpose, m, n, k int, a []float64, lda int, tau, c []float64, ldc int, work []float64, lwork int) {
48	nq := n
49	nw := m
50	if side == blas.Left {
51		nq = m
52		nw = n
53	}
54	applyQ := vect == lapack.ApplyQ
55	switch {
56	case !applyQ && vect != lapack.ApplyP:
57		panic(badApplyOrtho)
58	case side != blas.Left && side != blas.Right:
59		panic(badSide)
60	case trans != blas.NoTrans && trans != blas.Trans:
61		panic(badTrans)
62	case m < 0:
63		panic(mLT0)
64	case n < 0:
65		panic(nLT0)
66	case k < 0:
67		panic(kLT0)
68	case applyQ && lda < max(1, min(nq, k)):
69		panic(badLdA)
70	case !applyQ && lda < max(1, nq):
71		panic(badLdA)
72	case ldc < max(1, n):
73		panic(badLdC)
74	case lwork < max(1, nw) && lwork != -1:
75		panic(badLWork)
76	case len(work) < max(1, lwork):
77		panic(shortWork)
78	}
79
80	// Quick return if possible.
81	if m == 0 || n == 0 {
82		work[0] = 1
83		return
84	}
85
86	// The current implementation does not use opts, but a future change may
87	// use these options so construct them.
88	var opts string
89	if side == blas.Left {
90		opts = "L"
91	} else {
92		opts = "R"
93	}
94	if trans == blas.Trans {
95		opts += "T"
96	} else {
97		opts += "N"
98	}
99	var nb int
100	if applyQ {
101		if side == blas.Left {
102			nb = impl.Ilaenv(1, "DORMQR", opts, m-1, n, m-1, -1)
103		} else {
104			nb = impl.Ilaenv(1, "DORMQR", opts, m, n-1, n-1, -1)
105		}
106	} else {
107		if side == blas.Left {
108			nb = impl.Ilaenv(1, "DORMLQ", opts, m-1, n, m-1, -1)
109		} else {
110			nb = impl.Ilaenv(1, "DORMLQ", opts, m, n-1, n-1, -1)
111		}
112	}
113	lworkopt := max(1, nw) * nb
114	if lwork == -1 {
115		work[0] = float64(lworkopt)
116		return
117	}
118
119	minnqk := min(nq, k)
120	switch {
121	case applyQ && len(a) < (nq-1)*lda+minnqk:
122		panic(shortA)
123	case !applyQ && len(a) < (minnqk-1)*lda+nq:
124		panic(shortA)
125	case len(tau) < minnqk:
126		panic(shortTau)
127	case len(c) < (m-1)*ldc+n:
128		panic(shortC)
129	}
130
131	if applyQ {
132		// Change the operation to get Q depending on the size of the initial
133		// matrix to Dgebrd. The size matters due to the storage location of
134		// the off-diagonal elements.
135		if nq >= k {
136			impl.Dormqr(side, trans, m, n, k, a, lda, tau[:k], c, ldc, work, lwork)
137		} else if nq > 1 {
138			mi := m
139			ni := n - 1
140			i1 := 0
141			i2 := 1
142			if side == blas.Left {
143				mi = m - 1
144				ni = n
145				i1 = 1
146				i2 = 0
147			}
148			impl.Dormqr(side, trans, mi, ni, nq-1, a[1*lda:], lda, tau[:nq-1], c[i1*ldc+i2:], ldc, work, lwork)
149		}
150		work[0] = float64(lworkopt)
151		return
152	}
153
154	transt := blas.Trans
155	if trans == blas.Trans {
156		transt = blas.NoTrans
157	}
158
159	// Change the operation to get P depending on the size of the initial
160	// matrix to Dgebrd. The size matters due to the storage location of
161	// the off-diagonal elements.
162	if nq > k {
163		impl.Dormlq(side, transt, m, n, k, a, lda, tau, c, ldc, work, lwork)
164	} else if nq > 1 {
165		mi := m
166		ni := n - 1
167		i1 := 0
168		i2 := 1
169		if side == blas.Left {
170			mi = m - 1
171			ni = n
172			i1 = 1
173			i2 = 0
174		}
175		impl.Dormlq(side, transt, mi, ni, nq-1, a[1:], lda, tau, c[i1*ldc+i2:], ldc, work, lwork)
176	}
177	work[0] = float64(lworkopt)
178}
179