vendor/gonum.org/v1/gonum/blas/gonum/level1cmplx128.go

   1 // Copyright ©2017 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         "math"
   9
  10         "gonum.org/v1/gonum/internal/asm/c128"
  11 )
  12
  13 // Dzasum returns the sum of the absolute values of the elements of x
  14 //  \sum_i |Re(x[i])| + |Im(x[i])|
  15 // Dzasum returns 0 if incX is negative.
  16 func (Implementation) Dzasum(n int, x []complex128, incX int) float64 {
  17         if n < 0 {
  18                 panic(negativeN)
  19         }
  20         if incX < 1 {
  21                 if incX == 0 {
  22                         panic(zeroIncX)
  23                 }
  24                 return 0
  25         }
  26         var sum float64
  27         if incX == 1 {
  28                 if len(x) < n {
  29                         panic(badX)
  30                 }
  31                 for _, v := range x[:n] {
  32                         sum += dcabs1(v)
  33                 }
  34                 return sum
  35         }
  36         if (n-1)*incX >= len(x) {
  37                 panic(badX)
  38         }
  39         for i := 0; i < n; i++ {
  40                 v := x[i*incX]
  41                 sum += dcabs1(v)
  42         }
  43         return sum
  44 }
  45
  46 // Dznrm2 computes the Euclidean norm of the complex vector x,
  47 //  ‖x‖_2 = sqrt(\sum_i x[i] * conj(x[i])).
  48 // This function returns 0 if incX is negative.
  49 func (Implementation) Dznrm2(n int, x []complex128, incX int) float64 {
  50         if incX < 1 {
  51                 if incX == 0 {
  52                         panic(zeroIncX)
  53                 }
  54                 return 0
  55         }
  56         if n < 1 {
  57                 if n == 0 {
  58                         return 0
  59                 }
  60                 panic(negativeN)
  61         }
  62         if (n-1)*incX >= len(x) {
  63                 panic(badX)
  64         }
  65         var (
  66                 scale float64
  67                 ssq   float64 = 1
  68         )
  69         if incX == 1 {
  70                 for _, v := range x[:n] {
  71                         re, im := math.Abs(real(v)), math.Abs(imag(v))
  72                         if re != 0 {
  73                                 if re > scale {
  74                                         ssq = 1 + ssq*(scale/re)*(scale/re)
  75                                         scale = re
  76                                 } else {
  77                                         ssq += (re / scale) * (re / scale)
  78                                 }
  79                         }
  80                         if im != 0 {
  81                                 if im > scale {
  82                                         ssq = 1 + ssq*(scale/im)*(scale/im)
  83                                         scale = im
  84                                 } else {
  85                                         ssq += (im / scale) * (im / scale)
  86                                 }
  87                         }
  88                 }
  89                 if math.IsInf(scale, 1) {
  90                         return math.Inf(1)
  91                 }
  92                 return scale * math.Sqrt(ssq)
  93         }
  94         for ix := 0; ix < n*incX; ix += incX {
  95                 re, im := math.Abs(real(x[ix])), math.Abs(imag(x[ix]))
  96                 if re != 0 {
  97                         if re > scale {
  98                                 ssq = 1 + ssq*(scale/re)*(scale/re)
  99                                 scale = re
 100                         } else {
 101                                 ssq += (re / scale) * (re / scale)
 102                         }
 103                 }
 104                 if im != 0 {
 105                         if im > scale {
 106                                 ssq = 1 + ssq*(scale/im)*(scale/im)
 107                                 scale = im
 108                         } else {
 109                                 ssq += (im / scale) * (im / scale)
 110                         }
 111                 }
 112         }
 113         if math.IsInf(scale, 1) {
 114                 return math.Inf(1)
 115         }
 116         return scale * math.Sqrt(ssq)
 117 }
 118
 119 // Izamax returns the index of the first element of x having largest |Re(·)|+|Im(·)|.
 120 // Izamax returns -1 if n is 0 or incX is negative.
 121 func (Implementation) Izamax(n int, x []complex128, incX int) int {
 122         if incX < 1 {
 123                 if incX == 0 {
 124                         panic(zeroIncX)
 125                 }
 126                 // Return invalid index.
 127                 return -1
 128         }
 129         if n < 1 {
 130                 if n == 0 {
 131                         // Return invalid index.
 132                         return -1
 133                 }
 134                 panic(negativeN)
 135         }
 136         if len(x) <= (n-1)*incX {
 137                 panic(badX)
 138         }
 139         idx := 0
 140         max := dcabs1(x[0])
 141         if incX == 1 {
 142                 for i, v := range x[1:n] {
 143                         absV := dcabs1(v)
 144                         if absV > max {
 145                                 max = absV
 146                                 idx = i + 1
 147                         }
 148                 }
 149                 return idx
 150         }
 151         ix := incX
 152         for i := 1; i < n; i++ {
 153                 absV := dcabs1(x[ix])
 154                 if absV > max {
 155                         max = absV
 156                         idx = i
 157                 }
 158                 ix += incX
 159         }
 160         return idx
 161 }
 162
 163 // Zaxpy adds alpha times x to y:
 164 //  y[i] += alpha * x[i] for all i
 165 func (Implementation) Zaxpy(n int, alpha complex128, x []complex128, incX int, y []complex128, incY int) {
 166         if incX == 0 {
 167                 panic(zeroIncX)
 168         }
 169         if incY == 0 {
 170                 panic(zeroIncY)
 171         }
 172         if n < 1 {
 173                 if n == 0 {
 174                         return
 175                 }
 176                 panic(negativeN)
 177         }
 178         if (incX > 0 && (n-1)*incX >= len(x)) || (incX < 0 && (1-n)*incX >= len(x)) {
 179                 panic(badX)
 180         }
 181         if (incY > 0 && (n-1)*incY >= len(y)) || (incY < 0 && (1-n)*incY >= len(y)) {
 182                 panic(badY)
 183         }
 184         if alpha == 0 {
 185                 return
 186         }
 187         if incX == 1 && incY == 1 {
 188                 c128.AxpyUnitary(alpha, x[:n], y[:n])
 189                 return
 190         }
 191         var ix, iy int
 192         if incX < 0 {
 193                 ix = (1 - n) * incX
 194         }
 195         if incY < 0 {
 196                 iy = (1 - n) * incY
 197         }
 198         c128.AxpyInc(alpha, x, y, uintptr(n), uintptr(incX), uintptr(incY), uintptr(ix), uintptr(iy))
 199 }
 200
 201 // Zcopy copies the vector x to vector y.
 202 func (Implementation) Zcopy(n int, x []complex128, incX int, y []complex128, incY int) {
 203         if incX == 0 {
 204                 panic(zeroIncX)
 205         }
 206         if incY == 0 {
 207                 panic(zeroIncY)
 208         }
 209         if n < 1 {
 210                 if n == 0 {
 211                         return
 212                 }
 213                 panic(negativeN)
 214         }
 215         if (incX > 0 && (n-1)*incX >= len(x)) || (incX < 0 && (1-n)*incX >= len(x)) {
 216                 panic(badX)
 217         }
 218         if (incY > 0 && (n-1)*incY >= len(y)) || (incY < 0 && (1-n)*incY >= len(y)) {
 219                 panic(badY)
 220         }
 221         if incX == 1 && incY == 1 {
 222                 copy(y[:n], x[:n])
 223                 return
 224         }
 225         var ix, iy int
 226         if incX < 0 {
 227                 ix = (-n + 1) * incX
 228         }
 229         if incY < 0 {
 230                 iy = (-n + 1) * incY
 231         }
 232         for i := 0; i < n; i++ {
 233                 y[iy] = x[ix]
 234                 ix += incX
 235                 iy += incY
 236         }
 237 }
 238
 239 // Zdotc computes the dot product
 240 //  x^H · y
 241 // of two complex vectors x and y.
 242 func (Implementation) Zdotc(n int, x []complex128, incX int, y []complex128, incY int) complex128 {
 243         if incX == 0 {
 244                 panic(zeroIncX)
 245         }
 246         if incY == 0 {
 247                 panic(zeroIncY)
 248         }
 249         if n <= 0 {
 250                 if n == 0 {
 251                         return 0
 252                 }
 253                 panic(negativeN)
 254         }
 255         if incX == 1 && incY == 1 {
 256                 if len(x) < n {
 257                         panic(badX)
 258                 }
 259                 if len(y) < n {
 260                         panic(badY)
 261                 }
 262                 return c128.DotcUnitary(x[:n], y[:n])
 263         }
 264         var ix, iy int
 265         if incX < 0 {
 266                 ix = (-n + 1) * incX
 267         }
 268         if incY < 0 {
 269                 iy = (-n + 1) * incY
 270         }
 271         if ix >= len(x) || (n-1)*incX >= len(x) {
 272                 panic(badX)
 273         }
 274         if iy >= len(y) || (n-1)*incY >= len(y) {
 275                 panic(badY)
 276         }
 277         return c128.DotcInc(x, y, uintptr(n), uintptr(incX), uintptr(incY), uintptr(ix), uintptr(iy))
 278 }
 279
 280 // Zdotu computes the dot product
 281 //  x^T · y
 282 // of two complex vectors x and y.
 283 func (Implementation) Zdotu(n int, x []complex128, incX int, y []complex128, incY int) complex128 {
 284         if incX == 0 {
 285                 panic(zeroIncX)
 286         }
 287         if incY == 0 {
 288                 panic(zeroIncY)
 289         }
 290         if n <= 0 {
 291                 if n == 0 {
 292                         return 0
 293                 }
 294                 panic(negativeN)
 295         }
 296         if incX == 1 && incY == 1 {
 297                 if len(x) < n {
 298                         panic(badX)
 299                 }
 300                 if len(y) < n {
 301                         panic(badY)
 302                 }
 303                 return c128.DotuUnitary(x[:n], y[:n])
 304         }
 305         var ix, iy int
 306         if incX < 0 {
 307                 ix = (-n + 1) * incX
 308         }
 309         if incY < 0 {
 310                 iy = (-n + 1) * incY
 311         }
 312         if ix >= len(x) || (n-1)*incX >= len(x) {
 313                 panic(badX)
 314         }
 315         if iy >= len(y) || (n-1)*incY >= len(y) {
 316                 panic(badY)
 317         }
 318         return c128.DotuInc(x, y, uintptr(n), uintptr(incX), uintptr(incY), uintptr(ix), uintptr(iy))
 319 }
 320
 321 // Zdscal scales the vector x by a real scalar alpha.
 322 // Zdscal has no effect if incX < 0.
 323 func (Implementation) Zdscal(n int, alpha float64, x []complex128, incX int) {
 324         if incX < 1 {
 325                 if incX == 0 {
 326                         panic(zeroIncX)
 327                 }
 328                 return
 329         }
 330         if (n-1)*incX >= len(x) {
 331                 panic(badX)
 332         }
 333         if n < 1 {
 334                 if n == 0 {
 335                         return
 336                 }
 337                 panic(negativeN)
 338         }
 339         if alpha == 0 {
 340                 if incX == 1 {
 341                         x = x[:n]
 342                         for i := range x {
 343                                 x[i] = 0
 344                         }
 345                         return
 346                 }
 347                 for ix := 0; ix < n*incX; ix += incX {
 348                         x[ix] = 0
 349                 }
 350                 return
 351         }
 352         if incX == 1 {
 353                 x = x[:n]
 354                 for i, v := range x {
 355                         x[i] = complex(alpha*real(v), alpha*imag(v))
 356                 }
 357                 return
 358         }
 359         for ix := 0; ix < n*incX; ix += incX {
 360                 v := x[ix]
 361                 x[ix] = complex(alpha*real(v), alpha*imag(v))
 362         }
 363 }
 364
 365 // Zscal scales the vector x by a complex scalar alpha.
 366 // Zscal has no effect if incX < 0.
 367 func (Implementation) Zscal(n int, alpha complex128, x []complex128, incX int) {
 368         if incX < 1 {
 369                 if incX == 0 {
 370                         panic(zeroIncX)
 371                 }
 372                 return
 373         }
 374         if (n-1)*incX >= len(x) {
 375                 panic(badX)
 376         }
 377         if n < 1 {
 378                 if n == 0 {
 379                         return
 380                 }
 381                 panic(negativeN)
 382         }
 383         if alpha == 0 {
 384                 if incX == 1 {
 385                         x = x[:n]
 386                         for i := range x {
 387                                 x[i] = 0
 388                         }
 389                         return
 390                 }
 391                 for ix := 0; ix < n*incX; ix += incX {
 392                         x[ix] = 0
 393                 }
 394                 return
 395         }
 396         if incX == 1 {
 397                 c128.ScalUnitary(alpha, x[:n])
 398                 return
 399         }
 400         c128.ScalInc(alpha, x, uintptr(n), uintptr(incX))
 401 }
 402
 403 // Zswap exchanges the elements of two complex vectors x and y.
 404 func (Implementation) Zswap(n int, x []complex128, incX int, y []complex128, incY int) {
 405         if incX == 0 {
 406                 panic(zeroIncX)
 407         }
 408         if incY == 0 {
 409                 panic(zeroIncY)
 410         }
 411         if n < 1 {
 412                 if n == 0 {
 413                         return
 414                 }
 415                 panic(negativeN)
 416         }
 417         if (incX > 0 && (n-1)*incX >= len(x)) || (incX < 0 && (1-n)*incX >= len(x)) {
 418                 panic(badX)
 419         }
 420         if (incY > 0 && (n-1)*incY >= len(y)) || (incY < 0 && (1-n)*incY >= len(y)) {
 421                 panic(badY)
 422         }
 423         if incX == 1 && incY == 1 {
 424                 x = x[:n]
 425                 for i, v := range x {
 426                         x[i], y[i] = y[i], v
 427                 }
 428                 return
 429         }
 430         var ix, iy int
 431         if incX < 0 {
 432                 ix = (-n + 1) * incX
 433         }
 434         if incY < 0 {
 435                 iy = (-n + 1) * incY
 436         }
 437         for i := 0; i < n; i++ {
 438                 x[ix], y[iy] = y[iy], x[ix]
 439                 ix += incX
 440                 iy += incY
 441         }
 442 }