/[escript]/trunk/cusplibrary/cusp/complex.h
ViewVC logotype

Annotation of /trunk/cusplibrary/cusp/complex.h

Parent Directory Parent Directory | Revision Log Revision Log


Revision 5588 - (hide annotations)
Tue Apr 21 05:07:16 2015 UTC (5 years, 2 months ago) by caltinay
File MIME type: text/plain
File size: 35591 byte(s)
Update to copyright headers in cusp library.

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

  ViewVC Help
Powered by ViewVC 1.1.26