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.
8 "gonum.org/v1/gonum/blas"
9 "gonum.org/v1/gonum/lapack"
12 // Dlasr applies a sequence of plane rotations to the m×n matrix A. This series
13 // of plane rotations is implicitly represented by a matrix P. P is multiplied
14 // by a depending on the value of side -- A = P * A if side == lapack.Left,
15 // A = A * P^T if side == lapack.Right.
17 //The exact value of P depends on the value of pivot, but in all cases P is
18 // implicitly represented by a series of 2×2 rotation matrices. The entries of
19 // rotation matrix k are defined by s[k] and c[k]
20 // R(k) = [ c[k] s[k]]
22 // If direct == lapack.Forward, the rotation matrices are applied as
23 // P = P(z-1) * ... * P(2) * P(1), while if direct == lapack.Backward they are
24 // applied as P = P(1) * P(2) * ... * P(n).
26 // pivot defines the mapping of the elements in R(k) to P(k).
27 // If pivot == lapack.Variable, the rotation is performed for the (k, k+1) plane.
36 // if pivot == lapack.Top, the rotation is performed for the (1, k+1) plane,
37 // P(k) = [c[k] s[k] ]
45 // and if pivot == lapack.Bottom, the rotation is performed for the (k, z) plane.
54 // s and c have length m - 1 if side == blas.Left, and n - 1 if side == blas.Right.
56 // Dlasr is an internal routine. It is exported for testing purposes.
57 func (impl Implementation) Dlasr(side blas.Side, pivot lapack.Pivot, direct lapack.Direct, m, n int, c, s, a []float64, lda int) {
58 checkMatrix(m, n, a, lda)
59 if side != blas.Left && side != blas.Right {
62 if pivot != lapack.Variable && pivot != lapack.Top && pivot != lapack.Bottom {
65 if direct != lapack.Forward && direct != lapack.Backward {
68 if side == blas.Left {
86 if side == blas.Left {
87 if pivot == lapack.Variable {
88 if direct == lapack.Forward {
89 for j := 0; j < m-1; j++ {
92 if ctmp != 1 || stmp != 0 {
93 for i := 0; i < n; i++ {
96 a[(j+1)*lda+i] = ctmp*tmp - stmp*tmp2
97 a[j*lda+i] = stmp*tmp + ctmp*tmp2
103 for j := m - 2; j >= 0; j-- {
106 if ctmp != 1 || stmp != 0 {
107 for i := 0; i < n; i++ {
109 tmp := a[(j+1)*lda+i]
110 a[(j+1)*lda+i] = ctmp*tmp - stmp*tmp2
111 a[j*lda+i] = stmp*tmp + ctmp*tmp2
116 } else if pivot == lapack.Top {
117 if direct == lapack.Forward {
118 for j := 1; j < m; j++ {
121 if ctmp != 1 || stmp != 0 {
122 for i := 0; i < n; i++ {
125 a[j*lda+i] = ctmp*tmp - stmp*tmp2
126 a[i] = stmp*tmp + ctmp*tmp2
132 for j := m - 1; j >= 1; j-- {
135 if ctmp != 1 || stmp != 0 {
136 for i := 0; i < n; i++ {
139 if ctmp != 1 || stmp != 0 {
140 for i := 0; i < n; i++ {
143 a[j*lda+i] = ctmp*tmp - stmp*tmp2
144 a[i] = stmp*tmp + ctmp*tmp2
152 if direct == lapack.Forward {
153 for j := 0; j < m-1; j++ {
156 if ctmp != 1 || stmp != 0 {
157 for i := 0; i < n; i++ {
159 tmp2 := a[(m-1)*lda+i]
160 a[j*lda+i] = stmp*tmp2 + ctmp*tmp
161 a[(m-1)*lda+i] = ctmp*tmp2 - stmp*tmp
167 for j := m - 2; j >= 0; j-- {
170 if ctmp != 1 || stmp != 0 {
171 for i := 0; i < n; i++ {
173 tmp2 := a[(m-1)*lda+i]
174 a[j*lda+i] = stmp*tmp2 + ctmp*tmp
175 a[(m-1)*lda+i] = ctmp*tmp2 - stmp*tmp
181 if pivot == lapack.Variable {
182 if direct == lapack.Forward {
183 for j := 0; j < n-1; j++ {
186 if ctmp != 1 || stmp != 0 {
187 for i := 0; i < m; i++ {
190 a[i*lda+j+1] = ctmp*tmp - stmp*tmp2
191 a[i*lda+j] = stmp*tmp + ctmp*tmp2
197 for j := n - 2; j >= 0; j-- {
200 if ctmp != 1 || stmp != 0 {
201 for i := 0; i < m; i++ {
204 a[i*lda+j+1] = ctmp*tmp - stmp*tmp2
205 a[i*lda+j] = stmp*tmp + ctmp*tmp2
210 } else if pivot == lapack.Top {
211 if direct == lapack.Forward {
212 for j := 1; j < n; j++ {
215 if ctmp != 1 || stmp != 0 {
216 for i := 0; i < m; i++ {
219 a[i*lda+j] = ctmp*tmp - stmp*tmp2
220 a[i*lda] = stmp*tmp + ctmp*tmp2
226 for j := n - 1; j >= 1; j-- {
229 if ctmp != 1 || stmp != 0 {
230 for i := 0; i < m; i++ {
233 a[i*lda+j] = ctmp*tmp - stmp*tmp2
234 a[i*lda] = stmp*tmp + ctmp*tmp2
240 if direct == lapack.Forward {
241 for j := 0; j < n-1; j++ {
244 if ctmp != 1 || stmp != 0 {
245 for i := 0; i < m; i++ {
248 a[i*lda+j] = stmp*tmp2 + ctmp*tmp
249 a[i*lda+n-1] = ctmp*tmp2 - stmp*tmp
256 for j := n - 2; j >= 0; j-- {
259 if ctmp != 1 || stmp != 0 {
260 for i := 0; i < m; i++ {
263 a[i*lda+j] = stmp*tmp2 + ctmp*tmp
264 a[i*lda+n-1] = ctmp*tmp2 - stmp*tmp