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

Making some of my library code public.

Browse Source
This commit is contained in:
Aaron
2025-02-04 21:59:29 -05:00
commit 64a7bfc851
85 changed files with 9729 additions and 0 deletions
Download Patch File Download Diff File Expand all files Collapse all files

476
src/amsculib2/amscu_comp128.cu Normal file
View File

@ -0,0 +1,476 @@
#include <amsculib2/amsculib2.hpp>
namespace amscuda
{
namespace cmp
{
__host__ __device__ cucomp128::cucomp128()
{
real = 0.0;
imag = 0.0;
return;
}
__host__ __device__ cucomp128::~cucomp128()
{
real = 0.0;
imag = 0.0;
return;
}
__host__ __device__ cucomp128::cucomp128(const cucomp128 &other)
{
real = other.real;
imag = other.imag;
return;
}
__host__ __device__ cucomp128::cucomp128(const double &other)
{
real = other;
imag = 0.0;
return;
}
__host__ __device__ cucomp128& cucomp128::operator=(cucomp128& other)
{
real = other.real;
imag = other.imag;
return *this;
}
__host__ __device__ const cucomp128& cucomp128::operator=(const cucomp128& other)
{
this->real = other.real;
this->imag = other.imag;
return *this;
}
__host__ __device__ cucomp128& cucomp128::operator=(double& other)
{
real = other;
imag = 0.0;
return *this;
}
__host__ __device__ const cucomp128& cucomp128::operator=(const double& other)
{
this->real = other;
this->imag = 0.0;
return *this;
}
__host__ __device__ double& cucomp128::operator[](int& ind)
{
if(ind==0)
{
return this->real;
}
else
{
return this->imag;
}
}
__host__ __device__ const double& cucomp128::operator[](const int& ind) const
{
if(ind==0)
{
return this->real;
}
else
{
return this->imag;
}
}
__host__ __device__ cucomp128 cucomp128::operator+(const cucomp128& z)
{
cucomp128 ret;
ret.real = real + z.real;
ret.imag = imag + z.imag;
return ret;
}
__host__ __device__ cucomp128 cucomp128::operator-(const cucomp128& z)
{
cucomp128 ret;
ret.real = real - z.real;
ret.imag = imag - z.imag;
return ret;
}
__host__ __device__ cucomp128 cucomp128::operator*(const cucomp128& z)
{
cucomp128 ret;
ret.real = (real*z.real - imag*z.imag);
ret.imag = (imag*z.real + real*z.imag);
return ret;
}
__host__ __device__ cucomp128 cucomp128::operator/(const cucomp128& z)
{
cucomp128 ret;
double zm2 = z.real*z.real+z.imag*z.imag;
if(zm2>0.0)
{
ret.real = (this->real*z.real+this->imag*z.imag)/zm2;
ret.imag = (this->imag*z.real-this->real*z.imag)/zm2;
}
else
{
ret.real = (double) finf;
ret.imag = (double) finf;
}
return ret;
}
__host__ __device__ cucomp128 cucomp128::operator+(const double& z)
{
cucomp128 ret;
ret.real = this->real + z;
ret.imag = this->imag;
return ret;
}
__host__ __device__ cucomp128 cucomp128::operator-(const double& z)
{
cucomp128 ret;
ret.real = real-z;
ret.imag = imag;
return ret;
}
__host__ __device__ cucomp128 cucomp128::operator*(const double& z)
{
cucomp128 ret;
ret.real = real*z;
ret.imag = imag*z;
return ret;
}
__host__ __device__ cucomp128 cucomp128::operator/(const double& z)
{
cucomp128 ret;
if(z!=0.0f)
{
ret.real = real/z;
ret.imag = imag/z;
}
else
{
ret.real = finf;
ret.imag = finf;
}
return ret;
}
__host__ __device__ bool cucomp128::operator==(const cucomp128& z) const
{
bool ret = 0;
if(z.real == real && z.imag == imag)
{
ret = 1;
}
return ret;
}
__host__ __device__ bool cucomp128::operator!=(const cucomp128& z) const
{
return !(*this==z);
}
//sort first by real value, then by imaginary value
//this is done so that an ordering exists, as long as two values aren't equal
__host__ __device__ bool cucomp128::operator>(const cucomp128& z) const
{
bool ret = 0;
if(this->real>z.real)
{
ret = 1;
}
else if(this->real==z.real)
{
if(this->imag>z.imag)
{
ret = 1;
}
else
{
ret = 0;
}
}
else
{
ret = 0;
}
return ret;
}
__host__ __device__ bool cucomp128::operator<(const cucomp128& z) const
{
bool ret = 0;
if(this->real<z.real)
{
ret = 1;
}
else if(this->real==z.real)
{
if(this->imag<z.imag)
{
ret = 1;
}
else
{
ret = 0;
}
}
else
{
ret = 0;
}
return ret;
}
__host__ __device__ cucomp128 operator-(const cucomp128 &z)
{
cucomp128 ret;
ret.real = -z.real;
ret.imag = -z.imag;
return ret;
}
__host__ __device__ bool cucomp128::operator>=(const cucomp128& z) const
{
bool ret = (*this==z || *this>z);
return ret;
}
__host__ __device__ bool cucomp128::operator<=(const cucomp128& z) const
{
bool ret = (*this==z || *this<z);
return ret;
}
__host__ __device__ bool cucomp128::isnan() const
{
bool ret = 0;
if(::isnan(this->real) || ::isnan(this->imag))
{
ret = 1;
}
return ret;
}
__host__ __device__ bool cucomp128::isinf() const
{
bool ret = 0;
//calls math.h isinf()
if(::isinf(this->real) || ::isinf(this->imag))
{
ret = 1;
}
return ret;
}
__host__ __device__ bool cucomp128::isreal() const
{
bool ret = 1;
if(imag!=0.0f)
{
ret = 0;
}
return ret;
}
__host__ __device__ bool cucomp128::isimag() const
{
bool ret = 1;
if(real!=0.0f)
{
ret = 0;
}
return ret;
}
__host__ __device__ bool cucomp128::iszero() const
{
bool ret = 1;
if(real!=0.0f || imag!=0.0f)
{
ret = 0;
}
return ret;
}
__host__ __device__ double cucomp128::arg() const
{
double ret = 0.0;
ret = (double) amscuda::arg((double)real,(double)imag);
return ret;
}
__host__ __device__ double cucomp128::mag() const
{
double ret = 0.0;
ret = ::sqrt(real*real+imag*imag);
return ret;
}
__host__ __device__ cucomp128 cucomp128::conj() const
{
cucomp128 ret;
ret.real = real;
ret.imag = -imag;
return ret;
}
__host__ __device__ double arg(cucomp128 z)
{
return z.arg();
}
__host__ __device__ double abs(cucomp128 z)
{
return z.mag();
}
__host__ __device__ cucomp128 dtocomp(double _r, double _i)
{
cucomp128 ret;
ret.real = _r;
ret.imag = _i;
return ret;
}
__host__ __device__ double real(cucomp128 z)
{
return z.real;
}
__host__ __device__ double imag(cucomp128 z)
{
return z.imag;
}
__host__ __device__ cucomp128 sin(cucomp128 z)
{
cucomp128 ret;
cucomp128 im1 = dtocomp(0.0f,1.0f);
cucomp128 div = dtocomp(0.0f,2.0f);
ret = (exp(im1*z)-exp(-im1*z))/div;
return ret;
}
__host__ __device__ cucomp128 cos(cucomp128 z)
{
cucomp128 ret;
cucomp128 im1 = dtocomp(0.0f,1.0f);
cucomp128 div = dtocomp(2.0f,0.0f);
ret = (exp(im1*z)+exp(-im1*z))/div;
return ret;
}
__host__ __device__ cucomp128 tan(cucomp128 z)
{
return sin(z)/cos(z);
}
__host__ __device__ cucomp128 exp(cucomp128 z)
{
cucomp128 ret;
ret.real = ::exp(z.real)*::cos(z.imag);
ret.imag = ::exp(z.real)*::sin(z.imag);
return ret;
}
__host__ __device__ cucomp128 log(cucomp128 z)
{
cucomp128 ret;
ret.real = ::log(::sqrt(z.real*z.real+z.imag*z.imag));
ret.imag = amscuda::arg(z.real,z.imag);
return ret;
}
__host__ __device__ cucomp128 conj(cucomp128 z)
{
return z.conj();
}
__host__ __device__ cucomp128 cosh(cucomp128 z)
{
cucomp128 ret;
cucomp128 div = dtocomp(2.0f,0.0f);
ret = (exp(z)+exp(-z))/div;
return ret;
}
__host__ __device__ cucomp128 sinh(cucomp128 z)
{
cucomp128 ret;
cucomp128 div = dtocomp(2.0f,0.0f);
ret = (exp(z)-exp(-z))/div;
return ret;
}
__host__ __device__ cucomp128 tanh(cucomp128 z)
{
return sinh(z)/cosh(z);
}
__host__ __device__ cucomp128 pow(cucomp128 z1, cucomp128 z2)
{
cucomp128 ret;
if(z1.mag()>0.0)
ret = exp(z2*log(z1));
else
ret = dtocomp(0.0f,0.0f);
return ret;
}
void test_cucomp128_1()
{
cucomp128 z1;
cucomp128 a,b,c;
double d1;
double f1;
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)));
a = dtocomp(1.0,1.0);
b = dtocomp(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]);
c = a+b;
printf("c=a+b: c=%1.4f + %1.4fi\n",c[0],c[1]);
c = a-b;
printf("c=a-b: c=%1.4f + %1.4fi\n",c[0],c[1]);
c = a*b;
printf("c=a*b: c=%1.4f + %1.4fi\n",c[0],c[1]);
c = a/b;
printf("c=a/b: c=%1.4f + %1.4fi\n",c[0],c[1]);
f1 = abs(a);
printf("abs(a)=%1.4f\n",f1);
f1 = arg(a);
printf("abs(a)=%1.4f pi\n",f1/pi);
}
}; //end namespace cmp
}; //end namespace amscuda

