madrocketsci/amsculib3
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

working on comples64

Browse Source
This commit is contained in:
madrocketsci
2026-04-10 16:28:25 -04:00
parent a8ed51a904
commit f315dfb6f6
48 changed files with 551 additions and 2523 deletions
Download Patch File Download Diff File Expand all files Collapse all files

86
src/amsculib3/amscu_random.cu
View File

@ -163,59 +163,59 @@ __host__ __device__ float dpr64_randf(int64_t *seedinout)
void test_dprg64()
{
printf("Tests for dprg:\n");
long I;
int64_t seed = 133LL;
double d;
float f;
cuvect3 qv;
// printf("Tests for dprg:\n");
// long I;
// int64_t seed = 133LL;
// double d;
// float f;
// cuvec3 qv;
printf("dpr64_randd test\n");
seed = 133LL;
for(I=0;I<10;I++)
{
d = dpr64_randd(&seed);
printf("seed: %lld rand: %1.4f\n",(long long)seed,d);
}
// printf("dpr64_randd test\n");
// seed = 133LL;
// for(I=0;I<10;I++)
// {
// d = dpr64_randd(&seed);
// printf("seed: %lld rand: %1.4f\n",(long long)seed,d);
// }
printf("\n\n");
printf("dpr64_randf test\n");
seed = 133LL;
for(I=0;I<10;I++)
{
f = dpr64_randf(&seed);
printf("seed: %lld rand: %1.4f\n",(long long)seed,f);
}
// printf("\n\n");
// printf("dpr64_randf test\n");
// seed = 133LL;
// for(I=0;I<10;I++)
// {
// f = dpr64_randf(&seed);
// printf("seed: %lld rand: %1.4f\n",(long long)seed,f);
// }
return;
// return;
}
void test_dprg32()
{
printf("Tests for dprg:\n");
long I;
int32_t seed = 133;
double d;
float f;
cuvect3 qv;
// printf("Tests for dprg:\n");
// long I;
// int32_t seed = 133;
// double d;
// float f;
// cuvec3 qv;
printf("dpr32_randf test\n");
seed = 133;
for(I=0;I<10;I++)
{
f = dpr32_randf(&seed);
printf("seed: %lld rand: %1.4f\n",(long long)seed,f);
}
// printf("dpr32_randf test\n");
// seed = 133;
// for(I=0;I<10;I++)
// {
// f = dpr32_randf(&seed);
// printf("seed: %lld rand: %1.4f\n",(long long)seed,f);
// }
printf("\n\ndpr32_randnf test\n");
seed = 133;
for(I=0;I<10;I++)
{
f = dpr32_randnf(&seed);
printf("seed: %lld rand: %1.4f\n",(long long)seed,f);
}
// printf("\n\ndpr32_randnf test\n");
// seed = 133;
// for(I=0;I<10;I++)
// {
// f = dpr32_randnf(&seed);
// printf("seed: %lld rand: %1.4f\n",(long long)seed,f);
// }
return;
// return;
}

0
src/amsculib3/amscugeom.cu → src/amsculib3/geom/amscugeom.cu
View File

150
src/amsculib3/math/amscu_comp64.cu
View File

