complexfield.cuh
2.13 KB
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
#ifndef RTS_COMPLEXFIELD_H
#define RTS_COMPLEXFIELD_H
#include "cublas_v2.h"
#include <cuda_runtime.h>
#include "../math/field.cuh"
#include "../math/complex.h"
namespace rts{
/*This class stores functions for saving images of complex fields
*/
template<typename T, unsigned int D>
class complexfield : public field< rts::complex<T>, D >{
using field< rts::complex<T>, D >::R;
using field< rts::complex<T>, D >::X;
using field< rts::complex<T>, D >::shape;
public:
//find the maximum value of component n
rts::complex<T> find_max(unsigned int n){
cublasStatus_t stat;
cublasHandle_t handle;
//create a CUBLAS handle
stat = cublasCreate(&handle);
if(stat != CUBLAS_STATUS_SUCCESS){
std::cout<<"CUBLAS Error: initialization failed"<<std::endl;
exit(1);
}
int L = R[0] * R[1]; //compute the number of discrete points in a slice
int index; //result of the max operation
rts::complex<T> result;
if(sizeof(T) == 4)
stat = cublasIcamax(handle, L, (const cuComplex*)X[n], 1, &index);
else
stat = cublasIzamax(handle, L, (const cuDoubleComplex*)X[n], 1, &index);
index -= 1; //adjust for 1-based indexing
//if there was a GPU error, terminate
if(stat != CUBLAS_STATUS_SUCCESS){
std::cout<<"CUBLAS Error: failure finding maximum value."<<std::endl;
exit(1);
}
//retrieve the maximum value for this slice and store it in the maxVal array
std::cout<<X[n]<<std::endl;
HANDLE_ERROR(cudaMemcpy(&result, X[n] + index, sizeof(rts::complex<T>), cudaMemcpyDeviceToHost));
return result;
}
public:
//constructor (no parameters)
complexfield() : field<rts::complex<T>, D>(){};
//constructor (resolution specified)
complexfield(unsigned int r0, unsigned int r1) : field<rts::complex<T>, D>(r0, r1){};
//assignment operator (scalar value)
complexfield & operator= (const complex<T> rhs){
field< complex<T>, D >::operator=(rhs);
return *this;
}
//assignment operator (vector value)
complexfield & operator= (const vec< complex<T>, D > rhs){
field< complex<T>, D >::operator=(rhs);
return *this;
}
void toImage(std::string filename, unsigned int n){
}
};
} //end namespace rts
#endif