vendor/gonum.org/v1/gonum/lapack/gonum/dlasr.go

   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
   5 package gonum
   6
   7 import (
   8         "gonum.org/v1/gonum/blas"
   9         "gonum.org/v1/gonum/lapack"
  10 )
  11
  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.
  16 //
  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]]
  21 //         [-s[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).
  25 //
  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.
  28 //  P(k) = [1                    ]
  29 //         [    ...              ]
  30 //         [     1               ]
  31 //         [       c[k] s[k]     ]
  32 //         [      -s[k] c[k]     ]
  33 //         [                 1   ]
  34 //         [                ...  ]
  35 //         [                    1]
  36 // if pivot == lapack.Top, the rotation is performed for the (1, k+1) plane,
  37 //  P(k) = [c[k]        s[k]     ]
  38 //         [    1                ]
  39 //         [     ...             ]
  40 //         [         1           ]
  41 //         [-s[k]       c[k]     ]
  42 //         [                 1   ]
  43 //         [                ...  ]
  44 //         [                    1]
  45 // and if pivot == lapack.Bottom, the rotation is performed for the (k, z) plane.
  46 //  P(k) = [1                    ]
  47 //         [  ...                ]
  48 //         [      1              ]
  49 //         [        c[k]     s[k]]
  50 //         [           1         ]
  51 //         [            ...      ]
  52 //         [              1      ]
  53 //         [       -s[k]     c[k]]
  54 // s and c have length m - 1 if side == blas.Left, and n - 1 if side == blas.Right.
  55 //
  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 {
  60                 panic(badSide)
  61         }
  62         if pivot != lapack.Variable && pivot != lapack.Top && pivot != lapack.Bottom {
  63                 panic(badPivot)
  64         }
  65         if direct != lapack.Forward && direct != lapack.Backward {
  66                 panic(badDirect)
  67         }
  68         if side == blas.Left {
  69                 if len(c) < m-1 {
  70                         panic(badSlice)
  71                 }
  72                 if len(s) < m-1 {
  73                         panic(badSlice)
  74                 }
  75         } else {
  76                 if len(c) < n-1 {
  77                         panic(badSlice)
  78                 }
  79                 if len(s) < n-1 {
  80                         panic(badSlice)
  81                 }
  82         }
  83         if m == 0 || n == 0 {
  84                 return
  85         }
  86         if side == blas.Left {
  87                 if pivot == lapack.Variable {
  88                         if direct == lapack.Forward {
  89                                 for j := 0; j < m-1; j++ {
  90                                         ctmp := c[j]
  91                                         stmp := s[j]
  92                                         if ctmp != 1 || stmp != 0 {
  93                                                 for i := 0; i < n; i++ {
  94                                                         tmp2 := a[j*lda+i]
  95                                                         tmp := a[(j+1)*lda+i]
  96                                                         a[(j+1)*lda+i] = ctmp*tmp - stmp*tmp2
  97                                                         a[j*lda+i] = stmp*tmp + ctmp*tmp2
  98                                                 }
  99                                         }
 100                                 }
 101                                 return
 102                         }
 103                         for j := m - 2; j >= 0; j-- {
 104                                 ctmp := c[j]
 105                                 stmp := s[j]
 106                                 if ctmp != 1 || stmp != 0 {
 107                                         for i := 0; i < n; i++ {
 108                                                 tmp2 := a[j*lda+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
 112                                         }
 113                                 }
 114                         }
 115                         return
 116                 } else if pivot == lapack.Top {
 117                         if direct == lapack.Forward {
 118                                 for j := 1; j < m; j++ {
 119                                         ctmp := c[j-1]
 120                                         stmp := s[j-1]
 121                                         if ctmp != 1 || stmp != 0 {
 122                                                 for i := 0; i < n; i++ {
 123                                                         tmp := a[j*lda+i]
 124                                                         tmp2 := a[i]
 125                                                         a[j*lda+i] = ctmp*tmp - stmp*tmp2
 126                                                         a[i] = stmp*tmp + ctmp*tmp2
 127                                                 }
 128                                         }
 129                                 }
 130                                 return
 131                         }
 132                         for j := m - 1; j >= 1; j-- {
 133                                 ctmp := c[j-1]
 134                                 stmp := s[j-1]
 135                                 if ctmp != 1 || stmp != 0 {
 136                                         for i := 0; i < n; i++ {
 137                                                 ctmp := c[j-1]
 138                                                 stmp := s[j-1]
 139                                                 if ctmp != 1 || stmp != 0 {
 140                                                         for i := 0; i < n; i++ {
 141                                                                 tmp := a[j*lda+i]
 142                                                                 tmp2 := a[i]
 143                                                                 a[j*lda+i] = ctmp*tmp - stmp*tmp2
 144                                                                 a[i] = stmp*tmp + ctmp*tmp2
 145                                                         }
 146                                                 }
 147                                         }
 148                                 }
 149                         }
 150                         return
 151                 }
 152                 if direct == lapack.Forward {
 153                         for j := 0; j < m-1; j++ {
 154                                 ctmp := c[j]
 155                                 stmp := s[j]
 156                                 if ctmp != 1 || stmp != 0 {
 157                                         for i := 0; i < n; i++ {
 158                                                 tmp := a[j*lda+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
 162                                         }
 163                                 }
 164                         }
 165                         return
 166                 }
 167                 for j := m - 2; j >= 0; j-- {
 168                         ctmp := c[j]
 169                         stmp := s[j]
 170                         if ctmp != 1 || stmp != 0 {
 171                                 for i := 0; i < n; i++ {
 172                                         tmp := a[j*lda+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
 176                                 }
 177                         }
 178                 }
 179                 return
 180         }
 181         if pivot == lapack.Variable {
 182                 if direct == lapack.Forward {
 183                         for j := 0; j < n-1; j++ {
 184                                 ctmp := c[j]
 185                                 stmp := s[j]
 186                                 if ctmp != 1 || stmp != 0 {
 187                                         for i := 0; i < m; i++ {
 188                                                 tmp := a[i*lda+j+1]
 189                                                 tmp2 := a[i*lda+j]
 190                                                 a[i*lda+j+1] = ctmp*tmp - stmp*tmp2
 191                                                 a[i*lda+j] = stmp*tmp + ctmp*tmp2
 192                                         }
 193                                 }
 194                         }
 195                         return
 196                 }
 197                 for j := n - 2; j >= 0; j-- {
 198                         ctmp := c[j]
 199                         stmp := s[j]
 200                         if ctmp != 1 || stmp != 0 {
 201                                 for i := 0; i < m; i++ {
 202                                         tmp := a[i*lda+j+1]
 203                                         tmp2 := a[i*lda+j]
 204                                         a[i*lda+j+1] = ctmp*tmp - stmp*tmp2
 205                                         a[i*lda+j] = stmp*tmp + ctmp*tmp2
 206                                 }
 207                         }
 208                 }
 209                 return
 210         } else if pivot == lapack.Top {
 211                 if direct == lapack.Forward {
 212                         for j := 1; j < n; j++ {
 213                                 ctmp := c[j-1]
 214                                 stmp := s[j-1]
 215                                 if ctmp != 1 || stmp != 0 {
 216                                         for i := 0; i < m; i++ {
 217                                                 tmp := a[i*lda+j]
 218                                                 tmp2 := a[i*lda]
 219                                                 a[i*lda+j] = ctmp*tmp - stmp*tmp2
 220                                                 a[i*lda] = stmp*tmp + ctmp*tmp2
 221                                         }
 222                                 }
 223                         }
 224                         return
 225                 }
 226                 for j := n - 1; j >= 1; j-- {
 227                         ctmp := c[j-1]
 228                         stmp := s[j-1]
 229                         if ctmp != 1 || stmp != 0 {
 230                                 for i := 0; i < m; i++ {
 231                                         tmp := a[i*lda+j]
 232                                         tmp2 := a[i*lda]
 233                                         a[i*lda+j] = ctmp*tmp - stmp*tmp2
 234                                         a[i*lda] = stmp*tmp + ctmp*tmp2
 235                                 }
 236                         }
 237                 }
 238                 return
 239         }
 240         if direct == lapack.Forward {
 241                 for j := 0; j < n-1; j++ {
 242                         ctmp := c[j]
 243                         stmp := s[j]
 244                         if ctmp != 1 || stmp != 0 {
 245                                 for i := 0; i < m; i++ {
 246                                         tmp := a[i*lda+j]
 247                                         tmp2 := a[i*lda+n-1]
 248                                         a[i*lda+j] = stmp*tmp2 + ctmp*tmp
 249                                         a[i*lda+n-1] = ctmp*tmp2 - stmp*tmp
 250                                 }
 251
 252                         }
 253                 }
 254                 return
 255         }
 256         for j := n - 2; j >= 0; j-- {
 257                 ctmp := c[j]
 258                 stmp := s[j]
 259                 if ctmp != 1 || stmp != 0 {
 260                         for i := 0; i < m; i++ {
 261                                 tmp := a[i*lda+j]
 262                                 tmp2 := a[i*lda+n-1]
 263                                 a[i*lda+j] = stmp*tmp2 + ctmp*tmp
 264                                 a[i*lda+n-1] = ctmp*tmp2 - stmp*tmp
 265                         }
 266                 }
 267         }
 268 }