476
src/amsculib2/amscu_comp64.cu Normal file
View File

@ -0,0 +1,476 @@
#include <amsculib2/amsculib2.hpp>
namespace amscuda
{
namespace cmp
{
__host__ __device__ cucomp64::cucomp64()
{
real = 0.0;
imag = 0.0;
return;
}
__host__ __device__ cucomp64::~cucomp64()
{
real = 0.0;
imag = 0.0;
return;
}
__host__ __device__ cucomp64::cucomp64(const cucomp64 &other)
{
real = other.real;
imag = other.imag;
return;
}
__host__ __device__ cucomp64::cucomp64(const float &other)
{
real = other;
imag = 0.0;
return;
}
__host__ __device__ cucomp64& cucomp64::operator=(cucomp64& other)
{
real = other.real;
imag = other.imag;
return *this;
}
__host__ __device__ const cucomp64& cucomp64::operator=(const cucomp64& other)
{
this->real = other.real;
this->imag = other.imag;
return *this;
}
__host__ __device__ cucomp64& cucomp64::operator=(float& other)
{
real = other;
imag = 0.0;
return *this;
}
__host__ __device__ const cucomp64& cucomp64::operator=(const float& other)
{
this->real = other;
this->imag = 0.0;
return *this;
}
__host__ __device__ float& cucomp64::operator[](int& ind)
{
if(ind==0)
{
return this->real;
}
else
{
return this->imag;
}
}
__host__ __device__ const float& cucomp64::operator[](const int& ind) const
{
if(ind==0)
{
return this->real;
}
else
{
return this->imag;
}
}
__host__ __device__ cucomp64 cucomp64::operator+(const cucomp64& z)
{
cucomp64 ret;
ret.real = real + z.real;
ret.imag = imag + z.imag;
return ret;
}
__host__ __device__ cucomp64 cucomp64::operator-(const cucomp64& z)
{
cucomp64 ret;
ret.real = real - z.real;
ret.imag = imag - z.imag;
return ret;
}
__host__ __device__ cucomp64 cucomp64::operator*(const cucomp64& z)
{
cucomp64 ret;
ret.real = (real*z.real - imag*z.imag);
ret.imag = (imag*z.real + real*z.imag);
return ret;
}
__host__ __device__ cucomp64 cucomp64::operator/(const cucomp64& z)
{
cucomp64 ret;
float zm2 = z.real*z.real+z.imag*z.imag;
if(zm2>0.0)
{
ret.real = (this->real*z.real+this->imag*z.imag)/zm2;
ret.imag = (this->imag*z.real-this->real*z.imag)/zm2;
}
else
{
ret.real = (float) finf;
ret.imag = (float) finf;
}
return ret;
}
__host__ __device__ cucomp64 cucomp64::operator+(const float& z)
{
cucomp64 ret;
ret.real = this->real + z;
ret.imag = this->imag;
return ret;
}
__host__ __device__ cucomp64 cucomp64::operator-(const float& z)
{
cucomp64 ret;
ret.real = real-z;
ret.imag = imag;
return ret;
}
__host__ __device__ cucomp64 cucomp64::operator*(const float& z)
{
cucomp64 ret;
ret.real = real*z;
ret.imag = imag*z;
return ret;
}
__host__ __device__ cucomp64 cucomp64::operator/(const float& z)
{
cucomp64 ret;
if(z!=0.0f)
{
ret.real = real/z;
ret.imag = imag/z;
}
else
{
ret.real = finf;
ret.imag = finf;
}
return ret;
}
__host__ __device__ bool cucomp64::operator==(const cucomp64& z) const
{
bool ret = 0;
if(z.real == real && z.imag == imag)
{
ret = 1;
}
return ret;
}
__host__ __device__ bool cucomp64::operator!=(const cucomp64& z) const
{
return !(*this==z);
}
//sort first by real value, then by imaginary value
//this is done so that an ordering exists, as long as two values aren't equal
__host__ __device__ bool cucomp64::operator>(const cucomp64& z) const
{
bool ret = 0;
if(this->real>z.real)
{
ret = 1;
}
else if(this->real==z.real)
{
if(this->imag>z.imag)
{
ret = 1;
}
else
{
ret = 0;
}
}
else
{
ret = 0;
}
return ret;
}
__host__ __device__ bool cucomp64::operator<(const cucomp64& z) const
{
bool ret = 0;
if(this->real<z.real)
{
ret = 1;
}
else if(this->real==z.real)
{
if(this->imag<z.imag)
{
ret = 1;
}
else
{
ret = 0;
}
}
else
{
ret = 0;
}
return ret;
}
__host__ __device__ cucomp64 operator-(const cucomp64 &z)
{
cucomp64 ret;
ret.real = -z.real;
ret.imag = -z.imag;
return ret;
}
__host__ __device__ bool cucomp64::operator>=(const cucomp64& z) const
{
bool ret = (*this==z || *this>z);
return ret;
}
__host__ __device__ bool cucomp64::operator<=(const cucomp64& z) const
{
bool ret = (*this==z || *this<z);
return ret;
}
__host__ __device__ bool cucomp64::isnan() const
{
bool ret = 0;
if(::isnan(this->real) || ::isnan(this->imag))
{
ret = 1;
}
return ret;
}
__host__ __device__ bool cucomp64::isinf() const
{
bool ret = 0;
//calls math.h isinf()
if(::isinf(this->real) || ::isinf(this->imag))
{
ret = 1;
}
return ret;
}
__host__ __device__ bool cucomp64::isreal() const
{
bool ret = 1;
if(imag!=0.0f)
{
ret = 0;
}
return ret;
}
__host__ __device__ bool cucomp64::isimag() const
{
bool ret = 1;
if(real!=0.0f)
{
ret = 0;
}
return ret;
}
__host__ __device__ bool cucomp64::iszero() const
{
bool ret = 1;
if(real!=0.0f || imag!=0.0f)
{
ret = 0;
}
return ret;
}
__host__ __device__ float cucomp64::arg() const
{
float ret = 0.0;
ret = (float) amscuda::arg((double)real,(double)imag);
return ret;
}
__host__ __device__ float cucomp64::mag() const
{
float ret = 0.0;
ret = ::sqrt(real*real+imag*imag);
return ret;
}
__host__ __device__ cucomp64 cucomp64::conj() const
{
cucomp64 ret;
ret.real = real;
ret.imag = -imag;
return ret;
}
__host__ __device__ float arg(cucomp64 z)
{
return z.arg();
}
__host__ __device__ float abs(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)
{
return z.real;
}
__host__ __device__ float imag(cucomp64 z)
{
return z.imag;
}
__host__ __device__ cucomp64 sin(cucomp64 z)
{
cucomp64 ret;
cucomp64 im1 = dtocomp64(0.0f,1.0f);
cucomp64 div = dtocomp64(0.0f,2.0f);
ret = (exp(im1*z)-exp(-im1*z))/div;
return ret;
}
__host__ __device__ cucomp64 cos(cucomp64 z)
{
cucomp64 ret;
cucomp64 im1 = dtocomp64(0.0f,1.0f);
cucomp64 div = dtocomp64(2.0f,0.0f);
ret = (exp(im1*z)+exp(-im1*z))/div;
return ret;
}
__host__ __device__ cucomp64 tan(cucomp64 z)
{
return sin(z)/cos(z);
}
__host__ __device__ cucomp64 exp(cucomp64 z)
{
cucomp64 ret;
ret.real = ::exp(z.real)*::cos(z.imag);
ret.imag = ::exp(z.real)*::sin(z.imag);
return ret;
}
__host__ __device__ cucomp64 log(cucomp64 z)
{
cucomp64 ret;
ret.real = ::log(::sqrt(z.real*z.real+z.imag*z.imag));
ret.imag = amscuda::arg(z.real,z.imag);
return ret;
}
__host__ __device__ cucomp64 conj(cucomp64 z)
{
return z.conj();
}
__host__ __device__ cucomp64 cosh(cucomp64 z)
{
cucomp64 ret;
cucomp64 div = dtocomp64(2.0f,0.0f);
ret = (exp(z)+exp(-z))/div;
return ret;
}
__host__ __device__ cucomp64 sinh(cucomp64 z)
{
cucomp64 ret;
cucomp64 div = dtocomp64(2.0f,0.0f);
ret = (exp(z)-exp(-z))/div;
return ret;
}
__host__ __device__ cucomp64 tanh(cucomp64 z)
{
return sinh(z)/cosh(z);
}
__host__ __device__ cucomp64 pow(cucomp64 z1, cucomp64 z2)
{
cucomp64 ret;
if(z1.mag()>0.0)
ret = exp(z2*log(z1));
else
ret = dtocomp64(0.0f,0.0f);
return ret;
}
void test_cucomp64_1()
{
cucomp64 z1;
cucomp64 a,b,c;
double d1;
float f1;
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)));
a = dtocomp64(1.0,1.0);
b = dtocomp64(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]);
c = a+b;
printf("c=a+b: c=%1.4f + %1.4fi\n",c[0],c[1]);
c = a-b;
printf("c=a-b: c=%1.4f + %1.4fi\n",c[0],c[1]);
c = a*b;
printf("c=a*b: c=%1.4f + %1.4fi\n",c[0],c[1]);
c = a/b;
printf("c=a/b: c=%1.4f + %1.4fi\n",c[0],c[1]);
f1 = abs(a);
printf("abs(a)=%1.4f\n",f1);
f1 = arg(a);
printf("abs(a)=%1.4f pi\n",f1/pi);
}
}; //end namespace cmp
}; //end namespace amscuda

