1/*
2 Copyright (C) 2011 X. Andrade
3
4 This program is free software; you can redistribute it and/or modify
5 it under the terms of the GNU General Public License as published by
6 the Free Software Foundation; either version 2, or (at your option)
7 any later version.
8
9 This program is distributed in the hope that it will be useful,
10 but WITHOUT ANY WARRANTY; without even the implied warranty of
11 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
12 GNU General Public License for more details.
13
14 You should have received a copy of the GNU General Public License
15 along with this program; if not, write to the Free Software
16 Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
17 02110-1301, USA.
18
19*/
20
21#include <cl_global.h>
22#include <cl_complex.h>
23
24#if defined(EXT_KHR_FP64) || defined(CUDA)
25#define HAVE_SINCOS
26#endif
27
28__kernel void phase(const int offset,
29		    const __global double * restrict phase,
30		    __global double2 * restrict psi, const int ldpsi){
31
32  const int ist = get_global_id(0);
33  const int ip  = get_global_id(1);
34
35#ifdef HAVE_SINCOS
36  double cc;
37  double ss = sincos(phase[offset + ip], &cc);
38#else
39#warning "Using single-precision sincos functions."
40  float fcc;
41  double ss = (double) sincos((float) phase[offset + ip], &fcc);
42  double cc = (double) fcc;
43#endif
44
45  const double2 ff = psi[(ip<<ldpsi) + ist];
46
47#ifdef CUDA
48  psi[(ip<<ldpsi) + ist] = double2(cc*ff.x + ss*ff.y, cc*ff.y - ss*ff.x);
49#else
50  psi[(ip<<ldpsi) + ist] = (double2)(cc*ff.x + ss*ff.y, cc*ff.y - ss*ff.x);
51#endif
52
53}
54
55__kernel void phase_hamiltonian(const int conjugate,
56				const int offset,
57				const int np,
58				const __global double2 * restrict phase,
59				__global const double2 * src, const int ldsrc,
60				__global double2 * psi, const int ldpsi){
61
62  const int ist = get_global_id(0);
63  const int ip  = get_global_id(1);
64
65  if(ip >= np) return;
66
67  if(conjugate){
68    psi[(ip<<ldpsi) + ist] = complex_mul(complex_conj(phase[offset + ip]), src[(ip<<ldpsi) + ist]);
69  } else {
70    psi[(ip<<ldpsi) + ist] = complex_mul(phase[offset + ip], src[(ip<<ldpsi) + ist]);
71  }
72
73}
74/*
75 Local Variables:
76 mode: c
77 coding: utf-8
78 End:
79*/
80