@ -7,15 +7,15 @@ namespace cmp
__host__ __device__ cucomp64::cucomp64()
{
real = 0.0;
imag = 0.0;
real = 0.0f;
imag = 0.0f;
return;
}
__host__ __device__ cucomp64::~cucomp64()
{
real = 0.0;
imag = 0.0;
real = 0.0f;
imag = 0.0f;
return;
}
@ -30,7 +30,14 @@ __host__ __device__ cucomp64::cucomp64(const cucomp64 &other)
__host__ __device__ cucomp64::cucomp64(const float &other)
{
real = other;
imag = 0.0;
imag = 0.0f;
return;
}
__host__ __device__ cucomp64::cucomp64(const float &_real, const float &_imag)
{
real = _real;
imag = _imag;
return;
}
@ -51,14 +58,14 @@ __host__ __device__ const cucomp64& cucomp64::operator=(const cucomp64& other)
__host__ __device__ cucomp64& cucomp64::operator=(float& other)
{
real = other;
imag = 0.0;
imag = 0.0f;
return *this;
}
__host__ __device__ const cucomp64& cucomp64::operator=(const float& other)
{
this->real = other;
this->imag = 0.0;
this->imag = 0.0f;
return *this;
}
@ -307,14 +314,21 @@ __host__ __device__ bool cucomp64::iszero() const
__host__ __device__ float cucomp64::arg() const
{
float ret = 0.0;
float ret = 0.0f;
ret = (float) amscuda::arg((double)real,(double)imag);
return ret;
}
__host__ __device__ float cucomp64::mag() const
{
float ret = 0.0;
float ret = 0.0f;
ret = ::sqrt(real*real+imag*imag);
return ret;
}
__host__ __device__ float cucomp64::abs() const
{
float ret = 0.0f;
ret = ::sqrt(real*real+imag*imag);
return ret;
}
@ -327,62 +341,54 @@ __host__ __device__ cucomp64 cucomp64::conj() const
return ret;
}
__host__ __device__ float arg(cucomp64 z)
__host__ __device__ float arg(const cucomp64 &z)
{
return z.arg();
}
__host__ __device__ float abs(cucomp64 z)
__host__ __device__ float abs(const cucomp64 &z)
{
return z.mag();
}
__host__ __device__ cucomp64 dtocomp64(float _r, float _i)
{
cucomp64 ret;
ret.real = _r;
ret.imag = _i;
return ret;
}
__host__ __device__ float real(cucomp64 z)
__host__ __device__ float real(const cucomp64 &z)
{
return z.real;
}
__host__ __device__ float imag(cucomp64 z)
__host__ __device__ float imag(const cucomp64 &z)
{
return z.imag;
}
__host__ __device__ cucomp64 sin(cucomp64 z)
__host__ __device__ cucomp64 sin(const cucomp64 &z)
{
cucomp64 ret;
cucomp64 im1 = dtocomp64(0.0f,1.0f);
cucomp64 div = dtocomp64(0.0f,2.0f);
cucomp64 im1 = cucomp64(0.0f,1.0f);
cucomp64 div = cucomp64(0.0f,2.0f);
ret = (exp(im1*z)-exp(-im1*z))/div;
return ret;
}
__host__ __device__ cucomp64 cos(cucomp64 z)
__host__ __device__ cucomp64 cos(const cucomp64 &z)
{
cucomp64 ret;
cucomp64 im1 = dtocomp64(0.0f,1.0f);
cucomp64 div = dtocomp64(2.0f,0.0f);
cucomp64 im1 = cucomp64(0.0f,1.0f);
cucomp64 div = cucomp64(2.0f,0.0f);
ret = (exp(im1*z)+exp(-im1*z))/div;
return ret;
}
__host__ __device__ cucomp64 tan(cucomp64 z)
__host__ __device__ cucomp64 tan(const cucomp64 &z)
{
return sin(z)/cos(z);
}
__host__ __device__ cucomp64 exp(cucomp64 z)
__host__ __device__ cucomp64 exp(const cucomp64 &z)
{
cucomp64 ret;
ret.real = ::exp(z.real)*::cos(z.imag);
@ -390,7 +396,7 @@ __host__ __device__ cucomp64 exp(cucomp64 z)
return ret;
}
__host__ __device__ cucomp64 log(cucomp64 z)
__host__ __device__ cucomp64 log(const cucomp64 &z)
{
cucomp64 ret;
ret.real = ::log(::sqrt(z.real*z.real+z.imag*z.imag));
@ -398,46 +404,103 @@ __host__ __device__ cucomp64 log(cucomp64 z)
return ret;
}
__host__ __device__ cucomp64 conj(cucomp64 z)
__host__ __device__ cucomp64 conj(const cucomp64 &z)
{
return z.conj();
}
__host__ __device__ cucomp64 cosh(cucomp64 z)
__host__ __device__ cucomp64 cosh(const cucomp64 &z)
{
cucomp64 ret;
cucomp64 div = dtocomp64(2.0f,0.0f);
cucomp64 div = cucomp64(2.0f,0.0f);
ret = (exp(z)+exp(-z))/div;
return ret;
}
__host__ __device__ cucomp64 sinh(cucomp64 z)
__host__ __device__ cucomp64 sinh(const cucomp64 &z)
{
cucomp64 ret;
cucomp64 div = dtocomp64(2.0f,0.0f);
cucomp64 div = cucomp64(2.0f,0.0f);
ret = (exp(z)-exp(-z))/div;
return ret;
}
__host__ __device__ cucomp64 tanh(cucomp64 z)
__host__ __device__ cucomp64 tanh(const cucomp64 &z)
{
return sinh(z)/cosh(z);
}
__host__ __device__ cucomp64 pow(cucomp64 z1, cucomp64 z2)
__host__ __device__ cucomp64 pow(const cucomp64 &z1, const cucomp64 &z2)
{
cucomp64 ret;
if(z1.mag()>0.0)
ret = exp(z2*log(z1));
if(z1.mag()>0.0f)
ret = amscuda::cmp::exp(amscuda::cmp::log(z1)*z2);
else
ret = dtocomp64(0.0f,0.0f);
ret = cucomp64(0.0f,0.0f);
return ret;
}
// //returns "complex sign" of complex number - 0, or a unit number with same argument
__host__ __device__ cucomp64 csgn(const cucomp64 &z)
{
cucomp64 ret;
float mag = z.mag();
if(mag>0.0f)
{
ret = z;
ret = ret/mag;
}
return ret;
}
//////////////////////////
//accumulation operators//
//////////////////////////
__host__ __device__ cucomp64& cucomp64::operator+=(const cucomp64& z)
{
return (*this);
}
__host__ __device__ cucomp64& cucomp64::operator-=(const cucomp64& z)
{
return (*this);
}
__host__ __device__ cucomp64& cucomp64::operator*=(const cucomp64& z)
{
return (*this);
}
__host__ __device__ cucomp64& cucomp64::operator/=(const cucomp64& z)
{
return (*this);
}
__host__ __device__ cucomp64& cucomp64::operator+=(const float& z)
{
return (*this);
}
__host__ __device__ cucomp64& cucomp64::operator-=(const float& z)
{
return (*this);
}
__host__ __device__ cucomp64& cucomp64::operator*=(const float& z)
{
return (*this);
}
__host__ __device__ cucomp64& cucomp64::operator/=(const float& z)
{
return (*this);
}
void test_cucomp64_1()
{
@ -448,11 +511,10 @@ void test_cucomp64_1()
printf("sizeof double=%ld\n",(long)(8*sizeof(f1)));
printf("sizeof double=%ld\n",(long)(8*sizeof(d1)));
printf("sizeof complex=%ld\n",(long)(8*sizeof(z1)));
printf("sizeof cucomp128=%ld\n",(long)(8*sizeof(a)));
printf("sizeof complex64=%ld\n",(long)(8*sizeof(z1)));
a = dtocomp64(1.0,1.0);
b = dtocomp64(1.0,-1.0);
a = cucomp64(1.0,1.0);
b = cucomp64(1.0,-1.0);
printf("a=%1.4f + %1.4fi\n",a[0],a[1]);
printf("b=%1.4f + %1.4fi\n",b[0],b[1]);

0
src/amsculib3/amscumath.cu → src/amsculib3/math/amscumath.cu
View File

148
src/amsculib3/math/cuvect2f.cu → src/amsculib3/math/cuvec2f.cu
View File

@ -3,88 +3,88 @@
namespace amscuda
{
__host__ __device__ cuvect2f::cuvect2f()
__host__ __device__ cuvec2f::cuvec2f()
{
x = 0.0f; y = 0.0f;
return;
}
__host__ __device__ cuvect2f::~cuvect2f()
__host__ __device__ cuvec2f::~cuvec2f()
{
x = 0.0f; y = 0.0f;
return;
}
__host__ __device__ float& cuvect2f::operator[](const int &I)
__host__ __device__ float& cuvec2f::operator[](const int &I)
{
if(I==0) return x;
if(I==1) return y;
return x;
}
__host__ __device__ const float& cuvect2f::operator[](const int &I) const
__host__ __device__ const float& cuvec2f::operator[](const int &I) const
{
if(I==0) return x;
if(I==1) return y;
return x;
}
__host__ __device__ cuvect2f cuvect2f::operator+(const cuvect2f &rhs)
__host__ __device__ cuvec2f cuvec2f::operator+(const cuvec2f &rhs)
{
cuvect2f ret;
cuvec2f ret;
ret.x = x+rhs.x;
ret.y = y+rhs.y;
return ret;
}
__host__ __device__ cuvect2f cuvect2f::operator-(const cuvect2f &rhs)
__host__ __device__ cuvec2f cuvec2f::operator-(const cuvec2f &rhs)
{
cuvect2f ret;
cuvec2f ret;
ret.x = x-rhs.x;
ret.y = y-rhs.y;
return ret;
}
__host__ __device__ cuvect2f cuvect2f::operator*(const float &rhs)
__host__ __device__ cuvec2f cuvec2f::operator*(const float &rhs)
{
cuvect2f ret;
cuvec2f ret;
ret.x = x*rhs;
ret.y = y*rhs;
return ret;
}
__host__ __device__ cuvect2f cuvect2f::operator/(const float &rhs)
__host__ __device__ cuvec2f cuvec2f::operator/(const float &rhs)
{
cuvect2f ret;
cuvec2f ret;
ret.x = x/rhs;
ret.y = y/rhs;
return ret;
}
__host__ __device__ cuvect2f& cuvect2f::operator+=(const cuvect2f &rhs)
__host__ __device__ cuvec2f& cuvec2f::operator+=(const cuvec2f &rhs)
{
x = x + rhs.x;
y = y + rhs.y;
return *this;
}
__host__ __device__ cuvect2f& cuvect2f::operator-=(const cuvect2f &rhs)
__host__ __device__ cuvec2f& cuvec2f::operator-=(const cuvec2f &rhs)
{
x = x - rhs.x;
y = y - rhs.y;
return *this;
}
__host__ __device__ cuvect2f& cuvect2f::operator*=(const float &rhs)
__host__ __device__ cuvec2f& cuvec2f::operator*=(const float &rhs)
{
x = x * rhs;
y = y * rhs;
return *this;
}
__host__ __device__ cuvect2f& cuvect2f::operator/=(const float &rhs)
__host__ __device__ cuvec2f& cuvec2f::operator/=(const float &rhs)
{
x = x / rhs;
y = y / rhs;
@ -92,34 +92,34 @@ namespace amscuda
}
__host__ __device__ cuvect2f::cuvect2f(const float &_x, const float &_y)
__host__ __device__ cuvec2f::cuvec2f(const float &_x, const float &_y)
{
x = _x; y = _y;
return;
}
__host__ __device__ float cuvect2f_dot(const cuvect2f &a, const cuvect2f &b)
__host__ __device__ float cuvec2f_dot(const cuvec2f &a, const cuvec2f &b)
{
float ret = a.x*b.x+a.y*b.y;
return ret;
}
__host__ __device__ float cuvect2f_cross(const cuvect2f &a, const cuvect2f &b)
__host__ __device__ float cuvec2f_cross(const cuvec2f &a, const cuvec2f &b)
{
float ret = a.x*b.y-a.y*b.x;
return ret;
}
__host__ __device__ float cuvect2f_norm(const cuvect2f &a)
__host__ __device__ float cuvec2f_norm(const cuvec2f &a)
{
float ret = ::sqrtf(a.x*a.x+a.y*a.y);
return ret;
}
__host__ __device__ cuvect2f cuvect2f_normalize(const cuvect2f &a)
__host__ __device__ cuvec2f cuvec2f_normalize(const cuvec2f &a)
{
cuvect2f ret;
float m = cuvect2f_norm(a);
cuvec2f ret;
float m = cuvec2f_norm(a);
if(m>0.0)
{
ret.x = a.x/m; ret.y = a.y/m;
@ -131,11 +131,11 @@ namespace amscuda
return ret;
}
__host__ __device__ cuvect2f cuvect2f_proj(const cuvect2f &a, const cuvect2f &b)
__host__ __device__ cuvec2f cuvec2f_proj(const cuvec2f &a, const cuvec2f &b)
{
cuvect2f ret;
cuvect2f bn = cuvect2f_normalize(b);
float m = cuvect2f_dot(a,bn);
cuvec2f ret;
cuvec2f bn = cuvec2f_normalize(b);
float m = cuvec2f_dot(a,bn);
ret = bn*m;
return ret;
}
@ -264,9 +264,9 @@ namespace amscuda
return ret;
}
__host__ __device__ cuvect2f cumat2f::operator*(const cuvect2f &rhs)
__host__ __device__ cuvec2f cumat2f::operator*(const cuvec2f &rhs)
{
cuvect2f ret;
cuvec2f ret;
ret.x = m00*rhs.x + m01*rhs.y;
ret.y = m10*rhs.x + m11*rhs.y;
@ -440,94 +440,16 @@ namespace amscuda
return R;
}
///////////////////////////
//legacy array operations//
///////////////////////////
//copies src to dest
__host__ __device__ void mat2f_copy(float *mat2f_dest, const float *mat2f_src)
{
mat2f_dest[0+0*2] = mat2f_src[0+0*2];
mat2f_dest[1+0*2] = mat2f_src[1+0*2];
mat2f_dest[0+1*2] = mat2f_src[0+1*2];
mat2f_dest[1+1*2] = mat2f_src[1+1*2];
return;
}
//inverts mat?inout[4]
__host__ __device__ void mat2f_inverse(float *mat2inout)
{
float det = mat2inout[0+0*2]*mat2inout[1+1*2]-mat2inout[0+1*2]*mat2inout[1+0*2];
float mat2in[4];
mat2f_copy(mat2in,mat2inout);
mat2inout[0+0*2] = mat2inout[1+1*2]/det;
mat2inout[1+0*2] = -mat2inout[1+0*2]/det;
mat2inout[0+1*2] = -mat2inout[0+1*2]/det;
mat2inout[1+1*2] = mat2inout[0+0*2]/det;
return;
}
//rotatin matrix from angle
__host__ __device__ void mat2f_rot_from_angle(float angle, float *mat2)
{
mat2[0+0*2] = ::cosf(angle);
mat2[1+0*2] = ::sinf(angle);
mat2[0+1*2] = -::sinf(angle);
mat2[1+1*2] = ::cosf(angle);
return;
}
//multiplies c = a*b
__host__ __device__ void mat2f_mult(float *mat2a, float *mat2b, float *mat2c)
{
float mat2a_in[4];
float mat2b_in[4];
if(mat2a==NULL || mat2b==NULL || mat2c==NULL)
__host__ __device__ cuvec2f operator-(const cuvec2f &rhs)
{
return;
cuvec2f ret;
ret[0] = -rhs[0];
ret[1] = -rhs[1];
return ret;
}
mat2f_copy(mat2a_in,mat2a);
mat2f_copy(mat2b_in,mat2b);
mat2c[0+0*2] = mat2a_in[0+0*2]*mat2b_in[0+0*2] + mat2a_in[1+0*2]*mat2b_in[0+1*2];
mat2c[1+0*2] = mat2a_in[0+0*2]*mat2b_in[1+0*2] + mat2a_in[1+0*2]*mat2b_in[1+1*2];
mat2c[0+1*2] = mat2a_in[0+1*2]*mat2b_in[0+0*2] + mat2a_in[1+1*2]*mat2b_in[0+1*2];
mat2c[1+1*2] = mat2a_in[0+1*2]*mat2b_in[1+0*2] + mat2a_in[1+1*2]*mat2b_in[1+1*2];
return;
}
// ret = a*b
__host__ __device__ cuvect2f mat2f_mult(float *mat2a, const cuvect2f &b)
{
cuvect2f ret;
ret.x = b.x*mat2a[0+0*2] + b.y*mat2a[1+0*2];
ret.y = b.x*mat2a[0+1*2] + b.y*mat2a[1+1*2];
return ret;
}
__host__ __device__ cuvect2f operator-(cuvect2f rhs)
{
cuvect2f ret;
ret[0] = -rhs[0];
ret[1] = -rhs[1];
return ret;
}
__host__ __device__ cumat2f operator-(cumat2f rhs)
{
cumat2f ret;
int I;
for(I=0;I<4;I++) ret[I] = -rhs[I];
return ret;
}
void test_cuvect2f_1()
void test_cuvec2f_1()
{

94
src/amsculib3/math/cuvect3f.cu → src/amsculib3/math/cuvec3f.cu
View File

@ -3,19 +3,19 @@
namespace amscuda
{
__host__ __device__ cuvect3f::cuvect3f()
__host__ __device__ cuvec3f::cuvec3f()
{
x = 0.0f; y = 0.0f; z = 0.0f;
return;
}
__host__ __device__ cuvect3f::~cuvect3f()
__host__ __device__ cuvec3f::~cuvec3f()
{
x = 0.0f; y = 0.0f; z = 0.0f;
return;
}
__host__ __device__ float& cuvect3f::operator[](const int &I)
__host__ __device__ float& cuvec3f::operator[](const int &I)
{
if(I==0) return x;
if(I==1) return y;
@ -23,7 +23,7 @@ namespace amscuda
return x;
}
__host__ __device__ const float& cuvect3f::operator[](const int &I) const
__host__ __device__ const float& cuvec3f::operator[](const int &I) const
{
if(I==0) return x;
if(I==1) return y;
@ -31,9 +31,9 @@ namespace amscuda
return x;
}
__host__ __device__ cuvect3f cuvect3f::operator+(const cuvect3f &rhs)
__host__ __device__ cuvec3f cuvec3f::operator+(const cuvec3f &rhs)
{
cuvect3f ret;
cuvec3f ret;
ret.x = x+rhs.x;
ret.y = y+rhs.y;
ret.z = z+rhs.z;
@ -41,9 +41,9 @@ namespace amscuda
return ret;
}
__host__ __device__ cuvect3f cuvect3f::operator-(const cuvect3f &rhs)
__host__ __device__ cuvec3f cuvec3f::operator-(const cuvec3f &rhs)
{
cuvect3f ret;
cuvec3f ret;
ret.x = x-rhs.x;
ret.y = y-rhs.y;
ret.z = z-rhs.z;
@ -51,25 +51,25 @@ namespace amscuda
return ret;
}
__host__ __device__ cuvect3f cuvect3f::operator*(const float &rhs)
__host__ __device__ cuvec3f cuvec3f::operator*(const float &rhs)
{
cuvect3f ret;
cuvec3f ret;
ret.x = x*rhs;
ret.y = y*rhs;
ret.z = z*rhs;
return ret;
}
__host__ __device__ cuvect3f cuvect3f::operator/(const float &rhs)
__host__ __device__ cuvec3f cuvec3f::operator/(const float &rhs)
{
cuvect3f ret;
cuvec3f ret;
ret.x = x/rhs;
ret.y = y/rhs;
ret.z = z/rhs;
return ret;
}
__host__ __device__ cuvect3f& cuvect3f::operator+=(const cuvect3f &rhs)
__host__ __device__ cuvec3f& cuvec3f::operator+=(const cuvec3f &rhs)
{
x = x + rhs.x;
y = y + rhs.y;
@ -77,7 +77,7 @@ namespace amscuda
return *this;
}
__host__ __device__ cuvect3f& cuvect3f::operator-=(const cuvect3f &rhs)
__host__ __device__ cuvec3f& cuvec3f::operator-=(const cuvec3f &rhs)
{
x = x - rhs.x;
y = y - rhs.y;
@ -85,7 +85,7 @@ namespace amscuda
return *this;
}
__host__ __device__ cuvect3f& cuvect3f::operator*=(const float &rhs)
__host__ __device__ cuvec3f& cuvec3f::operator*=(const float &rhs)
{
x = x * rhs;
y = y * rhs;
@ -93,7 +93,7 @@ namespace amscuda
return *this;
}
__host__ __device__ cuvect3f& cuvect3f::operator/=(const float &rhs)
__host__ __device__ cuvec3f& cuvec3f::operator/=(const float &rhs)
{
x = x / rhs;
y = y / rhs;
@ -102,23 +102,23 @@ namespace amscuda
}
__host__ __device__ cuvect3f::cuvect3f(const float &_x, const float &_y, const float &_z)
__host__ __device__ cuvec3f::cuvec3f(const float &_x, const float &_y, const float &_z)
{
x = _x; y = _y; z = _z;
return;
}
__host__ __device__ float cuvect3f_dot(const cuvect3f &a, const cuvect3f &b)
__host__ __device__ float cuvec3f_dot(const cuvec3f &a, const cuvec3f &b)
{
float ret = a.x*b.x+a.y*b.y+a.z*b.z;
return ret;
}
__host__ __device__ cuvect3f cuvect3f_cross(const cuvect3f &a, const cuvect3f &b)
__host__ __device__ cuvec3f cuvec3f_cross(const cuvec3f &a, const cuvec3f &b)
{
cuvect3f ret;
cuvec3f ret;
ret[0] = a[1]*b[2]-a[2]*b[1];
ret[1] = a[2]*b[0]-a[0]*b[2];
ret[2] = a[0]*b[1]-a[1]*b[0];
@ -126,16 +126,16 @@ namespace amscuda
return ret;
}
__host__ __device__ float cuvect3f_norm(const cuvect3f &a)
__host__ __device__ float cuvec3f_norm(const cuvec3f &a)
{
float ret;
ret = ::sqrtf(a.x*a.x+a.y*a.y+a.z*a.z);
return ret;
}
__host__ __device__ cuvect3f cuvect3f_normalize(const cuvect3f &a)
__host__ __device__ cuvec3f cuvec3f_normalize(const cuvec3f &a)
{
cuvect3f ret;
cuvec3f ret;
float m;
m = ::sqrtf(a.x*a.x+a.y*a.y+a.z*a.z);
if(m>0.0)
@ -150,11 +150,11 @@ namespace amscuda
return ret;
}
__host__ __device__ cuvect3f cuvect3f_proj(const cuvect3f &a, const cuvect3f &b)
__host__ __device__ cuvec3f cuvec3f_proj(const cuvec3f &a, const cuvec3f &b)
{
cuvect3f ret;
cuvect3f bn = cuvect3f_normalize(b);
float m = cuvect3f_dot(a,bn);
cuvec3f ret;
cuvec3f bn = cuvec3f_normalize(b);
float m = cuvec3f_dot(a,bn);
ret = bn*m;
return ret;
}
@ -361,9 +361,9 @@ __host__ __device__ cumat3f cumat3f::operator/(const float &rhs)
return ret;
}
__host__ __device__ cuvect3f cumat3f::operator*(const cuvect3f &rhs)
__host__ __device__ cuvec3f cumat3f::operator*(const cuvec3f &rhs)
{
cuvect3f ret;
cuvec3f ret;
ret.x = m00*rhs.x + m01*rhs.y + m02*rhs.z;
ret.y = m10*rhs.x + m11*rhs.y + m12*rhs.z;
@ -597,7 +597,7 @@ __host__ void test_cudavectf_logic1()
// printf("3 dim vector and matrix functional tests on host side\n");
// cuvect3f a,b,c;
// cuvec3f a,b,c;
// float ma[9],mb[9],mc[9];
// int I,J;
@ -636,7 +636,7 @@ __host__ void test_cudavectf_logic1()
// }
// }
// a = cuvect3f(1,1,1);
// a = cuvec3f(1,1,1);
// b = mat3f_mult(ma,a);
// b = mat3f_mult(mb,b);
@ -645,8 +645,8 @@ __host__ void test_cudavectf_logic1()
// printf("a[%d] = %1.3f, b[%d] = %1.3f\n",I,a[I],I,b[I]);
// }
// a = cuvect3f(1,0,1);
// b = cuvect3f(0,1,-1);
// a = cuvec3f(1,0,1);
// b = cuvec3f(0,1,-1);
// c = a+b;
// for(I=0;I<3;I++)
@ -661,31 +661,31 @@ __host__ void test_cudavectf_logic1()
// printf("a[%d] = %1.3f, b[%d] = %1.3f, c[%d] = %1.3f\n",I,a[I],I,b[I],I,c[I]);
// }
// c = cuvect3f_cross(a,b);
// c = cuvec3f_cross(a,b);
// for(I=0;I<3;I++)
// {
// printf("a[%d] = %1.3f, b[%d] = %1.3f, c[%d] = %1.3f\n",I,a[I],I,b[I],I,c[I]);
// }
// printf("c dot a = %1.3f, c dot b = %1.3f\n",cuvect3f_dot(c,a),cuvect3f_dot(c,b));
// printf("c dot a = %1.3f, c dot b = %1.3f\n",cuvec3f_dot(c,a),cuvec3f_dot(c,b));
// printf("norm(a)=%1.3f, norm(b)=%1.3f, norm(c)=%1.3f\n",cuvect3f_norm(a),cuvect3f_norm(b),cuvect3f_norm(c));
// a = cuvect3f_normalize(a);
// b = cuvect3f_normalize(b);
// c = cuvect3f_normalize(c);
// printf("norm(a)=%1.3f, norm(b)=%1.3f, norm(c)=%1.3f\n",cuvec3f_norm(a),cuvec3f_norm(b),cuvec3f_norm(c));
// a = cuvec3f_normalize(a);
// b = cuvec3f_normalize(b);
// c = cuvec3f_normalize(c);
// for(I=0;I<3;I++)
// {
// printf("a[%d] = %1.3f, b[%d] = %1.3f, c[%d] = %1.3f\n",I,a[I],I,b[I],I,c[I]);
// }
// printf("c dot a = %1.3f, c dot b = %1.3f\n",cuvect3f_dot(c,a),cuvect3f_dot(c,b));
// printf("norm(a)=%1.3f, norm(b)=%1.3f, norm(c)=%1.3f\n",cuvect3f_norm(a),cuvect3f_norm(b),cuvect3f_norm(c));
// printf("c dot a = %1.3f, c dot b = %1.3f\n",cuvec3f_dot(c,a),cuvec3f_dot(c,b));
// printf("norm(a)=%1.3f, norm(b)=%1.3f, norm(c)=%1.3f\n",cuvec3f_norm(a),cuvec3f_norm(b),cuvec3f_norm(c));
// return;
}
__host__ __device__ cumat3f hodge_dual(const cuvect3f &vin)
__host__ __device__ cumat3f hodge_dual(const cuvec3f &vin)
{
cumat3f ret;
@ -704,9 +704,9 @@ __host__ __device__ cumat3f hodge_dual(const cuvect3f &vin)
return ret;
}
__host__ __device__ cuvect3f hodge_dual(const cumat3f &min)
__host__ __device__ cuvec3f hodge_dual(const cumat3f &min)
{
cuvect3f ret;
cuvec3f ret;
ret.x = 0.5f*(min.m12 - min.m21);
ret.y = 0.5f*(min.m20 - min.m02);
@ -731,13 +731,13 @@ __host__ __device__ const cumat3f cumat3f_zeros()
return ret;
}
__host__ __device__ cumat3f rotmat_from_axisangle(const cuvect3f &axis, const float &angle)
__host__ __device__ cumat3f rotmat_from_axisangle(const cuvec3f &axis, const float &angle)
{
cumat3f ret = cumat3f_zeros();
cumat3f H;
cuvect3f _naxis;
cuvec3f _naxis;
_naxis = cuvect3f_normalize(axis);
_naxis = cuvec3f_normalize(axis);
H = hodge_dual(_naxis);
ret += H*sinf(angle);

64
src/amsculib3/math/cuvect4f.cu → src/amsculib3/math/cuvec4f.cu
View File

@ -4,22 +4,22 @@ namespace amscuda
{
////////////
//cuvect4ff//
//cuvec4ff//
////////////
__host__ __device__ cuvect4f::cuvect4f()
__host__ __device__ cuvec4f::cuvec4f()
{
x = 0.0f; y = 0.0f; z = 0.0f; w = 0.0f;
return;
}
__host__ __device__ cuvect4f::~cuvect4f()
__host__ __device__ cuvec4f::~cuvec4f()
{
x = 0.0f; y = 0.0f; z = 0.0f; w = 0.0f;
return;
}
__host__ __device__ float& cuvect4f::operator[](const int &I)
__host__ __device__ float& cuvec4f::operator[](const int &I)
{
if(I==0) return x;
if(I==1) return y;
@ -28,7 +28,7 @@ __host__ __device__ float& cuvect4f::operator[](const int &I)
return x;
}
__host__ __device__ const float& cuvect4f::operator[](const int &I) const
__host__ __device__ const float& cuvec4f::operator[](const int &I) const
{
if(I==0) return x;
if(I==1) return y;
@ -37,9 +37,9 @@ __host__ __device__ const float& cuvect4f::operator[](const int &I) const
return x;
}
__host__ __device__ cuvect4f cuvect4f::operator+(const cuvect4f &rhs)
__host__ __device__ cuvec4f cuvec4f::operator+(const cuvec4f &rhs)
{
cuvect4f ret;
cuvec4f ret;
ret.x = x+rhs.x;
ret.y = y+rhs.y;
ret.z = z+rhs.z;
@ -48,9 +48,9 @@ __host__ __device__ cuvect4f cuvect4f::operator+(const cuvect4f &rhs)
return ret;
}
__host__ __device__ cuvect4f cuvect4f::operator-(const cuvect4f &rhs)
__host__ __device__ cuvec4f cuvec4f::operator-(const cuvec4f &rhs)
{
cuvect4f ret;
cuvec4f ret;
ret.x = x-rhs.x;
ret.y = y-rhs.y;
ret.z = z-rhs.z;
@ -59,9 +59,9 @@ __host__ __device__ cuvect4f cuvect4f::operator-(const cuvect4f &rhs)
return ret;
}
__host__ __device__ cuvect4f cuvect4f::operator*(const float &rhs)
__host__ __device__ cuvec4f cuvec4f::operator*(const float &rhs)
{
cuvect4f ret;
cuvec4f ret;
ret.x = x*rhs;
ret.y = y*rhs;
ret.z = z*rhs;
@ -69,9 +69,9 @@ __host__ __device__ cuvect4f cuvect4f::operator*(const float &rhs)
return ret;
}
__host__ __device__ cuvect4f cuvect4f::operator/(const float &rhs)
__host__ __device__ cuvec4f cuvec4f::operator/(const float &rhs)
{
cuvect4f ret;
cuvec4f ret;
ret.x = x/rhs;
ret.y = y/rhs;
ret.z = z/rhs;
@ -79,7 +79,7 @@ __host__ __device__ cuvect4f cuvect4f::operator/(const float &rhs)
return ret;
}
__host__ __device__ cuvect4f& cuvect4f::operator+=(const cuvect4f &rhs)
__host__ __device__ cuvec4f& cuvec4f::operator+=(const cuvec4f &rhs)
{
x = x + rhs.x;
y = y + rhs.y;
@ -88,7 +88,7 @@ __host__ __device__ cuvect4f& cuvect4f::operator+=(const cuvect4f &rhs)
return *this;
}
__host__ __device__ cuvect4f& cuvect4f::operator-=(const cuvect4f &rhs)
__host__ __device__ cuvec4f& cuvec4f::operator-=(const cuvec4f &rhs)
{
x = x - rhs.x;
y = y - rhs.y;
@ -97,7 +97,7 @@ __host__ __device__ cuvect4f& cuvect4f::operator-=(const cuvect4f &rhs)
return *this;
}
__host__ __device__ cuvect4f& cuvect4f::operator*=(const float &rhs)
__host__ __device__ cuvec4f& cuvec4f::operator*=(const float &rhs)
{
x = x * rhs;
y = y * rhs;
@ -106,7 +106,7 @@ __host__ __device__ cuvect4f& cuvect4f::operator*=(const float &rhs)
return *this;
}
__host__ __device__ cuvect4f& cuvect4f::operator/=(const float &rhs)
__host__ __device__ cuvec4f& cuvec4f::operator/=(const float &rhs)
{
x = x / rhs;
y = y / rhs;
@ -115,9 +115,9 @@ __host__ __device__ cuvect4f& cuvect4f::operator/=(const float &rhs)
return *this;
}
__host__ __device__ cuvect4f operator-(const cuvect4f &rhs)
__host__ __device__ cuvec4f operator-(const cuvec4f &rhs)
{
cuvect4f ret;
cuvec4f ret;
ret[0] = -rhs[0];
ret[1] = -rhs[1];
ret[2] = -rhs[2];
@ -126,7 +126,7 @@ __host__ __device__ cuvect4f operator-(const cuvect4f &rhs)
}
__host__ __device__ cuvect4f::cuvect4f(const float &_x, const float &_y, const float &_z, const float &_w)
__host__ __device__ cuvec4f::cuvec4f(const float &_x, const float &_y, const float &_z, const float &_w)
{
x = _x; y = _y; z = _z; w = _w;
return;
@ -404,9 +404,9 @@ __host__ __device__ cumat4f cumat4f::operator/(const float &rhs)
return ret;
}
__host__ __device__ cuvect4f cumat4f::operator*(const cuvect4f &rhs)
__host__ __device__ cuvec4f cumat4f::operator*(const cuvec4f &rhs)
{
cuvect4f ret;
cuvec4f ret;
ret.x = m00*rhs.x + m01*rhs.y + m02*rhs.z + m03*rhs.w;
ret.y = m10*rhs.x + m11*rhs.y + m12*rhs.z + m13*rhs.w;
ret.z = m20*rhs.x + m21*rhs.y + m22*rhs.z + m23*rhs.w;
@ -674,34 +674,34 @@ __host__ __device__ cumat4f& cumat4f::operator*=(const cumat4f &rhs)
}
__host__ __device__ float cuvect4f_dot(cuvect4f &a, cuvect4f &b)
__host__ __device__ float cuvec4f_dot(cuvec4f &a, cuvec4f &b)
{
float ret = 0.0f;
ret = a.x*b.x + a.y*b.y + a.z*b.z + a.w*b.w;
return ret;
}
__host__ __device__ float cuvect4f_norm(cuvect4f &a)
__host__ __device__ float cuvec4f_norm(cuvec4f &a)
{
float ret = 0.0f;
ret = ::sqrtf(cuvect4f_dot(a,a));
ret = ::sqrtf(cuvec4f_dot(a,a));
return ret;
}
__host__ __device__ cuvect4f cuvect4f_normalize(cuvect4f &a)
__host__ __device__ cuvec4f cuvec4f_normalize(cuvec4f &a)
{
cuvect4f ret = cuvect4f(0.0f,0.0f,0.0f,0.0f);
float nrm = cuvect4f_norm(a);
cuvec4f ret = cuvec4f(0.0f,0.0f,0.0f,0.0f);
float nrm = cuvec4f_norm(a);
if(nrm>0.0f)
ret = a/nrm;
return ret;
}
__host__ __device__ cuvect4f cuvect4f_proj(cuvect4f &a, cuvect4f &b)
__host__ __device__ cuvec4f cuvec4f_proj(cuvec4f &a, cuvec4f &b)
{
cuvect4f ret;
cuvect4f bn = cuvect4f_normalize(b);
float d = cuvect4f_dot(a,bn);
cuvec4f ret;
cuvec4f bn = cuvec4f_normalize(b);
float d = cuvec4f_dot(a,bn);
ret = bn*d;
return ret;
}

377
src/amsculib3/math/cuvect2.cu
View File

@ -1,377 +0,0 @@
#include <amsculib3/amsculib3.hpp>
namespace amscuda
{
__host__ __device__ cuvect2::cuvect2()
{
x = 0.0; y = 0.0;
return;
}
__host__ __device__ cuvect2::~cuvect2()
{
x = 0.0; y = 0.0;
return;
}
__host__ __device__ cuvect2::cuvect2(double _x, double _y)
{
x = _x; y = _y;
return;
}
__host__ __device__ double& cuvect2::operator[](const int I)
{
if(I==0) return x;
if(I==1) return y;
return x;
}
__host__ __device__ const double& cuvect2::operator[](const int I) const
{
if(I==0) return x;
if(I==1) return y;
return x;
}
__host__ __device__ cuvect2 cuvect2::operator+(cuvect2 lhs)
{
cuvect2 ret;
ret.x = x + lhs.x;
ret.y = y + lhs.y;
return ret;
}
__host__ __device__ cuvect2 cuvect2::operator-(cuvect2 lhs)
{
cuvect2 ret;
ret.x = x - lhs.x;
ret.y = y - lhs.y;
return ret;
}
__host__ __device__ cuvect2 cuvect2::operator*(double lhs)
{
cuvect2 ret;
ret.x = x*lhs;
ret.y = y*lhs;
return ret;
}
__host__ __device__ cuvect2 cuvect2::operator/(double lhs)
{
cuvect2 ret;
ret.x = x/lhs;
ret.y = y/lhs;
return ret;
}
__host__ __device__ double cuvect2_dot(cuvect2 a, cuvect2 b)
{
double ret = a.x*b.x+a.y*b.y;
return ret;
}
__host__ __device__ double cuvect2_cross(cuvect2 a, cuvect2 b)
{
double ret = a.x*b.y-a.y*b.x;
return ret;
}
__host__ __device__ double cuvect2_norm(cuvect2 a)
{
double ret = ::sqrt(a.x*a.x+a.y*a.y);
return ret;
}
__host__ __device__ cuvect2 cuvect2_normalize(cuvect2 a)
{
cuvect2 ret;
double m = cuvect2_norm(a);
if(m>0.0)
{
ret.x = a.x/m; ret.y = a.y/m;
}
else
{
ret.x = 0.0; ret.y = 0.0;
}
return ret;
}
__host__ __device__ cuvect2 cuvect2_proj(cuvect2 a, cuvect2 b)
{
cuvect2 ret;
cuvect2 bn = cuvect2_normalize(b);
double m = cuvect2_dot(a,bn);
ret = bn*m;
return ret;
}
__host__ __device__ cumat2::cumat2()
{
int I;
for(I=0;I<4;I++)
{
dat[I] = 0.0;
}
return;
}
__host__ __device__ cumat2::~cumat2()
{
int I;
for(I=0;I<4;I++)
{
dat[I] = 0.0;
}
return;
}
__host__ __device__ double& cumat2::operator[](const int I)
{
return dat[I];
}
__host__ __device__ double& cumat2::operator()(const int I, const int J)
{
return dat[I+2*J];
}
__host__ __device__ cumat2 cumat2::operator+(cumat2 lhs)
{
int I;
cumat2 ret;
for(I=0;I<4;I++)
{
ret.dat[I] = this->dat[I] + lhs.dat[I];
}
return ret;
}
__host__ __device__ cumat2 cumat2::operator-(cumat2 lhs)
{
int I;
cumat2 ret;
for(I=0;I<4;I++)
{
ret.dat[I] = this->dat[I] - lhs.dat[I];
}
return ret;
}
__host__ __device__ cumat2 cumat2::operator*(double lhs)
{
cumat2 ret;
int I;
for(I=0;I<4;I++)
{
ret.dat[I] = this->dat[I]*lhs;
}
return ret;
}
__host__ __device__ cumat2 cumat2::operator/(double lhs)
{
cumat2 ret;
int I;
for(I=0;I<4;I++)
{
ret.dat[I] = this->dat[I]/lhs;
}
return ret;
}
__host__ __device__ double& cumat2::at(const int I, const int J)
{
return dat[I+2*J];
}
__host__ __device__ cuvect2 cumat2::operator*(cuvect2 lhs)
{
cuvect2 ret = cuvect2(0.0,0.0);
int I,J;
for(I=0;I<2;I++)
{
for(J=0;J<2;J++)
{
ret[I] = ret[I] + this->at(I,J)*lhs[J];
}
}
return ret;
}
__host__ __device__ cumat2 cumat2::operator*(cumat2 lhs)
{
cumat2 ret;
int I,J,K;
for(I=0;I<2;I++)
{
for(J=0;J<2;J++)
{
ret(I,J) = 0.0;
for(K=0;K<2;K++)
{
ret(I,J) = ret(I,J) + this->at(I,K)*lhs(K,J);
}
}
}
return ret;
}
__host__ __device__ double cumat2::det()
{
double ret = 0.0;
ret = ret + this->at(0,0)*this->at(1,1);
ret = ret - this->at(1,0)*this->at(0,1);
return ret;
}
__host__ __device__ cumat2 cumat2::transpose()
{
cumat2 q;
int I,J;
for(I=0;I<2;I++)
{
for(J=0;J<2;J++)
{
q.at(I,J) = this->at(J,I);
}
}
return q;
}
__host__ __device__ cumat2 cumat2::inverse()
{
cumat2 q;
double dt = q.det();
if(dt!=0.0)
{
q(0,0) = this->at(1,1)/dt;
q(0,1) = -this->at(0,1)/dt;
q(1,0) = -this->at(1,0)/dt;
q(1,1) = this->at(0,0)/dt;
}
else
{
q(0,0) = inf;
q(0,1) = inf;
q(1,0) = inf;
q(1,1) = inf;
}
return q;
}
//2x2 matrix operations
//matrix order is assumed to be mat[I,J] = mat[I+3*J]
//transpose a 2x2 matrix in place
__host__ __device__ void mat2_transpose(double *mat2inout)
{
double mat2_in[4];
mat2_copy(mat2_in,mat2inout);
mat2inout[0+0*2] = mat2_in[0+0*2];
mat2inout[1+0*2] = mat2_in[0+1*2];
mat2inout[0+1*2] = mat2_in[1+0*2];
mat2inout[1+1*2] = mat2_in[1+1*2];
return;
}
//copies src to dest
__host__ __device__ void mat2_copy(double *mat2_dest, const double *mat2_src)
{
mat2_dest[0+0*2] = mat2_src[0+0*2];
mat2_dest[1+0*2] = mat2_src[1+0*2];
mat2_dest[0+1*2] = mat2_src[0+1*2];
mat2_dest[1+1*2] = mat2_src[1+1*2];
return;
}
//inverts mat?inout[4]
__host__ __device__ void mat2_inverse(double *mat2inout)
{
double det = mat2inout[0+0*2]*mat2inout[1+1*2]-mat2inout[0+1*2]*mat2inout[1+0*2];
double mat2in[4];
mat2_copy(mat2in,mat2inout);
mat2inout[0+0*2] = mat2inout[1+1*2]/det;
mat2inout[1+0*2] = -mat2inout[1+0*2]/det;
mat2inout[0+1*2] = -mat2inout[0+1*2]/det;
mat2inout[1+1*2] = mat2inout[0+0*2]/det;
return;
}
//rotatin matrix from angle
__host__ __device__ void mat2_rot_from_angle(double angle, double *mat2)
{
mat2[0+0*2] = ::cos(angle);
mat2[1+0*2] = ::sin(angle);
mat2[0+1*2] = -::sin(angle);
mat2[1+1*2] = ::cos(angle);
return;
}
//multiplies c = a*b
__host__ __device__ void mat2_mult(double *mat2a, double *mat2b, double *mat2c)
{
double mat2a_in[4];
double mat2b_in[4];
if(mat2a==NULL || mat2b==NULL || mat2c==NULL)
{
return;
}
mat2_copy(mat2a_in,mat2a);
mat2_copy(mat2b_in,mat2b);
mat2c[0+0*2] = mat2a_in[0+0*2]*mat2b_in[0+0*2] + mat2a_in[1+0*2]*mat2b_in[0+1*2];
mat2c[1+0*2] = mat2a_in[0+0*2]*mat2b_in[1+0*2] + mat2a_in[1+0*2]*mat2b_in[1+1*2];
mat2c[0+1*2] = mat2a_in[0+1*2]*mat2b_in[0+0*2] + mat2a_in[1+1*2]*mat2b_in[0+1*2];
mat2c[1+1*2] = mat2a_in[0+1*2]*mat2b_in[1+0*2] + mat2a_in[1+1*2]*mat2b_in[1+1*2];
return;
}
// ret = a*b
__host__ __device__ cuvect2 mat2_mult(double *mat2a, cuvect2 b)
{
cuvect2 ret;
ret.x = b.x*mat2a[0+0*2] + b.y*mat2a[1+0*2];
ret.y = b.x*mat2a[0+1*2] + b.y*mat2a[1+1*2];
return ret;
}
__host__ __device__ cuvect2 operator-(cuvect2 rhs)
{
cuvect2 ret;
ret[0] = -rhs[0];
ret[1] = -rhs[1];
return ret;
}
__host__ __device__ cumat2 operator-(cumat2 rhs)
{
cumat2 ret;
int I;
for(I=0;I<4;I++) ret[I] = -rhs[I];
return ret;
}
void test_cuvect2_1()
{
return;
}
};

860
src/amsculib3/math/cuvect3.cu
View File

@ -1,860 +0,0 @@
#include <amsculib3/amsculib3.hpp>
namespace amscuda
{
__host__ __device__ cuvect3::cuvect3()
{
x = 0.0; y = 0.0; z = 0.0;
return;
}
__host__ __device__ cuvect3::~cuvect3()
{
x = 0.0; y = 0.0; z = 0.0;
return;
}
__host__ __device__ double& cuvect3::operator[](const int &I)
{
if(I==0) return x;
if(I==1) return y;
if(I==2) return z;
return x;
}
__host__ __device__ const double& cuvect3::operator[](const int &I) const
{
if(I==0) return x;
if(I==1) return y;
if(I==2) return z;
return x;
}
__host__ __device__ cuvect3 cuvect3::operator+(const cuvect3 &rhs)
{
cuvect3 ret;
ret.x = x+rhs.x;
ret.y = y+rhs.y;
ret.z = z+rhs.z;
return ret;
}
__host__ __device__ cuvect3 cuvect3::operator-(const cuvect3 &rhs)
{
cuvect3 ret;
ret.x = x-rhs.x;
ret.y = y-rhs.y;
ret.z = z-rhs.z;
return ret;
}
__host__ __device__ cuvect3 cuvect3::operator*(const double &rhs)
{
cuvect3 ret;
ret.x = x*rhs;
ret.y = y*rhs;
ret.z = z*rhs;
return ret;
}
__host__ __device__ cuvect3 cuvect3::operator/(const double &rhs)
{
cuvect3 ret;
ret.x = x/rhs;
ret.y = y/rhs;
ret.z = z/rhs;
return ret;
}
__host__ __device__ cuvect3& cuvect3::operator+=(const cuvect3 &rhs)
{
x = x + rhs.x;
y = y + rhs.y;
z = z + rhs.z;
return *this;
}
__host__ __device__ cuvect3& cuvect3::operator-=(const cuvect3 &rhs)
{
x = x - rhs.x;
y = y - rhs.y;
z = z - rhs.z;
return *this;
}
__host__ __device__ cuvect3& cuvect3::operator*=(const double &rhs)
{
x = x * rhs;
y = y * rhs;
z = z * rhs;
return *this;
}
__host__ __device__ cuvect3& cuvect3::operator/=(const double &rhs)
{
x = x / rhs;
y = y / rhs;
z = z / rhs;
return *this;
}
__host__ __device__ cuvect3::cuvect3(const double &_x, const double &_y, const double &_z)
{
x = _x; y = _y; z = _z;
return;
}
__host__ __device__ double cuvect3_dot(const cuvect3 &a, const cuvect3 &b)
{
double ret = a.x*b.x+a.y*b.y+a.z*b.z;
return ret;
}
__host__ __device__ cuvect3 cuvect3_cross(const cuvect3 &a, const cuvect3 &b)
{
cuvect3 ret;
ret[0] = a[1]*b[2]-a[2]*b[1];
ret[1] = a[2]*b[0]-a[0]*b[2];
ret[2] = a[0]*b[1]-a[1]*b[0];
return ret;
}
__host__ __device__ double cuvect3_norm(const cuvect3 &a)
{
double ret;
ret = ::sqrtf(a.x*a.x+a.y*a.y+a.z*a.z);
return ret;
}
__host__ __device__ cuvect3 cuvect3_normalize(const cuvect3 &a)
{
cuvect3 ret;
double m;
m = ::sqrtf(a.x*a.x+a.y*a.y+a.z*a.z);
if(m>0.0)
{
ret.x = a.x/m; ret.y = a.y/m; ret.z = a.z/m;
}
else
{
ret.x = 0.0; ret.y = 0.0; ret.z = 0.0;
}
return ret;
}
__host__ __device__ cuvect3 cuvect3_proj(const cuvect3 &a, const cuvect3 &b)
{
cuvect3 ret;
cuvect3 bn = cuvect3_normalize(b);
double m = cuvect3_dot(a,bn);
ret = bn*m;
return ret;
}
__host__ __device__ cumat3::cumat3()
{
m00 = 0.0;
m01 = 0.0;
m02 = 0.0;
m10 = 0.0;
m11 = 0.0;
m12 = 0.0;
m20 = 0.0;
m21 = 0.0;
m22 = 0.0;
return;
}
__host__ __device__ cumat3::~cumat3()
{
m00 = 0.0;
m01 = 0.0;
m02 = 0.0;
m10 = 0.0;
m11 = 0.0;
m12 = 0.0;
m20 = 0.0;
m21 = 0.0;
m22 = 0.0;
return;
}
__host__ __device__ double& cumat3::operator[](const int &I)
{
if(I==0) return m00;
if(I==1) return m10;
if(I==2) return m20;
if(I==3) return m01;
if(I==4) return m11;
if(I==5) return m21;
if(I==6) return m02;
if(I==7) return m12;
if(I==8) return m22;
return m00;
}
__host__ __device__ double& cumat3::operator()(const int &I, const int &J)
{
if(I==0 && J==0) return m00;
if(I==1 && J==0) return m10;
if(I==2 && J==0) return m20;
if(I==0 && J==1) return m01;
if(I==1 && J==1) return m11;
if(I==2 && J==1) return m21;
if(I==0 && J==2) return m02;
if(I==1 && J==2) return m12;
if(I==2 && J==2) return m22;
return m00;
}
__host__ __device__ double& cumat3::at(const int &I, const int &J)
{
if(I==0 && J==0) return m00;
if(I==1 && J==0) return m10;
if(I==2 && J==0) return m20;
if(I==0 && J==1) return m01;
if(I==1 && J==1) return m11;
if(I==2 && J==1) return m21;
if(I==0 && J==2) return m02;
if(I==1 && J==2) return m12;
if(I==2 && J==2) return m22;
return m00;
}
__host__ __device__ const double& cumat3::operator[](const int &I) const
{
if(I==0) return m00;
if(I==1) return m10;
if(I==2) return m20;
if(I==3) return m01;
if(I==4) return m11;
if(I==5) return m21;
if(I==6) return m02;
if(I==7) return m12;
if(I==8) return m22;
return m00;
}
__host__ __device__ const double& cumat3::operator()(const int &I, const int &J) const
{
if(I==0 && J==0) return m00;
if(I==1 && J==0) return m10;
if(I==2 && J==0) return m20;
if(I==0 && J==1) return m01;
if(I==1 && J==1) return m11;
if(I==2 && J==1) return m21;
if(I==0 && J==2) return m02;
if(I==1 && J==2) return m12;
if(I==2 && J==2) return m22;
return m00;
}
__host__ __device__ const double& cumat3::at(const int &I, const int &J) const
{
if(I==0 && J==0) return m00;
if(I==1 && J==0) return m10;
if(I==2 && J==0) return m20;
if(I==0 && J==1) return m01;
if(I==1 && J==1) return m11;
if(I==2 && J==1) return m21;
if(I==0 && J==2) return m02;
if(I==1 && J==2) return m12;
if(I==2 && J==2) return m22;
return m00;
}
__host__ __device__ cumat3 cumat3::operator+(const cumat3 &rhs)
{
cumat3 ret;
ret.m00 = m00 + rhs.m00;
ret.m10 = m10 + rhs.m10;
ret.m20 = m20 + rhs.m20;
ret.m01 = m01 + rhs.m01;
ret.m11 = m11 + rhs.m11;
ret.m21 = m21 + rhs.m21;
ret.m02 = m02 + rhs.m02;
ret.m12 = m12 + rhs.m12;
ret.m22 = m22 + rhs.m22;
return ret;
}
__host__ __device__ cumat3 cumat3::operator-(const cumat3 &rhs)
{
cumat3 ret;
ret.m00 = m00 - rhs.m00;
ret.m10 = m10 - rhs.m10;
ret.m20 = m20 - rhs.m20;
ret.m01 = m01 - rhs.m01;
ret.m11 = m11 - rhs.m11;
ret.m21 = m21 - rhs.m21;
ret.m02 = m02 - rhs.m02;
ret.m12 = m12 - rhs.m12;
ret.m22 = m22 - rhs.m22;
return ret;
}
__host__ __device__ cumat3 cumat3::operator*(const double &rhs)
{
cumat3 ret;
ret.m00 = m00 * rhs;
ret.m10 = m10 * rhs;
ret.m20 = m20 * rhs;
ret.m01 = m01 * rhs;
ret.m11 = m11 * rhs;
ret.m21 = m21 * rhs;
ret.m02 = m02 * rhs;
ret.m12 = m12 * rhs;
ret.m22 = m22 * rhs;
return ret;
}
__host__ __device__ cumat3 cumat3::operator/(const double &rhs)
{
cumat3 ret;
ret.m00 = m00 / rhs;
ret.m10 = m10 / rhs;
ret.m20 = m20 / rhs;
ret.m01 = m01 / rhs;
ret.m11 = m11 / rhs;
ret.m21 = m21 / rhs;
ret.m02 = m02 / rhs;
ret.m12 = m12 / rhs;
ret.m22 = m22 / rhs;
return ret;
}
__host__ __device__ cuvect3 cumat3::operator*(const cuvect3 &rhs)
{
cuvect3 ret;
ret.x = m00*rhs.x + m01*rhs.y + m02*rhs.z;
ret.y = m10*rhs.x + m11*rhs.y + m12*rhs.z;
ret.z = m20*rhs.x + m21*rhs.y + m22*rhs.z;
return ret;
}
__host__ __device__ cumat3 cumat3::operator*(const cumat3 &rhs)
{
cumat3 ret;
ret.m00 = m00*rhs.m00 + m01*rhs.m10 + m02*rhs.m20;
ret.m01 = m00*rhs.m01 + m01*rhs.m11 + m02*rhs.m21;
ret.m02 = m00*rhs.m02 + m01*rhs.m12 + m02*rhs.m22;
ret.m10 = m10*rhs.m00 + m11*rhs.m10 + m12*rhs.m20;
ret.m11 = m10*rhs.m01 + m11*rhs.m11 + m12*rhs.m21;
ret.m12 = m10*rhs.m02 + m11*rhs.m12 + m12*rhs.m22;
ret.m20 = m20*rhs.m00 + m21*rhs.m10 + m22*rhs.m20;
ret.m21 = m20*rhs.m01 + m21*rhs.m11 + m22*rhs.m21;
ret.m22 = m20*rhs.m02 + m21*rhs.m12 + m22*rhs.m22;
return ret;
}
__host__ __device__ double cumat3::det()
{
double ret = 0.0;
ret += m00*m11*m22;
ret += m01*m12*m20;
ret += m02*m10*m21;
ret -= m20*m11*m02;
ret -= m21*m12*m00;
ret -= m22*m10*m01;
return ret;
}
__host__ __device__ cumat3 cumat3::transpose()
{
cumat3 ret;
ret.m00 = m00;
ret.m01 = m10;
ret.m02 = m20;
ret.m10 = m01;
ret.m11 = m11;
ret.m12 = m21;
ret.m20 = m02;
ret.m21 = m12;
ret.m22 = m22;
return ret;
}
__host__ __device__ cumat3 cumat3::inverse()
{
cumat3 q;
double dt = det();
if(dt!=0.0)
{
q(0,0) = (at(1,1)*at(2,2)-at(1,2)*at(2,1))/dt;
q(1,0) = -(at(1,0)*at(2,2)-at(1,2)*at(2,0))/dt;
q(2,0) = (at(1,0)*at(2,1)-at(1,1)*at(2,0))/dt;
q(0,1) = -(at(0,1)*at(2,2)-at(0,2)*at(2,1))/dt;
q(1,1) = (at(0,0)*at(2,2)-at(0,2)*at(2,0))/dt;
q(2,1) = -(at(0,0)*at(2,1)-at(0,1)*at(2,0))/dt;
q(0,2) = (at(0,1)*at(1,2)-at(0,2)*at(1,1))/dt;
q(1,2) = -(at(0,0)*at(1,2)-at(0,2)*at(1,0))/dt;
q(2,2) = (at(0,0)*at(1,1)-at(0,1)*at(1,0))/dt;
// q(0,0) = (at(1,1)*at(2,2)-at(1,2)*at(2,1))/dt;
// q(0,1) = -(at(1,0)*at(2,2)-at(1,2)*at(2,0))/dt;
// q(0,2) = (at(1,0)*at(2,1)-at(1,1)*at(2,0))/dt;
// q(1,0) = -(at(0,1)*at(2,2)-at(0,2)*at(2,1))/dt;
// q(1,1) = (at(0,0)*at(2,2)-at(0,2)*at(2,0))/dt;
// q(1,2) = -(at(0,0)*at(2,1)-at(0,1)*at(2,0))/dt;
// q(2,0) = (at(0,1)*at(1,2)-at(0,2)*at(1,1))/dt;
// q(2,1) = -(at(0,0)*at(1,2)-at(0,2)*at(1,0))/dt;
// q(2,2) = (at(0,0)*at(1,1)-at(0,1)*at(1,0))/dt;
// q = q.transpose();
}
else
{
q(0,0) = inf;
q(0,1) = inf;
q(0,2) = inf;
q(1,0) = inf;
q(1,1) = inf;
q(1,2) = inf;
q(2,0) = inf;
q(2,1) = inf;
q(2,2) = inf;
}
return q;
}
__host__ __device__ cumat3 operator-(const cumat3 &rhs)
{
cumat3 ret;
ret.m00 = -rhs.m00;
ret.m10 = -rhs.m10;
ret.m20 = -rhs.m20;
ret.m01 = -rhs.m01;
ret.m11 = -rhs.m11;
ret.m21 = -rhs.m21;
ret.m02 = -rhs.m02;
ret.m12 = -rhs.m12;
ret.m22 = -rhs.m22;
return ret;
}
__host__ __device__ cumat3& cumat3::operator+=(const cumat3 &rhs)
{
m00 += rhs.m00;
m10 += rhs.m10;
m20 += rhs.m20;
m01 += rhs.m01;
m11 += rhs.m11;
m21 += rhs.m21;
m02 += rhs.m02;
m12 += rhs.m12;
m22 += rhs.m22;
return *this;
}
__host__ __device__ cumat3& cumat3::operator-=(const cumat3 &rhs)
{
m00 -= rhs.m00;
m10 -= rhs.m10;
m20 -= rhs.m20;
m01 -= rhs.m01;
m11 -= rhs.m11;
m21 -= rhs.m21;
m02 -= rhs.m02;
m12 -= rhs.m12;
m22 -= rhs.m22;
return *this;
}
__host__ __device__ cumat3& cumat3::operator/=(const double &rhs)
{
m00 /= rhs;
m10 /= rhs;
m20 /= rhs;
m01 /= rhs;
m11 /= rhs;
m21 /= rhs;
m02 /= rhs;
m12 /= rhs;
m22 /= rhs;
return *this;
}
__host__ __device__ cumat3& cumat3::operator*=(const double &rhs)
{
m00 *= rhs;
m10 *= rhs;
m20 *= rhs;
m01 *= rhs;
m11 *= rhs;
m21 *= rhs;
m02 *= rhs;
m12 *= rhs;
m22 *= rhs;
return *this;
}
__host__ __device__ cumat3& cumat3::operator*=(const cumat3 &rhs)
{
cumat3 tmp;
tmp.m00 = m00*rhs.m00 + m01*rhs.m10 + m02*rhs.m20;
tmp.m01 = m00*rhs.m01 + m01*rhs.m11 + m02*rhs.m21;
tmp.m02 = m00*rhs.m02 + m01*rhs.m12 + m02*rhs.m22;
tmp.m10 = m10*rhs.m00 + m11*rhs.m10 + m12*rhs.m20;
tmp.m11 = m10*rhs.m01 + m11*rhs.m11 + m12*rhs.m21;
tmp.m12 = m10*rhs.m02 + m11*rhs.m12 + m12*rhs.m22;
tmp.m20 = m20*rhs.m00 + m21*rhs.m10 + m22*rhs.m20;
tmp.m21 = m20*rhs.m01 + m21*rhs.m11 + m22*rhs.m21;
tmp.m22 = m20*rhs.m02 + m21*rhs.m12 + m22*rhs.m22;
(*this) = tmp;
return *this;
}
__host__ __device__ cumat3::cumat3(
const double & _m00, const double & _m01, const double & _m02,
const double & _m10, const double & _m11, const double & _m12,
const double & _m20, const double & _m21, const double & _m22
)
{
m00 = _m00;
m01 = _m01;
m02 = _m02;
m10 = _m10;
m11 = _m11;
m12 = _m12;
m20 = _m20;
m21 = _m21;
m22 = _m22;
}
__host__ __device__ double* cumat3::data()
{
//pointer to double[9] representation of matrix
return (double*) this;
}
__host__ __device__ const double* cumat3::data() const
{
//pointer to double[9] representation of matrix
return (const double*) this;
}
//transposes a 3x3 (9 element) matrix
__host__ __device__ void mat3f_transpose(double *mat3inout)
{
int I,J;
double matint[9];
for(I=0;I<9;I++)
{
matint[I] = mat3inout[I];
}
for(I=0;I<3;I++)
{
for(J=0;J<3;J++)
{
mat3inout[I+J*3] = matint[J+I*3];
}
}
return;
}
//copies src to dest
__host__ __device__ void mat3f_copy(double *mat3f_dest, const double *mat3f_src)
{
int I;
if(mat3f_dest==NULL || mat3f_src==NULL)
return;
for(I=0;I<9;I++)
mat3f_dest[I] = mat3f_src[I];
return;
}
__host__ __device__ double mat3f_det(double *mat3in)
{
double ret = 0.0;
ret = ret + mat3in[0+0*3]*mat3in[1+1*3]*mat3in[2+2*3];
ret = ret + mat3in[0+1*3]*mat3in[1+2*3]*mat3in[2+0*3];
ret = ret + mat3in[0+2*3]*mat3in[1+0*3]*mat3in[2+1*3];
ret = ret - mat3in[0+0*3]*mat3in[1+2*3]*mat3in[2+1*3];
ret = ret - mat3in[0+1*3]*mat3in[1+0*3]*mat3in[2+2*3];
ret = ret - mat3in[0+2*3]*mat3in[1+1*3]*mat3in[2+0*3];
return ret;
}
//inverts a 3x3 (9 element) matrix
__host__ __device__ void mat3f_inverse(double *mat3inout)
{
int I;
double matint[9];
double det = mat3f_det(mat3inout);
for(I=0;I<9;I++)
{
matint[I] = mat3inout[I];
}
mat3inout[0+0*3] = (matint[1+1*3]*matint[2+2*3]-matint[1+2*3]*matint[2+1*3])/det;
mat3inout[0+1*3] = -(matint[1+0*3]*matint[2+2*3]-matint[1+2*3]*matint[2+0*3])/det;
mat3inout[0+2*3] = (matint[1+0*3]*matint[2+1*3]-matint[1+1*3]*matint[2+0*3])/det;
mat3inout[1+0*3] = -(matint[0+1*3]*matint[2+2*3]-matint[0+2*3]*matint[2+1*3])/det;
mat3inout[1+1*3] = (matint[0+0*3]*matint[2+2*3]-matint[0+2*3]*matint[2+0*3])/det;
mat3inout[1+2*3] = -(matint[0+0*3]*matint[2+1*3]-matint[0+1*3]*matint[2+0*3])/det;
mat3inout[2+0*3] = (matint[0+1*3]*matint[1+2*3]-matint[0+2*3]*matint[1+1*3])/det;
mat3inout[2+1*3] = -(matint[0+0*3]*matint[1+2*3]-matint[0+2*3]*matint[1+0*3])/det;
mat3inout[2+2*3] = (matint[0+0*3]*matint[1+1*3]-matint[0+1*3]*matint[1+0*3])/det;
mat3f_transpose(mat3inout);
return;
}
__host__ __device__ cuvect3 mat3f_mult(double *mat3in, const cuvect3 &cvin)
{
int I,J;
cuvect3 ret;
for(I=0;I<3;I++)
{
ret[I] = 0.0;
for(J=0;J<3;J++)
{
ret[I] = ret[I] + mat3in[I+3*J]*cvin[J];
}
}
return ret;
}
__host__ __device__ void mat3f_mult(double *matina, double *matinb, double *matout)
{
double wrk[9];
int I,J,K;
for(I=0;I<3;I++)
{
for(J=0;J<3;J++)
{
wrk[I+3*J] = 0.0;
}
}
for(I=0;I<3;I++)
{
for(J=0;J<3;J++)
{
for(K=0;K<3;K++)
{
wrk[I+3*K] = wrk[I+3*K] + matina[I+3*J]*matinb[J+3*K];
}
}
}
for(I=0;I<3;I++)
{
for(J=0;J<3;J++)
{
matout[I+3*J] = wrk[I+3*J];
}
}
return;
}
__host__ void test_cudavect_logic1()
{
//3 dim vector and matrix functional tests on host side
// printf("3 dim vector and matrix functional tests on host side\n");
// cuvect3 a,b,c;
// double ma[9],mb[9],mc[9];
// int I,J;
// for(I=0;I<3;I++)
// {
// for(J=0;J<3;J++)
// {
// ma[I+3*J] = ((double) rand())/((double) RAND_MAX);
// mb[I+3*J] = ma[I+3*J];
// }
// }
// mat3f_inverse(mb);
// mat3f_mult(ma,mb,mc);
// for(I=0;I<3;I++)
// {
// for(J=0;J<3;J++)
// {
// printf("ma[%d,%d] = %1.3f\n",I,J,ma[I+3*J]);
// }
// }
// for(I=0;I<3;I++)
// {
// for(J=0;J<3;J++)
// {
// printf("mb[%d,%d] = %1.3f\n",I,J,mb[I+3*J]);
// }
// }
// for(I=0;I<3;I++)
// {
// for(J=0;J<3;J++)
// {
// printf("mc[%d,%d] = %1.3f\n",I,J,mc[I+3*J]);
// }
// }
// a = cuvect3(1,1,1);
// b = mat3f_mult(ma,a);
// b = mat3f_mult(mb,b);
// for(I=0;I<3;I++)
// {
// printf("a[%d] = %1.3f, b[%d] = %1.3f\n",I,a[I],I,b[I]);
// }
// a = cuvect3(1,0,1);
// b = cuvect3(0,1,-1);
// c = a+b;
// for(I=0;I<3;I++)
// {
// printf("a[%d] = %1.3f, b[%d] = %1.3f, c[%d] = %1.3f\n",I,a[I],I,b[I],I,c[I]);
// }
// c = c/2.0;
// for(I=0;I<3;I++)
// {
// printf("a[%d] = %1.3f, b[%d] = %1.3f, c[%d] = %1.3f\n",I,a[I],I,b[I],I,c[I]);
// }
// c = cuvect3_cross(a,b);
// for(I=0;I<3;I++)
// {
// printf("a[%d] = %1.3f, b[%d] = %1.3f, c[%d] = %1.3f\n",I,a[I],I,b[I],I,c[I]);
// }
// printf("c dot a = %1.3f, c dot b = %1.3f\n",cuvect3_dot(c,a),cuvect3_dot(c,b));
// printf("norm(a)=%1.3f, norm(b)=%1.3f, norm(c)=%1.3f\n",cuvect3_norm(a),cuvect3_norm(b),cuvect3_norm(c));
// a = cuvect3_normalize(a);
// b = cuvect3_normalize(b);
// c = cuvect3_normalize(c);
// for(I=0;I<3;I++)
// {
// printf("a[%d] = %1.3f, b[%d] = %1.3f, c[%d] = %1.3f\n",I,a[I],I,b[I],I,c[I]);
// }
// printf("c dot a = %1.3f, c dot b = %1.3f\n",cuvect3_dot(c,a),cuvect3_dot(c,b));
// printf("norm(a)=%1.3f, norm(b)=%1.3f, norm(c)=%1.3f\n",cuvect3_norm(a),cuvect3_norm(b),cuvect3_norm(c));
return;
}
__host__ __device__ cumat3 hodge_dual(const cuvect3 &vin)
{
cumat3 ret;
ret.m00 = 0.0;
ret.m11 = 0.0;
ret.m22 = 0.0;
ret.m01 = vin.z;
ret.m12 = vin.x;
ret.m20 = vin.y;
ret.m10 = -vin.z;
ret.m21 = -vin.x;
ret.m02 = -vin.y;
return ret;
}
__host__ __device__ cuvect3 hodge_dual(const cumat3 &min)
{
cuvect3 ret;
ret.x = 0.5f*(min.m12 - min.m21);
ret.y = 0.5f*(min.m20 - min.m02);
ret.z = 0.5f*(min.m01 - min.m10);
return ret;
}
__host__ __device__ const cumat3 cumat3_eye()
{
cumat3 ret;
ret.m00 = 1.0f;
ret.m11 = 1.0f;
ret.m22 = 1.0f;
return ret;
}
__host__ __device__ const cumat3 cumat3_zeros()
{
cumat3 ret;
return ret;
}
__host__ __device__ cumat3 rotmat_from_axisangle(const cuvect3 &axis, const double &angle)
{
cumat3 ret = cumat3_zeros();
cumat3 H;
cuvect3 _naxis;
_naxis = cuvect3_normalize(axis);
H = hodge_dual(_naxis);
ret += H*sinf(angle);
H = H*H;
ret -= H*cosf(angle);
ret += (H + cumat3_eye());
return ret;
}
};

432
src/amsculib3/math/cuvect4.cu
View File

@ -1,432 +0,0 @@
#include <amsculib3/amsculib3.hpp>
namespace amscuda
{
__host__ __device__ cuvect4::cuvect4()
{
x = 0.0; y = 0.0; z = 0.0; w = 0.0;
return;
}
__host__ __device__ cuvect4::~cuvect4()
{
x = 0.0; y = 0.0; z = 0.0; w = 0.0;
return;
}
__host__ __device__ cuvect4::cuvect4(double _x, double _y, double _z, double _w)
{
x = _x; y = _y; z = _z; w = _w;
return;
}
__host__ __device__ double& cuvect4::operator[](const int I)
{
if(I==0) return x;
else if(I==1) return y;
else if(I==2) return z;
else if(I==3) return w;
return x;
}
__host__ __device__ const double& cuvect4::operator[](const int I) const
{
if(I==0) return x;
else if(I==1) return y;
else if(I==2) return z;
else if(I==3) return w;
return x;
}
__host__ __device__ cuvect4 cuvect4::operator+(cuvect4 lhs)
{
cuvect4 ret;
ret.x = this->x + lhs.x;
ret.y = this->y + lhs.y;
ret.z = this->z + lhs.z;
ret.w = this->w + lhs.w;
return ret;
}
__host__ __device__ cuvect4 cuvect4::operator-(cuvect4 lhs)
{
cuvect4 ret;
ret.x = this->x - lhs.x;
ret.y = this->y - lhs.y;
ret.z = this->z - lhs.z;
ret.w = this->w - lhs.w;
return ret;
}
__host__ __device__ cuvect4 cuvect4::operator*(double lhs)
{
cuvect4 ret;
ret.x = this->x*lhs;
ret.y = this->y*lhs;
ret.z = this->z*lhs;
ret.w = this->w*lhs;
return ret;
}
__host__ __device__ cuvect4 cuvect4::operator/(double lhs)
{
cuvect4 ret;
ret.x = this->x/lhs;
ret.y = this->y/lhs;
ret.z = this->z/lhs;
ret.w = this->w/lhs;
return ret;
}
__host__ __device__ cumat4::cumat4()
{
int I;
for(I=0;I<16;I++)
{
dat[I] = 0.0;
}
return;
}
__host__ __device__ cumat4::~cumat4()
{
int I;
for(I=0;I<16;I++)
{
dat[I] = 0.0;
}
return;
}
__host__ __device__ double& cumat4::operator[](const int I)
{
return dat[I];
}
__host__ __device__ double& cumat4::operator()(const int I, const int J)
{
return dat[I+4*J];
}
__host__ __device__ double& cumat4::at(const int I, const int J)
{
return dat[I+4*J];
}
__host__ __device__ cumat4 cumat4::operator+(cumat4 lhs)
{
cumat4 ret;
int I;
for(I=0;I<16;I++)
{
ret.dat[I] = this->dat[I] + lhs.dat[I];
}
return ret;
}
__host__ __device__ cumat4 cumat4::operator-(cumat4 lhs)
{
cumat4 ret;
int I;
for(I=0;I<16;I++)
{
ret.dat[I] = this->dat[I] - lhs.dat[I];
}
return ret;
}
__host__ __device__ cumat4 cumat4::operator*(double lhs)
{
cumat4 ret;
int I;
for(I=0;I<16;I++)
{
ret.dat[I] = this->dat[I]*lhs;
}
return ret;
}
__host__ __device__ cumat4 cumat4::operator/(double lhs)
{
cumat4 ret;
int I;
for(I=0;I<16;I++)
{
ret.dat[I] = this->dat[I]/lhs;
}
return ret;
}
__host__ __device__ cuvect4 cumat4::operator*(cuvect4 lhs)
{
cuvect4 ret = cuvect4(0.0,0.0,0.0,0.0);
int I,J;
for(I=0;I<4;I++)
{
for(J=0;J<4;J++)
{
ret[I] = ret[I] + this->at(I,J)*lhs[J];
}
}
return ret;
}
__host__ __device__ cumat4 cumat4::operator*(cumat4 lhs)
{
cumat4 ret;
int I,J,K;
for(I=0;I<4;I++)
{
for(J=0;J<4;J++)
{
ret(I,J) = 0;
for(K=0;K<4;K++)
{
ret(I,J) = ret(I,J) + this->at(I,K) * lhs(K,J);
}
}
}
return ret;
}
__host__ __device__ cumat4 cumat4::transpose()
{
cumat4 q;
int I,J;
for(I=0;I<4;I++)
{
for(J=0;J<4;J++)
{
q(I,J) = this->at(J,I);
}
}
return q;
}
__host__ __device__ double cumat4::det()
{
double a00,a01,a02,a03;
double a10,a11,a12,a13;
double a20,a21,a22,a23;
double a30,a31,a32,a33;
double det;
a00 = this->at(0,0);
a01 = this->at(0,1);
a02 = this->at(0,2);
a03 = this->at(0,3);
a10 = this->at(1,0);
a11 = this->at(1,1);
a12 = this->at(1,2);
a13 = this->at(1,3);
a20 = this->at(2,0);
a21 = this->at(2,1);
a22 = this->at(2,2);
a23 = this->at(2,3);
a30 = this->at(3,0);
a31 = this->at(3,1);
a32 = this->at(3,2);
a33 = this->at(3,3);
det = a03*a12*a21*a30 -
a02*a13*a21*a30 -
a03*a11*a22*a30 +
a01*a13*a22*a30 +
a02*a11*a23*a30 -
a01*a12*a23*a30 -
a03*a12*a20*a31 +
a02*a13*a20*a31 +
a03*a10*a22*a31 -
a00*a13*a22*a31 -
a02*a10*a23*a31 +
a00*a12*a23*a31 +
a03*a11*a20*a32 -
a01*a13*a20*a32 -
a03*a10*a21*a32 +
a00*a13*a21*a32 +
a01*a10*a23*a32 -
a00*a11*a23*a32 -
a02*a11*a20*a33 +
a01*a12*a20*a33 +
a02*a10*a21*a33 -
a00*a12*a21*a33 -
a01*a10*a22*a33 +
a00*a11*a22*a33;
return det;
}
__host__ __device__ cumat4 minverse(cumat4 ma)
{
cumat4 mb;
double a00,a01,a02,a03;
double a10,a11,a12,a13;
double a20,a21,a22,a23;
double a30,a31,a32,a33;
double b00,b01,b02,b03;
double b10,b11,b12,b13;
double b20,b21,b22,b23;
double b30,b31,b32,b33;
double det = 0.0;
a00 = ma.at(0,0);
a01 = ma.at(0,1);
a02 = ma.at(0,2);
a03 = ma.at(0,3);
a10 = ma.at(1,0);
a11 = ma.at(1,1);
a12 = ma.at(1,2);
a13 = ma.at(1,3);
a20 = ma.at(2,0);
a21 = ma.at(2,1);
a22 = ma.at(2,2);
a23 = ma.at(2,3);
a30 = ma.at(3,0);
a31 = ma.at(3,1);
a32 = ma.at(3,2);
a33 = ma.at(3,3);
det = a03*a12*a21*a30 -
a02*a13*a21*a30 -
a03*a11*a22*a30 +
a01*a13*a22*a30 +
a02*a11*a23*a30 -
a01*a12*a23*a30 -
a03*a12*a20*a31 +
a02*a13*a20*a31 +
a03*a10*a22*a31 -
a00*a13*a22*a31 -
a02*a10*a23*a31 +
a00*a12*a23*a31 +
a03*a11*a20*a32 -
a01*a13*a20*a32 -
a03*a10*a21*a32 +
a00*a13*a21*a32 +
a01*a10*a23*a32 -
a00*a11*a23*a32 -
a02*a11*a20*a33 +
a01*a12*a20*a33 +
a02*a10*a21*a33 -
a00*a12*a21*a33 -
a01*a10*a22*a33 +
a00*a11*a22*a33;
if(det*det>1.0E-30)
{
b00 = -a13*a22*a31 + a12*a23*a31 + a13*a21*a32 - a11*a23*a32 - a12*a21*a33 + a11*a22*a33;
b01 = a03*a22*a31 - a02*a23*a31 - a03*a21*a32 + a01*a23*a32 + a02*a21*a33 - a01*a22*a33;
b02 = -a03*a12*a31 + a02*a13*a31 + a03*a11*a32 - a01*a13*a32 - a02*a11*a33 + a01*a12*a33;
b03 = a03*a12*a21 - a02*a13*a21 - a03*a11*a22 + a01*a13*a22 + a02*a11*a23 - a01*a12*a23;
b10 = a13*a22*a30 - a12*a23*a30 - a13*a20*a32 + a10*a23*a32 + a12*a20*a33 - a10*a22*a33;
b11 = -a03*a22*a30 + a02*a23*a30 + a03*a20*a32 - a00*a23*a32 - a02*a20*a33 + a00*a22*a33;
b12 = a03*a12*a30 - a02*a13*a30 - a03*a10*a32 + a00*a13*a32 + a02*a10*a33 - a00*a12*a33;
b13 = -a03*a12*a20 + a02*a13*a20 + a03*a10*a22 - a00*a13*a22 - a02*a10*a23 + a00*a12*a23;
b20 = -a13*a21*a30 + a11*a23*a30 + a13*a20*a31 - a10*a23*a31 - a11*a20*a33 + a10*a21*a33;
b21 = a03*a21*a30 - a01*a23*a30 - a03*a20*a31 + a00*a23*a31 + a01*a20*a33 - a00*a21*a33;
b22 = -a03*a11*a30 + a01*a13*a30 + a03*a10*a31 - a00*a13*a31 - a01*a10*a33 + a00*a11*a33;
b23 = a03*a11*a20 - a01*a13*a20 - a03*a10*a21 + a00*a13*a21 + a01*a10*a23 - a00*a11*a23;
b30 = a12*a21*a30 - a11*a22*a30 - a12*a20*a31 + a10*a22*a31 + a11*a20*a32 - a10*a21*a32;
b31 = -a02*a21*a30 + a01*a22*a30 + a02*a20*a31 - a00*a22*a31 - a01*a20*a32 + a00*a21*a32;
b32 = a02*a11*a30 - a01*a12*a30 - a02*a10*a31 + a00*a12*a31 + a01*a10*a32 - a00*a11*a32;
b33 = -a02*a11*a20 + a01*a12*a20 + a02*a10*a21 - a00*a12*a21 - a01*a10*a22 + a00*a11*a22;
b00 = b00/det;
b01 = b01/det;
b02 = b02/det;
b03 = b03/det;
b10 = b10/det;
b11 = b11/det;
b12 = b12/det;
b13 = b13/det;
b20 = b20/det;
b21 = b21/det;
b22 = b22/det;
b23 = b23/det;
b30 = b30/det;
b31 = b31/det;
b32 = b32/det;
b33 = b33/det;
mb.at(0,0) = b00;
mb.at(0,1) = b01;
mb.at(0,2) = b02;
mb.at(0,3) = b03;
mb.at(1,0) = b10;
mb.at(1,1) = b11;
mb.at(1,2) = b12;
mb.at(1,3) = b13;
mb.at(2,0) = b20;
mb.at(2,1) = b21;
mb.at(2,2) = b22;
mb.at(2,3) = b23;
mb.at(3,0) = b30;
mb.at(3,1) = b31;
mb.at(3,2) = b32;
mb.at(3,3) = b33;
}
//this was STUPID. Gaah. Computer algebra system saves the day? I'd be surprised if this didn't end up *slower* than gaussian elimination. Don't do this again!
return mb;
}
__host__ __device__ cumat4 cumat4::inverse()
{
return minverse(*this);
}
__host__ __device__ double cuvect4_dot(cuvect4 a, cuvect4 b)
{
double ret = 0.0;
ret = a.x*b.x + a.y*b.y + a.z*b.z + a.w*b.w;
return ret;
}
__host__ __device__ double cuvect4_norm(cuvect4 a)
{
double ret = 0.0;
ret = ::sqrt(cuvect4_dot(a,a));
return ret;
}
__host__ __device__ cuvect4 cuvect4_normalize(cuvect4 a)
{
cuvect4 ret = cuvect4(0.0f,0.0f,0.0f,0.0f);
double nrm = cuvect4_norm(a);
if(nrm>0.0)
ret = a/nrm;
return ret;
}
__host__ __device__ cuvect4 cuvect4_proj(cuvect4 a, cuvect4 b)
{
cuvect4 ret;
cuvect4 bn = cuvect4_normalize(b);
double d = cuvect4_dot(a,bn);
ret = bn*d;
return ret;
}
__host__ __device__ cuvect4 operator-(cuvect4 rhs)
{
cuvect4 ret;
ret[0] = -rhs[0];
ret[1] = -rhs[1];
ret[2] = -rhs[2];
ret[3] = -rhs[3];
return ret;
}
__host__ __device__ cumat4 operator-(cumat4 rhs)
{
cumat4 ret;
int I;
for(I=0;I<16;I++) ret[I] = -rhs[I];
return ret;
}
};