dlasr.go

// Copyright ©2015 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"
	"gonum.org/v1/gonum/lapack"
)

// Dlasr applies a sequence of plane rotations to the m×n matrix A. This series
// of plane rotations is implicitly represented by a matrix P. P is multiplied
// by a depending on the value of side -- A = P * A if side == lapack.Left,
// A = A * Pᵀ if side == lapack.Right.
//
// The exact value of P depends on the value of pivot, but in all cases P is
// implicitly represented by a series of 2×2 rotation matrices. The entries of
// rotation matrix k are defined by s[k] and c[k]
//
//	R(k) = [ c[k] s[k]]
//	       [-s[k] s[k]]
//
// If direct == lapack.Forward, the rotation matrices are applied as
// P = P(z-1) * ... * P(2) * P(1), while if direct == lapack.Backward they are
// applied as P = P(1) * P(2) * ... * P(n).
//
// pivot defines the mapping of the elements in R(k) to P(k).
// If pivot == lapack.Variable, the rotation is performed for the (k, k+1) plane.
//
//	P(k) = [1                    ]
//	       [    ...              ]
//	       [     1               ]
//	       [       c[k] s[k]     ]
//	       [      -s[k] c[k]     ]
//	       [                 1   ]
//	       [                ...  ]
//	       [                    1]
//
// if pivot == lapack.Top, the rotation is performed for the (1, k+1) plane,
//
//	P(k) = [c[k]        s[k]     ]
//	       [    1                ]
//	       [     ...             ]
//	       [         1           ]
//	       [-s[k]       c[k]     ]
//	       [                 1   ]
//	       [                ...  ]
//	       [                    1]
//
// and if pivot == lapack.Bottom, the rotation is performed for the (k, z) plane.
//
//	P(k) = [1                    ]
//	       [  ...                ]
//	       [      1              ]
//	       [        c[k]     s[k]]
//	       [           1         ]
//	       [            ...      ]
//	       [              1      ]
//	       [       -s[k]     c[k]]
//
// s and c have length m - 1 if side == blas.Left, and n - 1 if side == blas.Right.
//
// Dlasr is an internal routine. It is exported for testing purposes.
func (impl Implementation) Dlasr(side blas.Side, pivot lapack.Pivot, direct lapack.Direct, m, n int, c, s, a []float64, lda int) {
	switch {
	case side != blas.Left && side != blas.Right:
		panic(badSide)
	case pivot != lapack.Variable && pivot != lapack.Top && pivot != lapack.Bottom:
		panic(badPivot)
	case direct != lapack.Forward && direct != lapack.Backward:
		panic(badDirect)
	case m < 0:
		panic(mLT0)
	case n < 0:
		panic(nLT0)
	case lda < max(1, n):
		panic(badLdA)
	}

	// Quick return if possible.
	if m == 0 || n == 0 {
		return
	}

	if side == blas.Left {
		if len(c) < m-1 {
			panic(shortC)
		}
		if len(s) < m-1 {
			panic(shortS)
		}
	} else {
		if len(c) < n-1 {
			panic(shortC)
		}
		if len(s) < n-1 {
			panic(shortS)
		}
	}
	if len(a) < (m-1)*lda+n {
		panic(shortA)
	}

	if side == blas.Left {
		if pivot == lapack.Variable {
			if direct == lapack.Forward {
				for j := 0; j < m-1; j++ {
					ctmp := c[j]
					stmp := s[j]
					if ctmp != 1 || stmp != 0 {
						for i := 0; i < n; i++ {
							tmp2 := a[j*lda+i]
							tmp := a[(j+1)*lda+i]
							a[(j+1)*lda+i] = ctmp*tmp - stmp*tmp2
							a[j*lda+i] = stmp*tmp + ctmp*tmp2
						}
					}
				}
				return
			}
			for j := m - 2; j >= 0; j-- {
				ctmp := c[j]
				stmp := s[j]
				if ctmp != 1 || stmp != 0 {
					for i := 0; i < n; i++ {
						tmp2 := a[j*lda+i]
						tmp := a[(j+1)*lda+i]
						a[(j+1)*lda+i] = ctmp*tmp - stmp*tmp2
						a[j*lda+i] = stmp*tmp + ctmp*tmp2
					}
				}
			}
			return
		} else if pivot == lapack.Top {
			if direct == lapack.Forward {
				for j := 1; j < m; j++ {
					ctmp := c[j-1]
					stmp := s[j-1]
					if ctmp != 1 || stmp != 0 {
						for i := 0; i < n; i++ {
							tmp := a[j*lda+i]
							tmp2 := a[i]
							a[j*lda+i] = ctmp*tmp - stmp*tmp2
							a[i] = stmp*tmp + ctmp*tmp2
						}
					}
				}
				return
			}
			for j := m - 1; j >= 1; j-- {
				ctmp := c[j-1]
				stmp := s[j-1]
				if ctmp != 1 || stmp != 0 {
					for i := 0; i < n; i++ {
						ctmp := c[j-1]
						stmp := s[j-1]
						if ctmp != 1 || stmp != 0 {
							for i := 0; i < n; i++ {
								tmp := a[j*lda+i]
								tmp2 := a[i]
								a[j*lda+i] = ctmp*tmp - stmp*tmp2
								a[i] = stmp*tmp + ctmp*tmp2
							}
						}
					}
				}
			}
			return
		}
		if direct == lapack.Forward {
			for j := 0; j < m-1; j++ {
				ctmp := c[j]
				stmp := s[j]
				if ctmp != 1 || stmp != 0 {
					for i := 0; i < n; i++ {
						tmp := a[j*lda+i]
						tmp2 := a[(m-1)*lda+i]
						a[j*lda+i] = stmp*tmp2 + ctmp*tmp
						a[(m-1)*lda+i] = ctmp*tmp2 - stmp*tmp
					}
				}
			}
			return
		}
		for j := m - 2; j >= 0; j-- {
			ctmp := c[j]
			stmp := s[j]
			if ctmp != 1 || stmp != 0 {
				for i := 0; i < n; i++ {
					tmp := a[j*lda+i]
					tmp2 := a[(m-1)*lda+i]
					a[j*lda+i] = stmp*tmp2 + ctmp*tmp
					a[(m-1)*lda+i] = ctmp*tmp2 - stmp*tmp
				}
			}
		}
		return
	}
	if pivot == lapack.Variable {
		if direct == lapack.Forward {
			for j := 0; j < n-1; j++ {
				ctmp := c[j]
				stmp := s[j]
				if ctmp != 1 || stmp != 0 {
					for i := 0; i < m; i++ {
						tmp := a[i*lda+j+1]
						tmp2 := a[i*lda+j]
						a[i*lda+j+1] = ctmp*tmp - stmp*tmp2
						a[i*lda+j] = stmp*tmp + ctmp*tmp2
					}
				}
			}
			return
		}
		for j := n - 2; j >= 0; j-- {
			ctmp := c[j]
			stmp := s[j]
			if ctmp != 1 || stmp != 0 {
				for i := 0; i < m; i++ {
					tmp := a[i*lda+j+1]
					tmp2 := a[i*lda+j]
					a[i*lda+j+1] = ctmp*tmp - stmp*tmp2
					a[i*lda+j] = stmp*tmp + ctmp*tmp2
				}
			}
		}
		return
	} else if pivot == lapack.Top {
		if direct == lapack.Forward {
			for j := 1; j < n; j++ {
				ctmp := c[j-1]
				stmp := s[j-1]
				if ctmp != 1 || stmp != 0 {
					for i := 0; i < m; i++ {
						tmp := a[i*lda+j]
						tmp2 := a[i*lda]
						a[i*lda+j] = ctmp*tmp - stmp*tmp2
						a[i*lda] = stmp*tmp + ctmp*tmp2
					}
				}
			}
			return
		}
		for j := n - 1; j >= 1; j-- {
			ctmp := c[j-1]
			stmp := s[j-1]
			if ctmp != 1 || stmp != 0 {
				for i := 0; i < m; i++ {
					tmp := a[i*lda+j]
					tmp2 := a[i*lda]
					a[i*lda+j] = ctmp*tmp - stmp*tmp2
					a[i*lda] = stmp*tmp + ctmp*tmp2
				}
			}
		}
		return
	}
	if direct == lapack.Forward {
		for j := 0; j < n-1; j++ {
			ctmp := c[j]
			stmp := s[j]
			if ctmp != 1 || stmp != 0 {
				for i := 0; i < m; i++ {
					tmp := a[i*lda+j]
					tmp2 := a[i*lda+n-1]
					a[i*lda+j] = stmp*tmp2 + ctmp*tmp
					a[i*lda+n-1] = ctmp*tmp2 - stmp*tmp
				}

			}
		}
		return
	}
	for j := n - 2; j >= 0; j-- {
		ctmp := c[j]
		stmp := s[j]
		if ctmp != 1 || stmp != 0 {
			for i := 0; i < m; i++ {
				tmp := a[i*lda+j]
				tmp2 := a[i*lda+n-1]
				a[i*lda+j] = stmp*tmp2 + ctmp*tmp
				a[i*lda+n-1] = ctmp*tmp2 - stmp*tmp
			}
		}
	}
}