21
src/amsculib2/amscu_cudafunctions.cu Normal file
View File

@ -0,0 +1,21 @@
#include <amsculib2/amsculib2.hpp>
namespace amscuda
{
int cuda_errortrap(const char *msgheader)
{
int ret = 0;
cudaError_t err = cudaSuccess;
err = cudaGetLastError();
if(err!=cudaSuccess)
{
printf("%s :%s\n",msgheader,cudaGetErrorString(err));
ret = 1;
}
return ret;
}
}; //end namespace amscuda

222
src/amsculib2/amscu_random.cu Normal file
View File

@ -0,0 +1,222 @@
#include <amsculib2/amsculib2.hpp>
namespace amscuda
{
__device__ __host__ float fhash1d_su(float x)
{
float ret;
ret = x*(x>0.0f) + -x*(x<0.0f); //sign without conditionals?
ret = fmodf(ret,10000.0f); //restrain domain
ret = fmodf(ret*(ret+3678.453f)+7890.453f,10000.0f);
ret = fmodf(ret*(ret+8927.2134f),10000.0f);
ret = fmodf(ret*(ret+3656.234f),10000.0f);
//ret = fmodf(ret*(ret+892.2134f),1000.0f);
//ret = fmodf(ret,1000.0f);
ret = ret/10000.0f;
return ret;
}
__device__ __host__ float fhash3d_su(float x, float y=0.0f, float z=0.0f)
{
float ret = 0.0f;
ret = fhash1d_su(z);
ret = fhash1d_su(1000.0f*ret*ret + y);
ret = fhash1d_su(1000.0f*ret*ret + x);
return ret;
}
__device__ __host__ float fhash4d_su(float x, float y=0.0f, float z=0.0f, float w=0.0f)
{
float ret = 0.0f;
ret = fhash1d_su(w);
ret = fhash1d_su(1000.0f*ret*ret + z);
ret = fhash1d_su(1000.0f*ret*ret + y);
ret = fhash1d_su(1000.0f*ret*ret + x);
return ret;
}
//////////////////////////////////////////////////
// Deterministic Pseudorandom int32_t Generator //
//////////////////////////////////////////////////
//Simple 32 bit integer deterministic pseudo-random generator
// *not* for cryptography
// Frequency of generated floats should be uniform [0,1)
AMSCU_CONST static const int32_t dpr32_mod = 1<<30-1;
AMSCU_CONST static const int32_t dpr32_mult = 25137;
//Next seed in simple 32 bit integer deterministic psuedo-rand generator
__host__ __device__ void dpr32_nextseed(int32_t *rseed_inout)
{
int32_t lseed;
if(rseed_inout!=NULL) lseed = *rseed_inout;
lseed = (lseed*dpr32_mult + 1)%dpr32_mod;
lseed = (lseed>=0)*(lseed)+(lseed<0)*(lseed+dpr32_mod); //ensure mod is positive
if(rseed_inout!=NULL) *rseed_inout = lseed;
return;
}
//Simple 32 bit integer deterministic pseudo-random generator
// *not* for cryptography
// Frequency of generated floats should be uniform [0,1)
__host__ __device__ float dpr32_randf(int32_t *rseed_inout)
{
int32_t lseed = 1;
float ret = 0.0f;
if(rseed_inout!=NULL) lseed = *rseed_inout;
dpr32_nextseed(&lseed);
ret = ((float)(lseed))/((float)dpr32_mod);
if(rseed_inout!=NULL) *rseed_inout = lseed;
return ret;
}
//box muller standard normal variable
__host__ __device__ float dpr32_randnf(int32_t *rseed_inout)
{
int32_t lseed = 1;
float ret = 0.0f;
float u1,u2;
if(rseed_inout!=NULL) lseed = *rseed_inout;
u1 = dpr32_randf(&lseed);
u2 = dpr32_randf(&lseed);
ret = ::sqrtf(-2.0f*::logf(u1))*::cosf(2.0f*pif*u2);
if(rseed_inout!=NULL) *rseed_inout = lseed;
return ret;
}
//////////////////////////////////////////////////
// Deterministic Pseudorandom int64_t Generator //
//////////////////////////////////////////////////
//"goodenough" deterministic pseudo-random number generator
//random enough for procedural applications, deterministic,
//operates without side-effects for thread safety
AMSCU_CONST const int64_t random_dpr64_mod = (2LL<<31LL)-1LL;
AMSCU_CONST const int64_t random_dpr64_mult = 1201633LL;
__host__ __device__ void dpr64_nextseed(int64_t *seedinout)
{
int64_t lseed = 0LL;
if(seedinout!=NULL) lseed = *seedinout;
lseed = (random_dpr64_mult*lseed+1LL)%random_dpr64_mod;
lseed = (lseed>=0)*(lseed)+(lseed<0)*(lseed+random_dpr64_mod);
if(seedinout!=NULL) *seedinout = lseed;
return;
}
__host__ __device__ double dpr64_randd(int64_t *seedinout)
{
double ret = 0.0;
int64_t lseed = 0LL;
if(seedinout!=NULL) lseed = *seedinout;
dpr64_nextseed(&lseed);
ret = ((double)lseed)/((double)(random_dpr64_mod-1LL));
if(seedinout!=NULL) *seedinout = lseed;
return ret;
}
__host__ __device__ float dpr64_randf(int64_t *seedinout)
{
float ret = 0.0f;
int64_t lseed = 0LL;
if(seedinout!=NULL) lseed = *seedinout;
dpr64_nextseed(&lseed);
ret = ((float)lseed)/((float)(random_dpr64_mod-1LL));
if(seedinout!=NULL) *seedinout = lseed;
return ret;
}
///////////
// Tests //
///////////
void test_dprg64()
{
printf("Tests for dprg:\n");
long I;
int64_t seed = 133LL;
double d;
float f;
cuvect3 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("\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;
}
void test_dprg32()
{
printf("Tests for dprg:\n");
long I;
int32_t seed = 133;
double d;
float f;
cuvect3 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("\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;
}
}; //namespace amscuda

63
src/amsculib2/amscuarray.cu Normal file
View File

@ -0,0 +1,63 @@
#include <amsculib2/amsculib2.hpp>
namespace amscuda
{
__global__ void test_cuarray_sum_kf(cuarray<float> *dq1, float *sum)
{
int I;
*sum = 0.0f;
for(I=0;I<dq1->length;I++)
{
*sum = *sum + dq1->data[I];
}
//*sum = (float)dq1->length;
return;
}
float test_cuarray_sum(cuarray<float> *q1)
{
float ret = 0.0f;
int res;
cuarray<float> *dq1 = NULL;
float *dsum;
cudaError_t err = cudaSuccess;
cudaMalloc(&dsum,sizeof(float));
res = q1->device_send(&dq1);
printf("error: res=%d\n",res);
test_cuarray_sum_kf<<<1,1>>>(dq1,dsum);
cudaDeviceSynchronize();
err = cudaGetLastError();
if(err!=cudaSuccess)
{
printf("test_cuarray_sum Error: %s\n",cudaGetErrorString(err));
}
cudaMemcpy(&ret,dsum,sizeof(float),cudaMemcpyDeviceToHost);
q1->device_free(&dq1);
cudaFree(dsum); dsum = NULL;
return ret;
}
void test_cuarray()
{
cuarray<float> q1;
int I;
q1.resize(100);
for(I=0;I<q1.length;I++)
{
q1.data[I] = I;
}
printf("q1.length=%d\n",q1.length);
printf("sum of array: %1.6g\n",test_cuarray_sum(&q1));
}
};

213
src/amsculib2/amscuarray_dops.cu Normal file
View File

@ -0,0 +1,213 @@
#include <amsculib2/amsculib2.hpp>
namespace amscuda
{
//template instantiations
template float dbuff_sum(float *devbuffer, int N);
template void dbuff_minmax(float *devbuffer, int N, float *min, float *max);
template void dbuff_setall(float *devbuffer, int N, float setto, int nblocks, int nthreads);
//fill devbuffer with random uniform numbers between 0 and 1 using int32_t based generator
__global__ void dbuff_rand_dpr32_kf(float *devbuffer, int N, int32_t *seeds)
{
int I0 = threadIdx.x + blockIdx.x*blockDim.x;
int Is = blockDim.x*gridDim.x;
int I;
int32_t lseed;
float f;
lseed = seeds[I0];
for(I=I0;I<N;I=I+Is)
{
f = dpr32_randf(&lseed);
devbuffer[I] = f;
}
return;
}
void dbuff_rand_dpr32(float *devbuffer, int N, int32_t *rseedinout, int nblocks, int nthreads)
{
cudaError_t err = cudaSuccess;
int I;
int32_t *seeds = NULL;
int32_t *devseeds = NULL;
int32_t lseed;
if(devbuffer==NULL || N<=0)
{
return;
}
seeds = new(std::nothrow) int32_t[nblocks*nthreads];
cudaMalloc(&devseeds,sizeof(int32_t)*nblocks*nthreads);
if(rseedinout!=NULL) lseed = *rseedinout; else lseed = 1;
for(I=0;I<nblocks*nthreads;I++)
{
lseed = lseed + I + 1;
dpr32_nextseed(&lseed);
seeds[I] = lseed;
}
cudaMemcpy(devseeds,seeds,sizeof(int32_t)*nblocks*nthreads,cudaMemcpyHostToDevice);
if(rseedinout!=NULL) *rseedinout = lseed;
dbuff_rand_dpr32_kf<<<nblocks,nthreads>>>(devbuffer,N,devseeds);
cudaDeviceSynchronize();
err = cudaGetLastError();
if(err!=cudaSuccess)
{
printf("amscu::dbuff_rand_dpr32 error: %s\n",cudaGetErrorString(err));
}
cudaFree(devseeds); devseeds = NULL;
delete[] seeds; seeds = NULL;
return;
}
__global__ void dbuff_rand_dpr32n_kf(float *devbuffer, int N, int32_t *seeds)
{
int I0 = threadIdx.x + blockIdx.x*blockDim.x;
int Is = blockDim.x*gridDim.x;
int I;
int32_t lseed;
float f;
lseed = seeds[I0];
for(I=I0;I<N;I=I+Is)
{
f = dpr32_randnf(&lseed);
devbuffer[I] = f;
}
return;
}
void dbuff_rand_dpr32n(float *devbuffer, int N, int32_t *rseedinout, int nblocks, int nthreads)
{
cudaError_t err = cudaSuccess;
int I;
int32_t *seeds = NULL;
int32_t *devseeds = NULL;
int32_t lseed;
if(devbuffer==NULL || N<=0)
{
return;
}
seeds = new(std::nothrow) int32_t[nblocks*nthreads];
cudaMalloc(&devseeds,sizeof(int32_t)*nblocks*nthreads);
if(rseedinout!=NULL) lseed = *rseedinout; else lseed = 1;
for(I=0;I<nblocks*nthreads;I++)
{
lseed = lseed + I + 1;
dpr32_nextseed(&lseed);
seeds[I] = lseed;
}
cudaMemcpy(devseeds,seeds,sizeof(int32_t)*nblocks*nthreads,cudaMemcpyHostToDevice);
if(rseedinout!=NULL) *rseedinout = lseed;
dbuff_rand_dpr32n_kf<<<nblocks,nthreads>>>(devbuffer,N,devseeds);
cudaDeviceSynchronize();
err = cudaGetLastError();
if(err!=cudaSuccess)
{
printf("amscu::dbuff_rand_dpr32 error: %s\n",cudaGetErrorString(err));
}
cudaFree(devseeds); devseeds = NULL;
delete[] seeds; seeds = NULL;
return;
}
void dbuff_rand_dpr64(float *devbuffer, int N, int64_t *rseedinout, int nblocks, int nthreads)
{
int I;
return;
}
///////////
// Tests //
///////////
void test_dbuff_rand_dpr32()
{
cuarray<float> data;
float *dev_data = NULL;
int Nx = 5000;
int Ny = 5000;
cuarray<int> dims;
int32_t rseed = 15;
FILE *fp = NULL;
const char *fname = "./test_scripts/test_dbuff_rand_dpr32.bin";
clock_t t0,t1,t2;
double dt0,dt1;
printf("Tests of dbuff_rand_dpr32...\n");
fp = fopen(fname,"w+");
if(fp==NULL)
{
printf("Error: Could not open %s for writing.\n",fname);
return;
}
data.resize(Nx*Ny);
dims.resize(2);
dims[0] = Nx; dims[1] = Ny;
cudaMalloc(&dev_data,Nx*Ny*sizeof(float));
t0 = clock();
dbuff_rand_dpr32(dev_data,Nx*Ny,&rseed,256,512);
t1 = clock();
cudaMemcpy(data.data,dev_data,Nx*Ny*sizeof(float),cudaMemcpyDeviceToHost);
t2 = clock();
dt0 = (double)(t1-t0)/(double)CLOCKS_PER_SEC*1000.0;
dt1 = (double)(t2-t1)/(double)CLOCKS_PER_SEC*1000.0;
printf("dbuff_rand_dpr32 exec time: %1.3f msec\n",dt0);
printf("copy (%d,%d) to host time: %1.3f msec\n",Nx,Ny,dt1);
fwrite_ndarray(fp,&dims,&data);
t0 = clock();
dbuff_rand_dpr32n(dev_data,Nx*Ny,&rseed,256,512);
t1 = clock();
cudaMemcpy(data.data,dev_data,Nx*Ny*sizeof(float),cudaMemcpyDeviceToHost);
t2 = clock();
dt0 = (double)(t1-t0)/(double)CLOCKS_PER_SEC*1000.0;
dt1 = (double)(t2-t1)/(double)CLOCKS_PER_SEC*1000.0;
printf("dbuff_rand_dpr32n exec time: %1.3f msec\n",dt0);
printf("copy (%d,%d) to host time: %1.3f msec\n",Nx,Ny,dt1);
fwrite_ndarray(fp,&dims,&data);
fclose(fp);
cudaFree(dev_data); dev_data = NULL;
}
};

6
src/amsculib2/amscugeom.cu Normal file
View File

@ -0,0 +1,6 @@
#include <amsculib2/amsculib2.hpp>
namespace amscuda
{
};

6
src/amsculib2/amsculib2.cu Normal file
View File

@ -0,0 +1,6 @@
#include <amsculib2/amsculib2.hpp>
namespace amscuda
{
};

269
src/amsculib2/amscumath.cu Normal file
View File

@ -0,0 +1,269 @@
#include <amsculib2/amsculib2.hpp>
namespace amscuda
{
__host__ __device__ double dabs(double x)
{
if(x<0.0)
{
x = -x;
}
return x;
}
__host__ __device__ float fabs(float x)
{
if(x<0.0f)
{
x = -x;
}
return x;
}
__host__ __device__ double mod(double x, double md)
{
x = fmod(x,md);
if(x<0.0)
{
x = x + md;
}
return x;
}
__host__ __device__ float mod(float x, float md)
{
x = fmodf(x,md);
if(x<0.0f)
{
x = x + md;
}
return x;
}
__host__ __device__ int mod(int x, int n)
{
x = x % n;
if(x<0)
{
x = x + n;
}
return x;
}
__host__ __device__ long mod(long x, long n)
{
x = x % n;
if(x<0)
{
x = x + n;
}
return x;
}
__host__ __device__ int truediv(int x, int y)
{
int z = 0;
if(x>=0 && y>0)
{
z = x/y;
}
else if(x<0 && y>0)
{
z = -((-x)/y) - 1;
}
else if(x>=0 && y<0)
{
z = -(x/(-y)) - 1;
}
else if(x<0 && y<0)
{
z = ((-x)/(-y));
}
return z;
}
__host__ __device__ long truediv(long x, long y)
{
int z = 0;
if(x>=0 && y>0)
{
z = x/y;
}
else if(x<0 && y>0)
{
z = -((-x)/y) - 1;
}
else if(x>=0 && y<0)
{
z = -(x/(-y)) - 1;
}
else if(x<0 && y<0)
{
z = ((-x)/(-y));
}
return z;
}
template<> __host__ __device__ double min(double a, double b)
{
if(isnan(a))
{
return b;
}
else if(isnan(b))
{
return a;
}
else if(a>b)
{
return b;
}
else
{
return a;
}
}
template<> __host__ __device__ float min(float a, float b)
{
if(isnan(a))
{
return b;
}
else if(isnan(b))
{
return a;
}
else if(a>b)
{
return b;
}
else
{
return a;
}
}
template<> __host__ __device__ double max(double a, double b)
{
if(isnan(a))
{
return b;
}
else if(isnan(b))
{
return a;
}
else if(a>b)
{
return a;
}
else
{
return b;
}
}
template<> __host__ __device__ float max(float a, float b)
{
if(isnan(a))
{
return b;
}
else if(isnan(b))
{
return a;
}
else if(a>b)
{
return a;
}
else
{
return b;
}
}
__device__ __host__ double arg(double x, double y)
{
double ret = 0.0;
double z = ::sqrt(x*x+y*y);
if(z>0.0)
{
if(y<=x && y>=-x)
{
ret = asin(y/z);
}
else if(y>=x && y>=-x)
{
ret = acos(x/z);
}
else if(y>=x && y<=-x)
{
ret = pi-asin(y/z);
}
else
{
ret = 2.0*pi-acos(x/z);
}
}
if(ret<0.0)
{
ret = 2.0*pi+ret;
}
return ret;
}
__device__ __host__ void get_azel(double x, double y, double z, double *az, double *el)
{
//int ret = -2; //should never see this return
double n, rp;
n = ::sqrt(x*x+y*y+z*z);
if(n>0.0)
{
rp = ::sqrt(x*x+y*y);
if(rp>0.0)
{
//ret = 1; //nonzero vector - should work
*az = arg(x,y);
*el = ::atan(z/rp);
}
else
{
//ret = 0; //straight up or straight down
if(z>0.0)
{
*az = 0.0;
*el = pi/2.0;
}
else
{
*az = 0.0;
*el = -pi/2.0;
}
}
}
else
{
*az = 0.0;
*el = 0.0;
//ret = -1; //zero vector - no real azimuth/elevation
}
return;
}
void test_amscumath1()
{
printf("pi = %1.16f\n",amscuda::pi);
}
};

85
src/amsculib2/amscurarray.cu Normal file
View File

@ -0,0 +1,85 @@
#include <amsculib2/amsculib2.hpp>
namespace amscuda
{
__global__ void test_amscurarray1_kf1(curarray<int> *q)
{
int I,J,Na;
int I0 = threadIdx.x + blockIdx.x*blockDim.x;
int Is = blockDim.x*gridDim.x;
int N = q->Narrays;
for(I=I0; I<N;I = I + Is)
{
Na = q->N[I];
//printf("I:%d Na: %d\n",I,Na);
//q->dev_resizearray(I, Na);
for(J=0;J<Na;J++)
{
q->devarrayptrs[I][J] = J;
}
}
}
void test_amscurarray1()
{
int I;
cudaError_t err = cudaSuccess;
printf("test_amscurarray1:\n");
curarray<int> *qarray = NULL;
curarray_new(&qarray,100);
err = cudaGetLastError();
if(err!=cudaSuccess)
{
printf("debug error trap 1: %s\n",cudaGetErrorString(err));
}
for(I=0;I<100;I++)
{
qarray->resizearray(I,5);
}
qarray->push();
qarray->pull();
cuda_errortrap("debug: error trap 2");
for(I=0;I<5;I++)
{
printf("array[%d], size %d\n",I,qarray->N[I]);
}
//
for(I=0;I<100;I++)
{
qarray->resizearray(I,I%5);
cuda_errortrap("debug: error trap resize2");
}
qarray->push();
qarray->pull();
test_amscurarray1_kf1<<<128,1>>>(qarray->devptr);
cuda_errortrap("debug: error trap kf1");
qarray->pull();
cuda_errortrap("debug: error trap pull2");
for(I=0;I<5;I++)
{
printf("array[%d], size %d\n",I,qarray->N[I]);
}
curarray_delete(&qarray);
cuda_errortrap("debug: error trap 3");
return;
}
};

361
src/amsculib2/cuvect2.cu Normal file
View File

@ -0,0 +1,361 @@
#include <amsculib2/amsculib2.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;
}
void test_cuvect2_1()
{
return;
}
};

