level1cmplx128.go 8.1 KB

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