/[escript]/branches/diaplayground/cusplibrary/cusp/complex.h
ViewVC logotype

Contents of /branches/diaplayground/cusplibrary/cusp/complex.h

Parent Directory Parent Directory | Revision Log Revision Log


Revision 4956 - (show annotations)
Tue May 20 12:34:41 2014 UTC (5 years, 1 month ago) by caltinay
File MIME type: text/plain
File size: 35441 byte(s)
minor changes to cusp to allow compilation with -pedantic and -Werror.
One bug fix in spmv_dia.

1 /*
2 * Copyright 2008-2009 NVIDIA Corporation
3 *
4 * Licensed under the Apache License, Version 2.0 (the "License");
5 * you may not use this file except in compliance with the License.
6 * You may obtain a copy of the License at
7 *
8 * http://www.apache.org/licenses/LICENSE-2.0
9 *
10 * Unless required by applicable law or agreed to in writing, software
11 * distributed under the License is distributed on an "AS IS" BASIS,
12 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
13 * See the License for the specific language governing permissions and
14 * limitations under the License.
15 */
16
17 /*
18 * Copyright (c) 2010, The Regents of the University of California,
19 * through Lawrence Berkeley National Laboratory (subject to receipt of
20 * any required approvals from U.S. Dept. of Energy) All rights reserved.
21 *
22 * Redistribution and use in source and binary forms, with or
23 * without modification, are permitted provided that the
24 * following conditions are met:
25 *
26 * * Redistributions of source code must retain the above
27 * copyright notice, this list of conditions and the following
28 * disclaimer.
29 *
30 * * Redistributions in binary form must reproduce the
31 * above copyright notice, this list of conditions and the
32 * following disclaimer in the documentation and/or other
33 * materials provided with the distribution.
34 *
35 * * Neither the name of the University of California,
36 * Berkeley, nor the names of its contributors may be used to
37 * endorse or promote products derived from this software
38 * without specific prior written permission.
39 *
40 * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND
41 * CONTRIBUTORS "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES,
42 * INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF
43 * MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
44 * DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR
45 * CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
46 * SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING,
47 * BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
48 * SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
49 * INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY,
50 * WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
51 * NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
52 * OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
53 * SUCH DAMAGE.
54 *
55 */
56
57 /*! \file complex.h
58 * \brief Complex numbers
59 */
60
61 #pragma once
62
63 #include <cusp/detail/config.h>
64
65 #if (defined THRUST_DEVICE_BACKEND && THRUST_DEVICE_BACKEND == THRUST_DEVICE_BACKEND_CUDA) || (defined THRUST_DEVICE_SYSTEM && THRUST_DEVICE_SYSTEM == THRUST_DEVICE_SYSTEM_CUDA)
66
67 #ifdef _WIN32
68 #define _USE_MATH_DEFINES 1 // make sure M_PI is defined
69 #endif
70
71 #include <math.h>
72 #include <complex>
73 #include <cuComplex.h>
74 #include <sstream>
75 #include <cusp/cmath.h>
76
77 namespace cusp
78 {
79
80 template <typename ValueType> struct complex;
81 template <> struct complex<float>;
82 template <> struct complex<double>;
83
84
85 /// Returns the magnitude of z.
86 template<typename ValueType>
87 __host__ __device__
88 ValueType abs(const complex<ValueType>& z);
89 /// Returns the phase angle of z.
90 template<typename ValueType>
91 __host__ __device__
92 ValueType arg(const complex<ValueType>& z);
93 /// Returns the magnitude of z squared.
94 template<typename ValueType>
95 __host__ __device__
96 ValueType norm(const complex<ValueType>& z);
97
98 /// Returns the complex conjugate of z.
99 template<typename ValueType>
100 __host__ __device__
101 complex<ValueType> conj(const complex<ValueType>& z);
102 /// Returns the complex with magnitude m and angle theta in radians.
103 template<typename ValueType>
104 __host__ __device__
105 complex<ValueType> polar(const ValueType& m, const ValueType& theta = 0);
106
107 // Arithmetic operators:
108 // Multiplication
109 template <typename ValueType>
110 __host__ __device__
111 inline complex<ValueType> operator*(const complex<ValueType>& lhs, const complex<ValueType>& rhs);
112 template <typename ValueType>
113 __host__ __device__
114 inline complex<ValueType> operator*(const complex<ValueType>& lhs, const ValueType & rhs);
115 template <typename ValueType>
116 __host__ __device__
117 inline complex<ValueType> operator*(const ValueType& lhs, const complex<ValueType>& rhs);
118 // Division
119 template <typename ValueType>
120 __host__ __device__
121 inline complex<ValueType> operator/(const complex<ValueType>& lhs, const complex<ValueType>& rhs);
122 template <>
123 __host__ __device__
124 inline complex<float> operator/(const complex<float>& lhs, const complex<float>& rhs);
125 template <>
126 __host__ __device__
127 inline complex<double> operator/(const complex<double>& lhs, const complex<double>& rhs);
128
129 // Addition
130 template <typename ValueType>
131 __host__ __device__
132 inline complex<ValueType> operator+(const complex<ValueType>& lhs, const complex<ValueType>& rhs);
133 template <typename ValueType>
134 __host__ __device__
135 inline complex<ValueType> operator+(const complex<ValueType>& lhs, const ValueType & rhs);
136 template <typename ValueType>
137 __host__ __device__
138 inline complex<ValueType> operator+(const ValueType& lhs, const complex<ValueType>& rhs);
139 // Subtraction
140 template <typename ValueType>
141 __host__ __device__
142 inline complex<ValueType> operator-(const complex<ValueType>& lhs, const complex<ValueType>& rhs);
143 template <typename ValueType>
144 __host__ __device__
145 inline complex<ValueType> operator-(const complex<ValueType>& lhs, const ValueType & rhs);
146 template <typename ValueType>
147 __host__ __device__
148 inline complex<ValueType> operator-(const ValueType& lhs, const complex<ValueType>& rhs);
149
150 // Unary plus and minus
151 template <typename ValueType>
152 __host__ __device__
153 inline complex<ValueType> operator+(const complex<ValueType>& rhs);
154 template <typename ValueType>
155 __host__ __device__
156 inline complex<ValueType> operator-(const complex<ValueType>& rhs);
157
158 // Transcendentals:
159 // Returns the complex cosine of z.
160 template <typename ValueType>
161 __host__ __device__
162 complex<ValueType> cos(const complex<ValueType>& z);
163 // Returns the complex hyperbolic cosine of z.
164 template <typename ValueType>
165 __host__ __device__
166 complex<ValueType> cosh(const complex<ValueType>& z);
167 // Returns the complex base e exponential of z.
168 template <typename ValueType>
169 __host__ __device__
170 complex<ValueType> exp(const complex<ValueType>& z);
171 // Returns the complex natural logarithm of z.
172 template <typename ValueType>
173 __host__ __device__
174 complex<ValueType> log(const complex<ValueType>& z);
175 // Returns the complex base 10 logarithm of z.
176 template <typename ValueType>
177 __host__ __device__
178 complex<ValueType> log10(const complex<ValueType>& z);
179 // Returns z to the n'th power.
180 template <typename ValueType>
181 __host__ __device__
182 complex<ValueType> pow(const complex<ValueType>& z, const int& n);
183 // Returns z to the x'th power.
184 template <typename ValueType>
185 __host__ __device__
186 complex<ValueType> pow(const complex<ValueType>&z, const ValueType&x);
187 // Returns z to the z2'th power.
188 template <typename ValueType>
189 __host__ __device__
190 complex<ValueType> pow(const complex<ValueType>&z,
191 const complex<ValueType>&z2);
192 // Returns x to the z'th power.
193 template <typename ValueType>
194 __host__ __device__
195 complex<ValueType> pow(const ValueType& x, const complex<ValueType>& z);
196 // Returns the complex sine of z.
197 template <typename ValueType>
198 __host__ __device__
199 complex<ValueType> sin(const complex<ValueType>&z);
200 // Returns the complex hyperbolic sine of z.
201 template <typename ValueType>
202 __host__ __device__
203 complex<ValueType> sinh(const complex<ValueType>&z);
204 // Returns the complex square root of z.
205 template <typename ValueType>
206 __host__ __device__
207 complex<ValueType> sqrt(const complex<ValueType>&z);
208 // Returns the complex tangent of z.
209 template <typename ValueType>
210 __host__ __device__
211 complex<ValueType> tan(const complex<ValueType>&z);
212 // Returns the complex hyperbolic tangent of z.
213 template <typename ValueType>
214 __host__ __device__
215 complex<ValueType> tanh(const complex<ValueType>&z);
216
217
218 // Inverse Trigonometric:
219 // Returns the complex arc cosine of z.
220 template <typename ValueType>
221 __host__ __device__
222 complex<ValueType> acos(const complex<ValueType>& z);
223 // Returns the complex arc sine of z.
224 template <typename ValueType>
225 __host__ __device__
226 complex<ValueType> asin(const complex<ValueType>& z);
227 // Returns the complex arc tangent of z.
228 template <typename ValueType>
229 __host__ __device__
230 complex<ValueType> atan(const complex<ValueType>& z);
231 // Returns the complex hyperbolic arc cosine of z.
232 template <typename ValueType>
233 __host__ __device__
234 complex<ValueType> acosh(const complex<ValueType>& z);
235 // Returns the complex hyperbolic arc sine of z.
236 template <typename ValueType>
237 __host__ __device__
238 complex<ValueType> asinh(const complex<ValueType>& z);
239 // Returns the complex hyperbolic arc tangent of z.
240 template <typename ValueType>
241 __host__ __device__
242 complex<ValueType> atanh(const complex<ValueType>& z);
243
244
245
246 // Stream operators:
247 template<typename ValueType,class charT, class traits>
248 std::basic_ostream<charT, traits>& operator<<(std::basic_ostream<charT, traits>& os, const complex<ValueType>& z);
249 template<typename ValueType, typename charT, class traits>
250 std::basic_istream<charT, traits>&
251 operator>>(std::basic_istream<charT, traits>& is, complex<ValueType>& z);
252
253
254 // Stream operators
255 template<typename ValueType,class charT, class traits>
256 std::basic_ostream<charT, traits>& operator<<(std::basic_ostream<charT, traits>& os, const complex<ValueType>& z)
257 {
258 os << '(' << z.real() << ',' << z.imag() << ')';
259 return os;
260 }
261
262 template<typename ValueType, typename charT, class traits>
263 std::basic_istream<charT, traits>&
264 operator>>(std::basic_istream<charT, traits>& is, complex<ValueType>& z)
265 {
266 ValueType re, im;
267
268 charT ch;
269 is >> ch;
270
271 if(ch == '(')
272 {
273 is >> re >> ch;
274 if (ch == ',')
275 {
276 is >> im >> ch;
277 if (ch == ')')
278 {
279 z = complex<ValueType>(re, im);
280 }
281 else
282 {
283 is.setstate(std::ios_base::failbit);
284 }
285 }
286 else if (ch == ')')
287 {
288 z = re;
289 }
290 else
291 {
292 is.setstate(std::ios_base::failbit);
293 }
294 }
295 else
296 {
297 is.putback(ch);
298 is >> re;
299 z = re;
300 }
301 return is;
302 }
303
304 template <typename T>
305 struct norm_type {
306 typedef T type;
307 };
308 template <typename T>
309 struct norm_type< complex<T> > {
310 typedef T type;
311 };
312
313 template <typename ValueType>
314 struct complex
315 {
316 public:
317 typedef ValueType value_type;
318
319 // Constructors
320 __host__ __device__
321 inline complex<ValueType>(const ValueType & re = ValueType(), const ValueType& im = ValueType())
322 {
323 real(re);
324 imag(im);
325 }
326
327 template <class X>
328 __host__ __device__
329 inline complex<ValueType>(const complex<X> & z)
330 {
331 real(z.real());
332 imag(z.imag());
333 }
334
335 template <class X>
336 __host__ __device__
337 inline complex<ValueType>(const std::complex<X> & z)
338 {
339 real(z.real());
340 imag(z.imag());
341 }
342
343 template <typename T>
344 __host__ __device__
345 inline complex<ValueType>& operator=(const complex<T> z)
346 {
347 real(z.real());
348 imag(z.imag());
349 return *this;
350 }
351
352 __host__ __device__
353 inline complex<ValueType>& operator+=(const complex<ValueType> z)
354 {
355 real(real()+z.real());
356 imag(imag()+z.imag());
357 return *this;
358 }
359
360 __host__ __device__
361 inline complex<ValueType>& operator-=(const complex<ValueType> z)
362 {
363 real(real()-z.real());
364 imag(imag()-z.imag());
365 return *this;
366 }
367
368 __host__ __device__
369 inline complex<ValueType>& operator*=(const complex<ValueType> z)
370 {
371 *this = *this * z;
372 return *this;
373 }
374
375 __host__ __device__
376 inline complex<ValueType>& operator/=(const complex<ValueType> z)
377 {
378 *this = *this / z;
379 return *this;
380 }
381
382 __host__ __device__ inline ValueType real() const volatile;
383 __host__ __device__ inline ValueType imag() const volatile;
384 __host__ __device__ inline ValueType real() const;
385 __host__ __device__ inline ValueType imag() const;
386 __host__ __device__ inline void real(ValueType) volatile;
387 __host__ __device__ inline void imag(ValueType) volatile;
388 __host__ __device__ inline void real(ValueType);
389 __host__ __device__ inline void imag(ValueType);
390 };
391
392 // TODO make cuFloatComplex and cuDoubleComplex protected
393 // TODO see if returning references is a perf hazard
394
395 template<>
396 struct complex <float> : public cuFloatComplex
397 {
398 public:
399 typedef float value_type;
400 __host__ __device__
401 inline complex<float>(){};
402 __host__ __device__
403 inline complex<float>(const float & re, const float& im = float())
404 {
405 real(re);
406 imag(im);
407 }
408
409 // For some reason having the following constructor
410 // explicitly makes things faster with at least g++
411 __host__ __device__
412 complex<float>(const complex<float> & z)
413 : cuFloatComplex(z){}
414
415 __host__ __device__
416 complex<float>(cuFloatComplex z)
417 : cuFloatComplex(z){}
418
419 template <class X>
420 inline complex<float>(const std::complex<X> & z)
421 {
422 real(z.real());
423 imag(z.imag());
424 }
425
426 // Member operators
427 template <typename T>
428 __host__ __device__
429 inline volatile complex<float>& operator=(const complex<T> z) volatile
430 {
431 real(z.real());
432 imag(z.imag());
433 return *this;
434 }
435
436 template <typename T>
437 __host__ __device__
438 inline complex<float>& operator=(const complex<T> z)
439 {
440 real(z.real());
441 imag(z.imag());
442 return *this;
443 }
444
445 __host__ __device__
446 inline complex<float>& operator+=(const complex<float> z)
447 {
448 real(real()+z.real());
449 imag(imag()+z.imag());
450 return *this;
451 }
452
453 __host__ __device__
454 inline complex<float>& operator-=(const complex<float> z)
455 {
456 real(real()-z.real());
457 imag(imag()-z.imag());
458 return *this;
459 }
460
461 __host__ __device__
462 inline complex<float>& operator*=(const complex<float> z)
463 {
464 *this = *this * z;
465 return *this;
466 }
467
468 __host__ __device__
469 inline complex<float>& operator/=(const complex<float> z)
470 {
471 *this = *this / z;
472 return *this;
473 }
474
475 // Let the compiler synthesize the copy and assignment operators.
476 __host__ __device__ inline complex<float>(const volatile complex<float> & z)
477 {
478 real(z.real());
479 imag(z.imag());
480 }
481
482 __host__ __device__ inline float real() const volatile{ return x; }
483 __host__ __device__ inline float imag() const volatile{ return y; }
484 __host__ __device__ inline float real() const{ return x; }
485 __host__ __device__ inline float imag() const{ return y; }
486 __host__ __device__ inline void real(float re)volatile{ x = re; }
487 __host__ __device__ inline void imag(float im)volatile{ y = im; }
488 __host__ __device__ inline void real(float re){ x = re; }
489 __host__ __device__ inline void imag(float im){ y = im; }
490
491 // cast operators
492 inline operator std::complex<float>() const { return std::complex<float>(real(),imag()); }
493 // inline operator float() const { return real(); }
494 };
495
496 template<>
497 struct complex <double> : public cuDoubleComplex
498 {
499 public:
500 typedef double value_type;
501 __host__ __device__
502 inline complex<double>(){};
503 __host__ __device__
504 inline complex<double>(const double & re, const double& im = double())
505 {
506 real(re);
507 imag(im);
508 }
509
510 // For some reason having the following constructor
511 // explicitly makes things faster with at least g++
512 __host__ __device__
513 inline complex<double>(const complex<double> & z)
514 : cuDoubleComplex(z) {}
515
516 __host__ __device__
517 inline complex<double>(cuDoubleComplex z)
518 : cuDoubleComplex(z) {}
519
520 template <class X>
521 inline complex<double>(const std::complex<X> & z)
522 {
523 real(z.real());
524 imag(z.imag());
525 }
526
527 // Member operators
528 template <typename T>
529 __host__ __device__
530 inline volatile complex<double>& operator=(const complex<T> z) volatile
531 {
532 real(z.real());
533 imag(z.imag());
534 return *this;
535 }
536
537 template <typename T>
538 __host__ __device__
539 inline complex<double>& operator=(const complex<T> z)
540 {
541 real(z.real());
542 imag(z.imag());
543 return *this;
544 }
545
546 __host__ __device__
547 inline complex<double>& operator+=(const complex<double> z)
548 {
549 real(real()+z.real());
550 imag(imag()+z.imag());
551 return *this;
552 }
553
554 __host__ __device__
555 inline complex<double>& operator-=(const complex<double> z)
556 {
557 real(real()-z.real());
558 imag(imag()-z.imag());
559 return *this;
560 }
561
562 __host__ __device__
563 inline complex<double>& operator*=(const complex<double> z)
564 {
565 *this = *this * z;
566 return *this;
567 }
568
569 __host__ __device__
570 inline complex<double>& operator/=(const complex<double> z)
571 {
572 *this = *this / z;
573 return *this;
574 }
575
576 __host__ __device__ inline complex<double>(const volatile complex<double> & z)
577 {
578 real(z.real());
579 imag(z.imag());
580 }
581
582 // Let the compiler synthesize the copy and assignment operators.
583 __host__ __device__ inline double real() const volatile { return x; }
584 __host__ __device__ inline double imag() const volatile { return y; }
585 __host__ __device__ inline double real() const { return x; }
586 __host__ __device__ inline double imag() const { return y; }
587 __host__ __device__ inline void real(double re)volatile{ x = re; }
588 __host__ __device__ inline void imag(double im)volatile{ y = im; }
589 __host__ __device__ inline void real(double re){ x = re; }
590 __host__ __device__ inline void imag(double im){ y = im; }
591
592 // cast operators
593 inline operator std::complex<double>() const { return std::complex<double>(real(),imag()); }
594 // inline operator double() { return real(); }
595 };
596
597
598
599 // Binary arithmetic operations
600 // At the moment I'm implementing the basic functions, and the
601 // corresponding cuComplex calls are commented.
602
603 template<typename ValueType>
604 __host__ __device__
605 inline complex<ValueType> operator+(const complex<ValueType>& lhs,
606 const complex<ValueType>& rhs){
607 return complex<ValueType>(lhs.real()+rhs.real(),lhs.imag()+rhs.imag());
608 // return cuCaddf(lhs,rhs);
609 }
610
611 template<typename ValueType>
612 __host__ __device__
613 inline complex<ValueType> operator+(const volatile complex<ValueType>& lhs,
614 const volatile complex<ValueType>& rhs){
615 return complex<ValueType>(lhs.real()+rhs.real(),lhs.imag()+rhs.imag());
616 // return cuCaddf(lhs,rhs);
617 }
618
619 template <typename ValueType>
620 __host__ __device__
621 inline complex<ValueType> operator+(const complex<ValueType>& lhs, const ValueType & rhs){
622 return complex<ValueType>(lhs.real()+rhs,lhs.imag());
623 // return cuCaddf(lhs,complex<ValueType>(rhs));
624 }
625 template <typename ValueType>
626 __host__ __device__
627 inline complex<ValueType> operator+(const ValueType& lhs, const complex<ValueType>& rhs){
628 return complex<ValueType>(rhs.real()+lhs,rhs.imag());
629 // return cuCaddf(complex<float>(lhs),rhs);
630 }
631
632 template <typename ValueType>
633 __host__ __device__
634 inline complex<ValueType> operator-(const complex<ValueType>& lhs, const complex<ValueType>& rhs){
635 return complex<ValueType>(lhs.real()-rhs.real(),lhs.imag()-rhs.imag());
636 // return cuCsubf(lhs,rhs);
637 }
638 template <typename ValueType>
639 __host__ __device__
640 inline complex<ValueType> operator-(const complex<ValueType>& lhs, const ValueType & rhs){
641 return complex<ValueType>(lhs.real()-rhs,lhs.imag());
642 // return cuCsubf(lhs,complex<float>(rhs));
643 }
644 template <typename ValueType>
645 __host__ __device__
646 inline complex<ValueType> operator-(const ValueType& lhs, const complex<ValueType>& rhs){
647 return complex<ValueType>(lhs-rhs.real(),-rhs.imag());
648 // return cuCsubf(complex<float>(lhs),rhs);
649 }
650
651 template <typename ValueType>
652 __host__ __device__
653 inline complex<ValueType> operator*(const complex<ValueType>& lhs,
654 const complex<ValueType>& rhs){
655 return complex<ValueType>(lhs.real()*rhs.real()-lhs.imag()*rhs.imag(),
656 lhs.real()*rhs.imag()+lhs.imag()*rhs.real());
657 // return cuCmulf(lhs,rhs);
658 }
659
660 template <typename ValueType>
661 __host__ __device__
662 inline complex<ValueType> operator*(const complex<ValueType>& lhs, const ValueType & rhs){
663 return complex<ValueType>(lhs.real()*rhs,lhs.imag()*rhs);
664 // return cuCmulf(lhs,complex<float>(rhs));
665 }
666
667 template <typename ValueType>
668 __host__ __device__
669 inline complex<ValueType> operator*(const ValueType& lhs, const complex<ValueType>& rhs){
670 return complex<ValueType>(rhs.real()*lhs,rhs.imag()*lhs);
671 // return cuCmulf(complex<float>(lhs),rhs);
672 }
673
674
675 template <typename ValueType>
676 __host__ __device__
677 inline complex<ValueType> operator/(const complex<ValueType>& lhs, const complex<ValueType>& rhs){
678 const ValueType cross_norm = lhs.real() * rhs.real() + lhs.imag() * rhs.imag();
679 const ValueType rhs_norm = norm(rhs);
680 return complex<ValueType>(cross_norm/rhs_norm,
681 (lhs.imag() * rhs.real() - lhs.real() * rhs.imag()) / rhs_norm);
682 }
683
684 template <>
685 __host__ __device__
686 inline complex<float> operator/(const complex<float>& lhs, const complex<float>& rhs){
687 return cuCdivf(lhs,rhs);
688 }
689
690 template <>
691 __host__ __device__
692 inline complex<double> operator/(const complex<double>& lhs, const complex<double>& rhs){
693 return cuCdiv(lhs,rhs);
694 }
695
696 template <typename ValueType>
697 __host__ __device__
698 inline complex<ValueType> operator/(const complex<ValueType>& lhs, const ValueType & rhs){
699 return complex<ValueType>(lhs.real()/rhs,lhs.imag()/rhs);
700 // return cuCdivf(lhs,complex<float>(rhs));
701 }
702
703 template <typename ValueType>
704 __host__ __device__
705 inline complex<ValueType> operator/(const ValueType& lhs, const complex<ValueType>& rhs){
706 const ValueType cross_norm = lhs * rhs.real();
707 const ValueType rhs_norm = norm(rhs);
708 return complex<ValueType>(cross_norm/rhs_norm,(-lhs.real() * rhs.imag()) / rhs_norm);
709 }
710
711 template <>
712 __host__ __device__
713 inline complex<float> operator/(const float& lhs, const complex<float>& rhs){
714 return cuCdivf(complex<float>(lhs),rhs);
715 }
716 template <>
717 __host__ __device__
718 inline complex<double> operator/(const double& lhs, const complex<double>& rhs){
719 return cuCdiv(complex<double>(lhs),rhs);
720 }
721
722
723 // Unary arithmetic operations
724 template <typename ValueType>
725 __host__ __device__
726 inline complex<ValueType> operator+(const complex<ValueType>& rhs){
727 return rhs;
728 }
729 template <typename ValueType>
730 __host__ __device__
731 inline complex<ValueType> operator-(const complex<ValueType>& rhs){
732 return rhs*-ValueType(1);
733 }
734
735 // Equality operators
736 template <typename ValueType>
737 __host__ __device__
738 inline bool operator==(const complex<ValueType>& lhs, const complex<ValueType>& rhs){
739 if(lhs.real() == rhs.real() && lhs.imag() == rhs.imag()){
740 return true;
741 }
742 return false;
743 }
744 template <typename ValueType>
745 __host__ __device__
746 inline bool operator==(const ValueType & lhs, const complex<ValueType>& rhs){
747 if(lhs == rhs.real() && rhs.imag() == 0){
748 return true;
749 }
750 return false;
751 }
752 template <typename ValueType>
753 __host__ __device__
754 inline bool operator==(const complex<ValueType> & lhs, const ValueType& rhs){
755 if(lhs.real() == rhs && lhs.imag() == 0){
756 return true;
757 }
758 return false;
759 }
760
761 template <typename ValueType>
762 __host__ __device__
763 inline bool operator!=(const complex<ValueType>& lhs, const complex<ValueType>& rhs){
764 return !(lhs == rhs);
765 }
766
767 template <typename ValueType>
768 __host__ __device__
769 inline bool operator!=(const ValueType & lhs, const complex<ValueType>& rhs){
770 return !(lhs == rhs);
771 }
772
773 template <typename ValueType>
774 __host__ __device__
775 inline bool operator!=(const complex<ValueType> & lhs, const ValueType& rhs){
776 return !(lhs == rhs);
777 }
778
779
780 template <typename ValueType>
781 __host__ __device__
782 inline complex<ValueType> conj(const complex<ValueType>& z){
783 return complex<ValueType>(z.real(),-z.imag());
784 }
785
786 template <typename ValueType>
787 __host__ __device__
788 inline ValueType abs(const complex<ValueType>& z){
789 return ::hypot(z.real(),z.imag());
790 }
791 template <>
792 __host__ __device__
793 inline float abs(const complex<float>& z){
794 return ::hypotf(z.real(),z.imag());
795 }
796 template<>
797 __host__ __device__
798 inline double abs(const complex<double>& z){
799 return ::hypot(z.real(),z.imag());
800 }
801
802 template <typename ValueType>
803 __host__ __device__
804 inline ValueType arg(const complex<ValueType>& z){
805 return atan2(z.imag(),z.real());
806 }
807 template<>
808 __host__ __device__
809 inline float arg(const complex<float>& z){
810 return atan2f(z.imag(),z.real());
811 }
812 template<>
813 __host__ __device__
814 inline double arg(const complex<double>& z){
815 return atan2(z.imag(),z.real());
816 }
817
818 template <typename ValueType>
819 __host__ __device__
820 inline ValueType norm(const complex<ValueType>& z){
821 return z.real()*z.real() + z.imag()*z.imag();
822 }
823
824 template <typename ValueType>
825 __host__ __device__
826 inline complex<ValueType> polar(const ValueType & m, const ValueType & theta){
827 return complex<ValueType>(m * ::cos(theta),m * ::sin(theta));
828 }
829
830 template <>
831 __host__ __device__
832 inline complex<float> polar(const float & magnitude, const float & angle){
833 return complex<float>(magnitude * ::cosf(angle),magnitude * ::sinf(angle));
834 }
835
836 template <>
837 __host__ __device__
838 inline complex<double> polar(const double & magnitude, const double & angle){
839 return complex<double>(magnitude * ::cos(angle),magnitude * ::sin(angle));
840 }
841
842 // Transcendental functions implementation
843 template <typename ValueType>
844 __host__ __device__
845 inline complex<ValueType> cos(const complex<ValueType>& z){
846 const ValueType re = z.real();
847 const ValueType im = z.imag();
848 return complex<ValueType>(::cos(re) * ::cosh(im), -::sin(re) * ::sinh(im));
849 }
850
851 template <>
852 __host__ __device__
853 inline complex<float> cos(const complex<float>& z){
854 const float re = z.real();
855 const float im = z.imag();
856 return complex<float>(cosf(re) * coshf(im), -sinf(re) * sinhf(im));
857 }
858
859 template <typename ValueType>
860 __host__ __device__
861 inline complex<ValueType> cosh(const complex<ValueType>& z){
862 const ValueType re = z.real();
863 const ValueType im = z.imag();
864 return complex<ValueType>(::cosh(re) * ::cos(im), ::sinh(re) * ::sin(im));
865 }
866
867 template <>
868 __host__ __device__
869 inline complex<float> cosh(const complex<float>& z){
870 const float re = z.real();
871 const float im = z.imag();
872 return complex<float>(::coshf(re) * ::cosf(im), ::sinhf(re) * ::sinf(im));
873 }
874
875
876 template <typename ValueType>
877 __host__ __device__
878 inline complex<ValueType> exp(const complex<ValueType>& z){
879 return polar(::exp(z.real()),z.imag());
880 }
881
882 template <>
883 __host__ __device__
884 inline complex<float> exp(const complex<float>& z){
885 return polar(::expf(z.real()),z.imag());
886 }
887
888 template <typename ValueType>
889 __host__ __device__
890 inline complex<ValueType> log(const complex<ValueType>& z){
891 return complex<ValueType>(::log(abs(z)),arg(z));
892 }
893
894 template <>
895 __host__ __device__
896 inline complex<float> log(const complex<float>& z){
897 return complex<float>(::logf(abs(z)),arg(z));
898 }
899
900
901 template <typename ValueType>
902 __host__ __device__
903 inline complex<ValueType> log10(const complex<ValueType>& z){
904 // Using the explicit literal prevents compile time warnings in
905 // devices that don't support doubles
906 return log(z)/ValueType(2.30258509299404568402);
907 // return log(z)/ValueType(::log(10.0));
908 }
909
910 template <typename ValueType>
911 __host__ __device__
912 inline complex<ValueType> pow(const complex<ValueType>& z, const ValueType & exponent){
913 return exp(log(z)*exponent);
914 }
915
916 template <typename ValueType>
917 __host__ __device__
918 inline complex<ValueType> pow(const complex<ValueType>& z, const complex<ValueType> & exponent){
919 return exp(log(z)*exponent);
920 }
921
922 template <typename ValueType>
923 __host__ __device__
924 inline complex<ValueType> pow(const ValueType & x, const complex<ValueType> & exponent){
925 return exp(::log(x)*exponent);
926 }
927
928 template <>
929 __host__ __device__
930 inline complex<float> pow(const float & x, const complex<float> & exponent){
931 return exp(::logf(x)*exponent);
932 }
933
934 template <typename ValueType>
935 __host__ __device__
936 inline complex<ValueType> pow(const complex<ValueType>& z,const int & exponent){
937 return exp(log(z)*ValueType(exponent));
938 }
939
940 template <typename ValueType>
941 __host__ __device__
942 inline complex<ValueType> sin(const complex<ValueType>& z){
943 const ValueType re = z.real();
944 const ValueType im = z.imag();
945 return complex<ValueType>(::sin(re) * ::cosh(im), ::cos(re) * ::sinh(im));
946 }
947
948 template <>
949 __host__ __device__
950 inline complex<float> sin(const complex<float>& z){
951 const float re = z.real();
952 const float im = z.imag();
953 return complex<float>(::sinf(re) * ::coshf(im), ::cosf(re) * ::sinhf(im));
954 }
955
956 template <typename ValueType>
957 __host__ __device__
958 inline complex<ValueType> sinh(const complex<ValueType>& z){
959 const ValueType re = z.real();
960 const ValueType im = z.imag();
961 return complex<ValueType>(::sinh(re) * ::cos(im), ::cosh(re) * ::sin(im));
962 }
963
964 template <>
965 __host__ __device__
966 inline complex<float> sinh(const complex<float>& z){
967 const float re = z.real();
968 const float im = z.imag();
969 return complex<float>(::sinhf(re) * ::cosf(im), ::coshf(re) * ::sinf(im));
970 }
971
972 template <typename ValueType>
973 __host__ __device__
974 inline complex<ValueType> sqrt(const complex<ValueType>& z){
975 return polar(::sqrt(abs(z)),arg(z)/ValueType(2));
976 }
977
978 template <typename ValueType>
979 __host__ __device__
980 inline complex<float> sqrt(const complex<float>& z){
981 return polar(::sqrtf(abs(z)),arg(z)/float(2));
982 }
983
984 template <typename ValueType>
985 __host__ __device__
986 inline complex<ValueType> tan(const complex<ValueType>& z){
987 return sin(z)/cos(z);
988 }
989
990 template <typename ValueType>
991 __host__ __device__
992 inline complex<ValueType> tanh(const complex<ValueType>& z){
993 // This implementation seems better than the simple sin/cos
994 return (exp(ValueType(2)*z)-ValueType(1))/(exp(ValueType(2)*z)+ValueType(1));
995 // return sinh(z)/cosh(z);
996 }
997
998 // Inverse trigonometric functions implementation
999
1000 template <typename ValueType>
1001 __host__ __device__
1002 inline complex<ValueType> acos(const complex<ValueType>& z){
1003 const complex<ValueType> ret = asin(z);
1004 return complex<ValueType>(ValueType(M_PI/2.0) - ret.real(),-ret.imag());
1005 }
1006
1007 template <typename ValueType>
1008 __host__ __device__
1009 inline complex<ValueType> asin(const complex<ValueType>& z){
1010 const complex<ValueType> i(0,1);
1011 return -i*asinh(i*z);
1012 }
1013
1014 template <typename ValueType>
1015 __host__ __device__
1016 inline complex<ValueType> atan(const complex<ValueType>& z){
1017 const complex<ValueType> i(0,1);
1018 return -i*atanh(i*z);
1019 }
1020
1021 template <typename ValueType>
1022 __host__ __device__
1023 inline complex<ValueType> acosh(const complex<ValueType>& z){
1024 cusp::complex<ValueType> ret((z.real() - z.imag()) * (z.real() + z.imag()) - ValueType(1.0),
1025 ValueType(2.0) * z.real() * z.imag());
1026 ret = sqrt(ret);
1027 if (z.real() < ValueType(0.0)){
1028 ret = -ret;
1029 }
1030 ret += z;
1031 ret = log(ret);
1032 if (ret.real() < ValueType(0.0)){
1033 ret = -ret;
1034 }
1035 return ret;
1036
1037 /*
1038 cusp::complex<ValueType> ret = log(sqrt(z*z-ValueType(1))+z);
1039 if(ret.real() < 0){
1040 ret.real(-ret.real());
1041 }
1042 return ret;
1043 */
1044 }
1045
1046 template <typename ValueType>
1047 __host__ __device__
1048 inline complex<ValueType> asinh(const complex<ValueType>& z){
1049 return log(sqrt(z*z+ValueType(1))+z);
1050 }
1051
1052 template <typename ValueType>
1053 __host__ __device__
1054 inline complex<ValueType> atanh(const complex<ValueType>& z){
1055 ValueType imag2 = z.imag() * z.imag();
1056 ValueType n = ValueType(1.0) + z.real();
1057 n = imag2 + n * n;
1058
1059 ValueType d = ValueType(1.0) - z.real();
1060 d = imag2 + d * d;
1061 complex<ValueType> ret(ValueType(0.25) * (::log(n) - ::log(d)),0);
1062
1063 d = ValueType(1.0) - z.real() * z.real() - imag2;
1064
1065 ret.imag(ValueType(0.5) * ::atan2(ValueType(2.0) * z.imag(), d));
1066 return ret;
1067 //return (log(ValueType(1)+z)-log(ValueType(1)-z))/ValueType(2);
1068 }
1069
1070 template <typename ValueType>
1071 __host__ __device__
1072 inline complex<float> atanh(const complex<float>& z){
1073 float imag2 = z.imag() * z.imag();
1074 float n = float(1.0) + z.real();
1075 n = imag2 + n * n;
1076
1077 float d = float(1.0) - z.real();
1078 d = imag2 + d * d;
1079 complex<float> ret(float(0.25) * (::logf(n) - ::logf(d)),0);
1080
1081 d = float(1.0) - z.real() * z.real() - imag2;
1082
1083 ret.imag(float(0.5) * ::atan2f(float(2.0) * z.imag(), d));
1084 return ret;
1085 //return (log(ValueType(1)+z)-log(ValueType(1)-z))/ValueType(2);
1086
1087 }
1088
1089 } // end namespace cusp
1090
1091 #else
1092 #include <complex>
1093
1094 namespace cusp
1095 {
1096 using std::complex;
1097 using std::conj;
1098 using std::abs;
1099 using std::arg;
1100 using std::norm;
1101 using std::polar;
1102 using std::cos;
1103 using std::cosh;
1104 using std::exp;
1105 using std::log;
1106 using std::log10;
1107 using std::pow;
1108 using std::sin;
1109 using std::sinh;
1110 using std::sqrt;
1111 using std::tan;
1112 using std::tanh;
1113
1114 using std::acos;
1115 using std::asin;
1116 using std::atan;
1117 #if __cplusplus >= 201103L
1118 using std::acosh;
1119 using std::asinh;
1120 using std::atanh;
1121 #else
1122 template <typename ValueType>
1123 inline complex<ValueType> acosh(const complex<ValueType>& z){
1124 cusp::complex<ValueType> ret((z.real() - z.imag()) * (z.real() + z.imag()) - ValueType(1.0),
1125 ValueType(2.0) * z.real() * z.imag());
1126 ret = sqrt(ret);
1127 if (z.real() < ValueType(0.0)){
1128 ret = -ret;
1129 }
1130 ret += z;
1131 ret = log(ret);
1132 if (ret.real() < ValueType(0.0)){
1133 ret = -ret;
1134 }
1135 return ret;
1136 }
1137 template <typename ValueType>
1138 inline complex<ValueType> asinh(const complex<ValueType>& z){
1139 return log(sqrt(z*z+ValueType(1))+z);
1140 }
1141
1142 template <typename ValueType>
1143 inline complex<ValueType> atanh(const complex<ValueType>& z){
1144 ValueType imag2 = z.imag() * z.imag();
1145 ValueType n = ValueType(1.0) + z.real();
1146 n = imag2 + n * n;
1147
1148 ValueType d = ValueType(1.0) - z.real();
1149 d = imag2 + d * d;
1150 complex<ValueType> ret(ValueType(0.25) * (::log(n) - ::log(d)),0);
1151
1152 d = ValueType(1.0) - z.real() * z.real() - imag2;
1153
1154 ret.imag(ValueType(0.5) * ::atan2(ValueType(2.0) * z.imag(), d));
1155 return ret;
1156 }
1157 #endif
1158
1159
1160 template <typename T>
1161 struct norm_type {
1162 typedef T type;
1163 };
1164
1165 template <typename T>
1166 struct norm_type< complex<T> > {
1167 typedef T type;
1168 };
1169 }
1170 #endif

  ViewVC Help
Powered by ViewVC 1.1.26