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