361
src/amsculib2/cuvect2f.cu Normal file
View File

@ -0,0 +1,361 @@
#include <amsculib2/amsculib2.hpp>
namespace amscuda
{
__host__ __device__ cuvect2f::cuvect2f()
{
x = 0.0; y = 0.0;
return;
}
__host__ __device__ cuvect2f::~cuvect2f()
{
x = 0.0; y = 0.0;
return;
}
__host__ __device__ cuvect2f::cuvect2f(float _x, float _y)
{
x = _x; y = _y;
return;
}
__host__ __device__ float& cuvect2f::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
{
if(I==0) return x;
if(I==1) return y;
return x;
}
__host__ __device__ cuvect2f cuvect2f::operator+(cuvect2f lhs)
{
cuvect2f ret;
ret.x = x + lhs.x;
ret.y = y + lhs.y;
return ret;
}
__host__ __device__ cuvect2f cuvect2f::operator-(cuvect2f lhs)
{
cuvect2f ret;
ret.x = x - lhs.x;
ret.y = y - lhs.y;
return ret;
}
__host__ __device__ cuvect2f cuvect2f::operator*(float lhs)
{
cuvect2f ret;
ret.x = x*lhs;
ret.y = y*lhs;
return ret;
}
__host__ __device__ cuvect2f cuvect2f::operator/(float lhs)
{
cuvect2f ret;
ret.x = x/lhs;
ret.y = y/lhs;
return ret;
}
__host__ __device__ float cuvect2f_dot(cuvect2f a, cuvect2f b)
{
float ret = a.x*b.x+a.y*b.y;
return ret;
}
__host__ __device__ float cuvect2f_cross(cuvect2f a, cuvect2f b)
{
float ret = a.x*b.y-a.y*b.x;
return ret;
}
__host__ __device__ float cuvect2f_norm(cuvect2f a)
{
float ret = ::sqrtf(a.x*a.x+a.y*a.y);
return ret;
}
__host__ __device__ cuvect2f cuvect2f_normalize(cuvect2f a)
{
cuvect2f ret;
float m = cuvect2f_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__ cuvect2f cuvect2f_proj(cuvect2f a, cuvect2f b)
{
cuvect2f ret;
cuvect2f bn = cuvect2f_normalize(b);
float m = cuvect2f_dot(a,bn);
ret = bn*m;
return ret;
}
__host__ __device__ cumat2f::cumat2f()
{
int I;
for(I=0;I<4;I++)
{
dat[I] = 0.0;
}
return;
}
__host__ __device__ cumat2f::~cumat2f()
{
int I;
for(I=0;I<4;I++)
{
dat[I] = 0.0;
}
return;
}
__host__ __device__ float& cumat2f::operator[](const int I)
{
return dat[I];
}
__host__ __device__ float& cumat2f::operator()(const int I, const int J)
{
return dat[I+2*J];
}
__host__ __device__ cumat2f cumat2f::operator+(cumat2f lhs)
{
int I;
cumat2f ret;
for(I=0;I<4;I++)
{
ret.dat[I] = this->dat[I] + lhs.dat[I];
}
return ret;
}
__host__ __device__ cumat2f cumat2f::operator-(cumat2f lhs)
{
int I;
cumat2f ret;
for(I=0;I<4;I++)
{
ret.dat[I] = this->dat[I] - lhs.dat[I];
}
return ret;
}
__host__ __device__ cumat2f cumat2f::operator*(float lhs)
{
cumat2f ret;
int I;
for(I=0;I<4;I++)
{
ret.dat[I] = this->dat[I]*lhs;
}
return ret;
}
__host__ __device__ cumat2f cumat2f::operator/(float lhs)
{
cumat2f ret;
int I;
for(I=0;I<4;I++)
{
ret.dat[I] = this->dat[I]/lhs;
}
return ret;
}
__host__ __device__ float& cumat2f::at(const int I, const int J)
{
return dat[I+2*J];
}
__host__ __device__ cuvect2f cumat2f::operator*(cuvect2f lhs)
{
cuvect2f ret = cuvect2f(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__ cumat2f cumat2f::operator*(cumat2f lhs)
{
cumat2f 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__ float cumat2f::det()
{
float 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__ cumat2f cumat2f::transpose()
{
cumat2f 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__ cumat2f cumat2f::inverse()
{
cumat2f q;
float 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 mat2f_transpose(float *mat2inout)
{
float mat2f_in[4];
mat2f_copy(mat2f_in,mat2inout);
mat2inout[0+0*2] = mat2f_in[0+0*2];
mat2inout[1+0*2] = mat2f_in[0+1*2];
mat2inout[0+1*2] = mat2f_in[1+0*2];
mat2inout[1+1*2] = mat2f_in[1+1*2];
return;
}
//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)
{
return;
}
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, 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;
}
void test_cuvect2f_1()
{
return;
}
};

581
src/amsculib2/cuvect3.cu Normal file
View File

@ -0,0 +1,581 @@
#include <amsculib2/amsculib2.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+(cuvect3 lhs)
{
cuvect3 ret;
ret.x = x+lhs.x;
ret.y = y+lhs.y;
ret.z = z+lhs.z;
return ret;
}
__host__ __device__ cuvect3 cuvect3::operator-(cuvect3 lhs)
{
cuvect3 ret;
ret.x = x-lhs.x;
ret.y = y-lhs.y;
ret.z = z-lhs.z;
return ret;
}
__host__ __device__ cuvect3 cuvect3::operator*(double lhs)
{
cuvect3 ret;
ret.x = x*lhs;
ret.y = y*lhs;
ret.z = z*lhs;
return ret;
}
__host__ __device__ cuvect3 cuvect3::operator/(double lhs)
{
cuvect3 ret;
ret.x = x/lhs;
ret.y = y/lhs;
ret.z = z/lhs;
return ret;
}
__host__ __device__ cuvect3::cuvect3(double _x, double _y, double _z)
{
x = _x; y = _y; z = _z;
return;
}
__host__ __device__ double cuvect3_dot(cuvect3 a, cuvect3 b)
{
double ret = a.x*b.x+a.y*b.y+a.z*b.z;
return ret;
}
__host__ __device__ cuvect3 cuvect3_cross(cuvect3 a, 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(cuvect3 a)
{
double ret;
ret = ::sqrt(a.x*a.x+a.y*a.y+a.z*a.z);
return ret;
}
__host__ __device__ cuvect3 cuvect3_normalize(cuvect3 a)
{
cuvect3 ret;
double m;
m = ::sqrt(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(cuvect3 a, cuvect3 b)
{
cuvect3 ret;
cuvect3 bn = cuvect3_normalize(b);
double m = cuvect3_dot(a,bn);
ret = bn*m;
return ret;
}
//transposes a 3x3 (9 element) matrix
__host__ __device__ void mat3_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 mat3_copy(double *mat3_dest, const double *mat3_src)
{
int I;
if(mat3_dest==NULL || mat3_src==NULL)
return;
for(I=0;I<9;I++)
mat3_dest[I] = mat3_src[I];
return;
}
__host__ __device__ double mat3_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 mat3_inverse(double *mat3inout)
{
int I;
double matint[9];
double det = mat3_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;
mat3_transpose(mat3inout);
return;
}
__host__ __device__ cuvect3 mat3_mult(double *mat3in, 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 mat3_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__ __device__ void mat3_hodgedual(cuvect3 vecin, double *matout)
{
matout[0 + 0*3] = 0.0f;
matout[1 + 0*3] = -vecin[2];
matout[2 + 0*3] = vecin[1];
matout[0 + 1*3] = vecin[2];
matout[1 + 1*3] = 0.0f;
matout[2 + 1*3] = -vecin[0];
matout[0 + 2*3] = -vecin[1];
matout[1 + 2*3] = vecin[0];
matout[2 + 2*3] = 0.0f;
return;
}
__host__ __device__ void mat3_hodgedual(double *matin, cuvect3 vecout)
{
vecout[0] = 0.5*(matin[1 + 2*3] - matin[2 + 1*3]);
vecout[1] = 0.5*(matin[2 + 0*3] - matin[0 + 2*3]);
vecout[2] = 0.5*(matin[0 + 1*3] - matin[1 + 0*3]);
return;
}
//returns direction cosine rotation matrix from axis and angle
__host__ __device__ void mat3_rot_from_axisangle(cuvect3 axis, double angle, double *matout)
{
int I;
double H[9];
double Hsq[9];
double II[9];
for(I=0;I<9;I++) II[I] = 0.0;
II[0+0*3] = 1.0;
II[1+1*3] = 1.0;
II[2+2*3] = 1.0;
axis = cuvect3_normalize(axis);
mat3_hodgedual(axis,H);
mat3_mult(H,H,Hsq);
for(I=0;I<9;I++)
{
matout[I] = (II[I] + Hsq[I]) + H[I]*sin(angle) - Hsq[I]*cos(angle);
}
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];
}
}
mat3_inverse(mb);
mat3_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 = mat3_mult(ma,a);
b = mat3_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::cumat3()
{
int I;
for(I=0;I<9;I++)
{
dat[I] = 0.0;
}
return;
}
__host__ __device__ cumat3::~cumat3()
{
int I;
for(I=0;I<9;I++)
{
dat[I] = 0.0;
}
return;
}
__host__ __device__ double& cumat3::operator[](const int I)
{
return dat[I];
}
__host__ __device__ double& cumat3::operator()(const int I, const int J)
{
return dat[I+3*J];
}
__host__ __device__ cumat3 cumat3::operator+(cumat3 lhs)
{
int I;
cumat3 ret;
for(I=0;I<9;I++)
{
ret.dat[I] = this->dat[I] + lhs.dat[I];
}
return ret;
}
__host__ __device__ cumat3 cumat3::operator-(cumat3 lhs)
{
int I;
cumat3 ret;
for(I=0;I<9;I++)
{
ret.dat[I] = this->dat[I] - lhs.dat[I];
}
return ret;
}
__host__ __device__ cumat3 cumat3::operator*(double lhs)
{
cumat3 ret;
int I;
for(I=0;I<9;I++)
{
ret.dat[I] = this->dat[I]*lhs;
}
return ret;
}
__host__ __device__ cumat3 cumat3::operator/(double lhs)
{
cumat3 ret;
int I;
for(I=0;I<9;I++)
{
ret.dat[I] = this->dat[I]/lhs;
}
return ret;
}
__host__ __device__ double& cumat3::at(const int I, const int J)
{
return dat[I+3*J];
}
__host__ __device__ cuvect3 cumat3::operator*(cuvect3 lhs)
{
cuvect3 ret = cuvect3(0.0,0.0,0.0);
int I,J;
for(I=0;I<3;I++)
{
for(J=0;J<3;J++)
{
ret[I] = ret[I] + this->at(I,J)*lhs[J];
}
}
return ret;
}
__host__ __device__ cumat3 cumat3::operator*(cumat3 lhs)
{
cumat3 ret;
int I,J,K;
for(I=0;I<3;I++)
{
for(J=0;J<3;J++)
{
ret(I,J) = 0.0;
for(K=0;K<3;K++)
{
ret(I,J) = ret(I,J) + this->at(I,K)*lhs(K,J);
}
}
}
return ret;
}
__host__ __device__ double cumat3::det()
{
double ret = 0.0;
ret = ret + this->at(0,0)*this->at(1,1)*this->at(2,2);
ret = ret + this->at(0,1)*this->at(1,2)*this->at(2,0);
ret = ret + this->at(0,2)*this->at(1,0)*this->at(2,1);
ret = ret - this->at(0,0)*this->at(1,2)*this->at(2,1);
ret = ret - this->at(0,1)*this->at(1,0)*this->at(2,2);
ret = ret - this->at(0,2)*this->at(1,1)*this->at(2,0);
return ret;
}
__host__ __device__ cumat3 cumat3::transpose()
{
cumat3 q;
int I,J;
for(I=0;I<3;I++)
{
for(J=0;J<3;J++)
{
q.at(I,J) = this->at(J,I);
}
}
return q;
}
__host__ __device__ cumat3 cumat3::inverse()
{
cumat3 q;
double dt = q.det();
if(dt!=0.0)
{
q(0,0) = (this->at(1,1)*this->at(2,2)-this->at(1,2)*this->at(2,1))/dt;
q(0,1) = -(this->at(1,0)*this->at(2,2)-this->at(1,2)*this->at(2,0))/dt;
q(0,2) = (this->at(1,0)*this->at(2,1)-this->at(1,1)*this->at(2,0))/dt;
q(1,0) = -(this->at(0,1)*this->at(2,2)-this->at(0,2)*this->at(2,1))/dt;
q(1,1) = (this->at(0,0)*this->at(2,2)-this->at(0,2)*this->at(2,0))/dt;
q(1,2) = -(this->at(0,0)*this->at(2,1)-this->at(0,1)*this->at(2,0))/dt;
q(2,0) = (this->at(0,1)*this->at(1,2)-this->at(0,2)*this->at(1,1))/dt;
q(2,1) = -(this->at(0,0)*this->at(1,2)-this->at(0,2)*this->at(1,0))/dt;
q(2,2) = (this->at(0,0)*this->at(1,1)-this->at(0,1)*this->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;
}
};

581
src/amsculib2/cuvect3f.cu Normal file
View File

@ -0,0 +1,581 @@
#include <amsculib2/amsculib2.hpp>
namespace amscuda
{
__host__ __device__ cuvect3f::cuvect3f()
{
x = 0.0; y = 0.0; z = 0.0;
return;
}
__host__ __device__ cuvect3f::~cuvect3f()
{
x = 0.0; y = 0.0; z = 0.0;
return;
}
__host__ __device__ float& cuvect3f::operator[](const int I)
{
if(I==0) return x;
if(I==1) return y;
if(I==2) return z;
return x;
}
__host__ __device__ const float& cuvect3f::operator[](const int I) const
{
if(I==0) return x;
if(I==1) return y;
if(I==2) return z;
return x;
}
__host__ __device__ cuvect3f cuvect3f::operator+(cuvect3f lhs)
{
cuvect3f ret;
ret.x = x+lhs.x;
ret.y = y+lhs.y;
ret.z = z+lhs.z;
return ret;
}
__host__ __device__ cuvect3f cuvect3f::operator-(cuvect3f lhs)
{
cuvect3f ret;
ret.x = x-lhs.x;
ret.y = y-lhs.y;
ret.z = z-lhs.z;
return ret;
}
__host__ __device__ cuvect3f cuvect3f::operator*(float lhs)
{
cuvect3f ret;
ret.x = x*lhs;
ret.y = y*lhs;
ret.z = z*lhs;
return ret;
}
__host__ __device__ cuvect3f cuvect3f::operator/(float lhs)
{
cuvect3f ret;
ret.x = x/lhs;
ret.y = y/lhs;
ret.z = z/lhs;
return ret;
}
__host__ __device__ cuvect3f::cuvect3f(float _x, float _y, float _z)
{
x = _x; y = _y; z = _z;
return;
}
__host__ __device__ float cuvect3f_dot(cuvect3f a, cuvect3f b)
{
float ret = a.x*b.x+a.y*b.y+a.z*b.z;
return ret;
}
__host__ __device__ cuvect3f cuvect3f_cross(cuvect3f a, cuvect3f b)
{
cuvect3f 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__ float cuvect3f_norm(cuvect3f a)
{
float ret;
ret = ::sqrtf(a.x*a.x+a.y*a.y+a.z*a.z);
return ret;
}
__host__ __device__ cuvect3f cuvect3f_normalize(cuvect3f a)
{
cuvect3f ret;
float 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__ cuvect3f cuvect3f_proj(cuvect3f a, cuvect3f b)
{
cuvect3f ret;
cuvect3f bn = cuvect3f_normalize(b);
float m = cuvect3f_dot(a,bn);
ret = bn*m;
return ret;
}
//transposes a 3x3 (9 element) matrix
__host__ __device__ void mat3f_transpose(float *mat3inout)
{
int I,J;
float 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(float *mat3f_dest, const float *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__ float mat3f_det(float *mat3in)
{
float 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(float *mat3inout)
{
int I;
float matint[9];
float 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__ cuvect3f mat3f_mult(float *mat3in, cuvect3f cvin)
{
int I,J;
cuvect3f 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(float *matina, float *matinb, float *matout)
{
float 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__ __device__ void mat3f_hodgedual(cuvect3f vecin, float *matout)
{
matout[0 + 0*3] = 0.0f;
matout[1 + 0*3] = -vecin[2];
matout[2 + 0*3] = vecin[1];
matout[0 + 1*3] = vecin[2];
matout[1 + 1*3] = 0.0f;
matout[2 + 1*3] = -vecin[0];
matout[0 + 2*3] = -vecin[1];
matout[1 + 2*3] = vecin[0];
matout[2 + 2*3] = 0.0f;
return;
}
__host__ __device__ void mat3f_hodgedual(float *matin, cuvect3f vecout)
{
vecout[0] = 0.5*(matin[1 + 2*3] - matin[2 + 1*3]);
vecout[1] = 0.5*(matin[2 + 0*3] - matin[0 + 2*3]);
vecout[2] = 0.5*(matin[0 + 1*3] - matin[1 + 0*3]);
return;
}
//returns direction cosine rotation matrix from axis and angle
__host__ __device__ void mat3f_rot_from_axisangle(cuvect3f axis, float angle, float *matout)
{
int I;
float H[9];
float Hsq[9];
float II[9];
for(I=0;I<9;I++) II[I] = 0.0;
II[0+0*3] = 1.0;
II[1+1*3] = 1.0;
II[2+2*3] = 1.0;
axis = cuvect3f_normalize(axis);
mat3f_hodgedual(axis,H);
mat3f_mult(H,H,Hsq);
for(I=0;I<9;I++)
{
matout[I] = (II[I] + Hsq[I]) + H[I]*sinf(angle) - Hsq[I]*cosf(angle);
}
return;
}
__host__ void test_cudavectf_logic1()
{
//3 dim vector and matrix functional tests on host side
printf("3 dim vector and matrix functional tests on host side\n");
cuvect3f a,b,c;
float 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] = ((float) rand())/((float) 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 = cuvect3f(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 = cuvect3f(1,0,1);
b = cuvect3f(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 = cuvect3f_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("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);
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));
return;
}
__host__ __device__ cumat3f::cumat3f()
{
int I;
for(I=0;I<9;I++)
{
dat[I] = 0.0;
}
return;
}
__host__ __device__ cumat3f::~cumat3f()
{
int I;
for(I=0;I<9;I++)
{
dat[I] = 0.0;
}
return;
}
__host__ __device__ float& cumat3f::operator[](const int I)
{
return dat[I];
}
__host__ __device__ float& cumat3f::operator()(const int I, const int J)
{
return dat[I+3*J];
}
__host__ __device__ cumat3f cumat3f::operator+(cumat3f lhs)
{
int I;
cumat3f ret;
for(I=0;I<9;I++)
{
ret.dat[I] = this->dat[I] + lhs.dat[I];
}
return ret;
}
__host__ __device__ cumat3f cumat3f::operator-(cumat3f lhs)
{
int I;
cumat3f ret;
for(I=0;I<9;I++)
{
ret.dat[I] = this->dat[I] - lhs.dat[I];
}
return ret;
}
__host__ __device__ cumat3f cumat3f::operator*(float lhs)
{
cumat3f ret;
int I;
for(I=0;I<9;I++)
{
ret.dat[I] = this->dat[I]*lhs;
}
return ret;
}
__host__ __device__ cumat3f cumat3f::operator/(float lhs)
{
cumat3f ret;
int I;
for(I=0;I<9;I++)
{
ret.dat[I] = this->dat[I]/lhs;
}
return ret;
}
__host__ __device__ float& cumat3f::at(const int I, const int J)
{
return dat[I+3*J];
}
__host__ __device__ cuvect3f cumat3f::operator*(cuvect3f lhs)
{
cuvect3f ret = cuvect3f(0.0,0.0,0.0);
int I,J;
for(I=0;I<3;I++)
{
for(J=0;J<3;J++)
{
ret[I] = ret[I] + this->at(I,J)*lhs[J];
}
}
return ret;
}
__host__ __device__ cumat3f cumat3f::operator*(cumat3f lhs)
{
cumat3f ret;
int I,J,K;
for(I=0;I<3;I++)
{
for(J=0;J<3;J++)
{
ret(I,J) = 0.0;
for(K=0;K<3;K++)
{
ret(I,J) = ret(I,J) + this->at(I,K)*lhs(K,J);
}
}
}
return ret;
}
__host__ __device__ float cumat3f::det()
{
float ret = 0.0;
ret = ret + this->at(0,0)*this->at(1,1)*this->at(2,2);
ret = ret + this->at(0,1)*this->at(1,2)*this->at(2,0);
ret = ret + this->at(0,2)*this->at(1,0)*this->at(2,1);
ret = ret - this->at(0,0)*this->at(1,2)*this->at(2,1);
ret = ret - this->at(0,1)*this->at(1,0)*this->at(2,2);
ret = ret - this->at(0,2)*this->at(1,1)*this->at(2,0);
return ret;
}
__host__ __device__ cumat3f cumat3f::transpose()
{
cumat3f q;
int I,J;
for(I=0;I<3;I++)
{
for(J=0;J<3;J++)
{
q.at(I,J) = this->at(J,I);
}
}
return q;
}
__host__ __device__ cumat3f cumat3f::inverse()
{
cumat3f q;
float dt = q.det();
if(dt!=0.0)
{
q(0,0) = (this->at(1,1)*this->at(2,2)-this->at(1,2)*this->at(2,1))/dt;
q(0,1) = -(this->at(1,0)*this->at(2,2)-this->at(1,2)*this->at(2,0))/dt;
q(0,2) = (this->at(1,0)*this->at(2,1)-this->at(1,1)*this->at(2,0))/dt;
q(1,0) = -(this->at(0,1)*this->at(2,2)-this->at(0,2)*this->at(2,1))/dt;
q(1,1) = (this->at(0,0)*this->at(2,2)-this->at(0,2)*this->at(2,0))/dt;
q(1,2) = -(this->at(0,0)*this->at(2,1)-this->at(0,1)*this->at(2,0))/dt;
q(2,0) = (this->at(0,1)*this->at(1,2)-this->at(0,2)*this->at(1,1))/dt;
q(2,1) = -(this->at(0,0)*this->at(1,2)-this->at(0,2)*this->at(1,0))/dt;
q(2,2) = (this->at(0,0)*this->at(1,1)-this->at(0,1)*this->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;
}
};

414
src/amsculib2/cuvect4.cu Normal file
View File

@ -0,0 +1,414 @@
#include <amsculib2/amsculib2.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;
}
};

417
src/amsculib2/cuvect4f.cu Normal file
View File

@ -0,0 +1,417 @@
#include <amsculib2/amsculib2.hpp>
namespace amscuda
{
////////////
//cuvect4ff//
////////////
__host__ __device__ cuvect4f::cuvect4f()
{
x = 0.0; y = 0.0; z = 0.0; w = 0.0;
return;
}
__host__ __device__ cuvect4f::~cuvect4f()
{
x = 0.0; y = 0.0; z = 0.0; w = 0.0;
return;
}
__host__ __device__ cuvect4f::cuvect4f(float _x, float _y, float _z, float _w)
{
x = _x; y = _y; z = _z; w = _w;
return;
}
__host__ __device__ float& cuvect4f::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 float& cuvect4f::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__ cuvect4f cuvect4f::operator+(cuvect4f lhs)
{
cuvect4f 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__ cuvect4f cuvect4f::operator-(cuvect4f lhs)
{
cuvect4f 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__ cuvect4f cuvect4f::operator*(float lhs)
{
cuvect4f 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__ cuvect4f cuvect4f::operator/(float lhs)
{
cuvect4f 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__ cumat4f::cumat4f()
{
int I;
for(I=0;I<16;I++)
{
dat[I] = 0.0;
}
return;
}
__host__ __device__ cumat4f::~cumat4f()
{
int I;
for(I=0;I<16;I++)
{
dat[I] = 0.0;
}
return;
}
__host__ __device__ float& cumat4f::operator[](const int I)
{
return dat[I];
}
__host__ __device__ float& cumat4f::operator()(const int I, const int J)
{
return dat[I+4*J];
}
__host__ __device__ float& cumat4f::at(const int I, const int J)
{
return dat[I+4*J];
}
__host__ __device__ cumat4f cumat4f::operator+(cumat4f lhs)
{
cumat4f ret;
int I;
for(I=0;I<16;I++)
{
ret.dat[I] = this->dat[I] + lhs.dat[I];
}
return ret;
}
__host__ __device__ cumat4f cumat4f::operator-(cumat4f lhs)
{
cumat4f ret;
int I;
for(I=0;I<16;I++)
{
ret.dat[I] = this->dat[I] - lhs.dat[I];
}
return ret;
}
__host__ __device__ cumat4f cumat4f::operator*(float lhs)
{
cumat4f ret;
int I;
for(I=0;I<16;I++)
{
ret.dat[I] = this->dat[I]*lhs;
}
return ret;
}
__host__ __device__ cumat4f cumat4f::operator/(float lhs)
{
cumat4f ret;
int I;
for(I=0;I<16;I++)
{
ret.dat[I] = this->dat[I]/lhs;
}
return ret;
}
__host__ __device__ cuvect4f cumat4f::operator*(cuvect4f lhs)
{
cuvect4f ret = cuvect4f(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__ cumat4f cumat4f::operator*(cumat4f lhs)
{
cumat4f 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__ cumat4f cumat4f::transpose()
{
cumat4f 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__ float cumat4f::det()
{
float a00,a01,a02,a03;
float a10,a11,a12,a13;
float a20,a21,a22,a23;
float a30,a31,a32,a33;
float 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__ cumat4f minverse(cumat4f ma)
{
cumat4f mb;
float a00,a01,a02,a03;
float a10,a11,a12,a13;
float a20,a21,a22,a23;
float a30,a31,a32,a33;
float b00,b01,b02,b03;
float b10,b11,b12,b13;
float b20,b21,b22,b23;
float b30,b31,b32,b33;
float 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__ cumat4f cumat4f::inverse()
{
return minverse(*this);
}
__host__ __device__ float cuvect4f_dot(cuvect4f a, cuvect4f b)
{
float ret = 0.0;
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)
{
float ret = 0.0;
ret = ::sqrtf(cuvect4f_dot(a,a));
return ret;
}
__host__ __device__ cuvect4f cuvect4f_normalize(cuvect4f a)
{
cuvect4f ret = cuvect4f(0.0f,0.0f,0.0f,0.0f);
float nrm = cuvect4f_norm(a);
if(nrm>0.0)
ret = a/nrm;
return ret;
}
__host__ __device__ cuvect4f cuvect4f_proj(cuvect4f a, cuvect4f b)
{
cuvect4f ret;
cuvect4f bn = cuvect4f_normalize(b);
float d = cuvect4f_dot(a,bn);
ret = bn*d;
return ret;
}
}; //namespace amscuda

27
src/main.cu Normal file
View File

@ -0,0 +1,27 @@
#include <amsculib2/amsculib2.hpp>
//using namespace ams;
using namespace amscuda;
int main(int argc, char* argv[])
{
printf("AMSCULIB2: Cuda Library Tests.\n");
//test_amscuarray_1();
//test_amscumath1();
//cmp::test_cucomp64_1();
//cmp::test_cucomp128_1();
//test_amscuarray_2();
//test_dprg64();
//printf("\n");
//test_dprg32();
//test_dbuff_rand_dpr32();
test_amscurarray1();
return 0;
}

BIN
src/main.o Normal file
View File

Binary file not shown.