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