mirror of https://github.com/jlizier/jidt
Initial commit of GPU code.
This commit is contained in:
Pedro Martínez Mediano
parent
13952302c1
commit
5074c1b6f7
Binary file not shown.
42
build.xml
42
build.xml
|
|
@ -1,5 +1,5 @@
|
|||
<?xml version="1.0" encoding="UTF-8"?>
|
||||
<project basedir="." default="build" name="Java Information Dynamics Toolkit">
|
||||
<project basedir="." default="build" name="Java Information Dynamics Toolkit" xmlns:if="ant:if">
|
||||
<description>
|
||||
Build file for the Java Information Dynamics Toolkit
|
||||
</description>
|
||||
|
|
@ -19,18 +19,25 @@
|
|||
<property name="reports.tests" location="${unittestsouttoplevel}/reports"/>
|
||||
<property name="javadocsdir" location="javadocs"/>
|
||||
<property name="versionfile" value="version-${version}.txt"/>
|
||||
<!-- To enable GPU code, uncomment from here... -->
|
||||
<property name="enablegpu" value="true"/>
|
||||
<!-- ...to here -->
|
||||
|
||||
<path id="project.classpath">
|
||||
<pathelement location="bin"/>
|
||||
</path>
|
||||
|
||||
|
||||
<path id="apache-classpath">
|
||||
<pathelement path="./share/commons-math3-3.5.jar"/>
|
||||
</path>
|
||||
|
||||
<!-- Make required directories -->
|
||||
<target name="init" description="Create the compiled code directories">
|
||||
<mkdir dir="${bin}"/>
|
||||
<mkdir dir="${unittestsbin}"/>
|
||||
<mkdir dir="${reports.tests}"/>
|
||||
</target>
|
||||
|
||||
|
||||
<!-- Compile the java toolkit -->
|
||||
<target name="compile" depends="init" description="compile the source and unittests">
|
||||
<!-- Compile to Java 6 to provide compatibility for users with older JREs.
|
||||
|
|
@ -43,7 +50,13 @@
|
|||
ignore the warning, and I don't want to bother installing Java 6 just to compile
|
||||
like this. I'll endeavour not to use JDK 7 libraries so as not to cause
|
||||
any issues here ... -->
|
||||
<javac srcdir="${src}" destdir="${bin}" includeAntRuntime="false" target="1.6" source="1.6"/>
|
||||
<javac srcdir="${src}" destdir="${bin}" includeAntRuntime="false" target="1.6" source="1.6">
|
||||
<classpath refid="apache-classpath"/>
|
||||
</javac>
|
||||
|
||||
<!-- Compiling Cpp code -->
|
||||
<exec executable="make" dir="${src}/infodynamics/measures/continuous/kraskov/cuda/" if:set="enablegpu"/>
|
||||
|
||||
</target>
|
||||
|
||||
<!-- Jar the toolkit -->
|
||||
|
|
@ -55,17 +68,19 @@
|
|||
<!-- Compile and run the JUnit tests -->
|
||||
<target name="junit" depends="compile" description="Run the junit tests and make sure they compile">
|
||||
<!-- Compile the junit tests first -->
|
||||
<javac destdir="${unittestsbin}" includeAntRuntime="true">
|
||||
<javac destdir="${unittestsbin}" includeAntRuntime="true" debug="true">
|
||||
<!-- Need includeAntRuntime=true if you want to pick up junit.jar and ant-junit.jar in ANT_HOME/lib;
|
||||
Otherwise you will need to set the classpath to include these here. -->
|
||||
<src path="${unittestssrc}"/>
|
||||
<classpath refid="project.classpath"/>
|
||||
<classpath refid="apache-classpath"/>
|
||||
</javac>
|
||||
<!-- Run the junit tests and make sure they complete ok -->
|
||||
<junit printsummary="yes" showoutput="yes" haltonfailure="yes" haltonerror="yes" includeantruntime="true">
|
||||
<junit printsummary="yes" showoutput="yes" fork="true" forkmode="once" haltonfailure="yes" haltonerror="yes" includeantruntime="true">
|
||||
<classpath>
|
||||
<path refid="project.classpath"/>
|
||||
<pathelement path="${unittestsbin}"/>
|
||||
<pathelement location="${basedir}/clover.jar"/>
|
||||
</classpath>
|
||||
<formatter type="plain"/>
|
||||
<batchtest todir="${reports.tests}"> <!-- Writes full reports with stdout and stderr to ${reports.tests} -->
|
||||
|
|
@ -138,6 +153,21 @@
|
|||
***********************************
|
||||
-->
|
||||
|
||||
<!-- Compile and jar the toolkit with debug symbols -->
|
||||
<target name="debug" depends="init" description="compile and jar with debug symbols">
|
||||
<echo message="Compiling for debug"/>
|
||||
<javac srcdir="${src}" destdir="${bin}" includeAntRuntime="false" target="1.6" source="1.6" debug="true">
|
||||
<classpath refid="apache-classpath"/>
|
||||
</javac>
|
||||
|
||||
<!-- Compiling Cpp code -->
|
||||
<exec executable="make" dir="${src}/infodynamics/measures/continuous/kraskov/cuda/" if:set="enablegpu">
|
||||
<env key="DEBUG" value="1"/>
|
||||
</exec>
|
||||
|
||||
<jar jarfile="${jarplainname}" basedir="${bin}"/>
|
||||
</target>
|
||||
|
||||
<!-- Developer build - creates readme and version files -->
|
||||
<target name="readmefiles" description="developer: create the readme and version files from templates">
|
||||
<tstamp>
|
||||
|
|
|
|||
|
|
@ -0,0 +1,18 @@
|
|||
#!/bin/bash
|
||||
|
||||
# Make sure the latest example source file is compiled:
|
||||
# Note: this does not nulify our claim that you can compile the code
|
||||
# once and then only change the properties file - we only compile here
|
||||
# every time this script is run so that we can capture any changes the user
|
||||
# made to the demo source code. The demo as written should be compiled once
|
||||
# then one can make dynamic changes to the props file and simply
|
||||
# run the class file without recompiling.
|
||||
javac -classpath "../../infodynamics.jar" "infodynamics/demos/Example10GPUBenchmark.java"
|
||||
|
||||
# Run the example, feeding in the file names as the command line argument
|
||||
file1=WhiteNoise.txt
|
||||
file2=Correlated2D.txt
|
||||
java -classpath ".:../../infodynamics.jar" infodynamics.demos.Example10GPUBenchmark $file1 $file2
|
||||
|
||||
# Plot the results
|
||||
python plotExample10BenchmarkResults.py $file1 $file2
|
||||
|
|
@ -0,0 +1,216 @@
|
|||
/*
|
||||
* Java Information Dynamics Toolkit (JIDT)
|
||||
* Copyright (C) 2012, Joseph T. Lizier
|
||||
*
|
||||
* This program is free software: you can redistribute it and/or modify
|
||||
* it under the terms of the GNU General Public License as published by
|
||||
* the Free Software Foundation, either version 3 of the License, or
|
||||
* (at your option) any later version.
|
||||
*
|
||||
* This program is distributed in the hope that it will be useful,
|
||||
* but WITHOUT ANY WARRANTY; without even the implied warranty of
|
||||
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
|
||||
* GNU General Public License for more details.
|
||||
*
|
||||
* You should have received a copy of the GNU General Public License
|
||||
* along with this program. If not, see <http://www.gnu.org/licenses/>.
|
||||
*/
|
||||
|
||||
package infodynamics.demos;
|
||||
|
||||
import infodynamics.measures.continuous.MutualInfoCalculatorMultiVariate;
|
||||
import infodynamics.measures.continuous.kraskov.MutualInfoCalculatorMultiVariateKraskov1;
|
||||
import infodynamics.measures.continuous.kraskov.MutualInfoCalculatorMultiVariateKraskov2;
|
||||
import infodynamics.utils.MatrixUtils;
|
||||
import infodynamics.utils.RandomGenerator;
|
||||
import infodynamics.utils.ParsedProperties;
|
||||
|
||||
import java.io.FileWriter;
|
||||
import java.io.PrintWriter;
|
||||
import java.io.IOException;
|
||||
|
||||
import java.lang.Math;
|
||||
|
||||
/**
|
||||
* = Example 10 - GPU benchmark script =
|
||||
*
|
||||
* This class is used to demonstrate how the GPU module is activated and
|
||||
* includes a benchmark to test how fast it is, in comparison with its CPU
|
||||
* counterpart.
|
||||
*
|
||||
* If run from Linux, the script <code>plotExample10BenchmarkResults.py</code>
|
||||
* will use the Python library matplotlib to plot the results of the benchmark.
|
||||
*
|
||||
* @author Pedro AM Mediano
|
||||
*
|
||||
*/
|
||||
public class Example10GPUBenchmark {
|
||||
|
||||
|
||||
/**
|
||||
* Calculate MI for the given multivariate time series using the Kraskov
|
||||
* algorithm, either with CPU or GPU implementation.
|
||||
*
|
||||
* @param src Source variable for MI calculation
|
||||
* @param tgt Target variable for MI calculation
|
||||
* @param useGPU Whether or not to use the GPU for this calculation
|
||||
* @return An array with two doubles. The first component is the time in ms
|
||||
* it took to run the calculation and the second component is the actual
|
||||
* value of the calculated MI.
|
||||
*/
|
||||
public static double[] runEstimator(double[][] src, double[][] tgt, boolean useGPU) throws Exception {
|
||||
MutualInfoCalculatorMultiVariate miCalc = new MutualInfoCalculatorMultiVariateKraskov1();
|
||||
if (useGPU)
|
||||
miCalc.setProperty("USE_GPU", "true");
|
||||
miCalc.setProperty("NOISE_LEVEL_TO_ADD", "0");
|
||||
miCalc.setProperty("NORMALISE", "false");
|
||||
miCalc.setProperty("k", "4");
|
||||
miCalc.initialise(src[0].length, tgt[0].length);
|
||||
miCalc.setObservations(src, tgt);
|
||||
double[] timeAndValue = new double[2];
|
||||
long startTime = System.nanoTime();
|
||||
timeAndValue[1] = miCalc.computeAverageLocalOfObservations();
|
||||
timeAndValue[0] = (System.nanoTime() - startTime)/1000000.0;
|
||||
return timeAndValue;
|
||||
}
|
||||
|
||||
/**
|
||||
* Benchmark on a small predefined dataset, for debugging purposes.
|
||||
*/
|
||||
public static void Benchmark0() throws Exception {
|
||||
double[][] src = new double[][] {{0}, {1}, {2}, {3}, {4}, {5}, {6}};
|
||||
double[][] tgt = new double[][] {{0}, {0}, {0}, {0}, {0}, {0}, {0}};
|
||||
System.out.printf("CPU value: %f\n", runEstimator(src, tgt, false)[1]);
|
||||
System.out.printf("GPU value: %f\n", runEstimator(src, tgt, true)[1]);
|
||||
|
||||
}
|
||||
|
||||
/**
|
||||
* Benchmark on white Gaussian data of arbitrary dimension.
|
||||
*/
|
||||
public static void Benchmark1(String filename) throws Exception {
|
||||
|
||||
boolean append_to_file = false;
|
||||
FileWriter write = new FileWriter(filename, append_to_file);
|
||||
PrintWriter print_line = new PrintWriter(write);
|
||||
|
||||
int[] N_vec = new int[] {500, 1000, 2000, 4000};
|
||||
int[] D_vec = new int[] {1, 3, 5};
|
||||
|
||||
RandomGenerator rg = new RandomGenerator();
|
||||
int nb_repetitions = 5;
|
||||
|
||||
for (int i = 0; i < N_vec.length; i++) {
|
||||
for (int j = 0; j < D_vec.length; j++) {
|
||||
|
||||
int N = N_vec[i];
|
||||
int D = D_vec[j];
|
||||
|
||||
double[][] src = rg.generateNormalData(N, D, 0, 1);
|
||||
double[][] tgt = rg.generateNormalData(N, D, 0, 1);
|
||||
|
||||
double cpuTime = 0, gpuTime = 0;
|
||||
for (int r = 0; r < nb_repetitions; r++) {
|
||||
// CPU calculation
|
||||
double[] cpuVals = runEstimator(src, tgt, false);
|
||||
cpuTime += cpuVals[0];
|
||||
|
||||
// GPU calculation
|
||||
double[] gpuVals = runEstimator(src, tgt, true);
|
||||
gpuTime += gpuVals[0];
|
||||
|
||||
if (Math.abs(cpuVals[1] - gpuVals[1]) > 1e-4) {
|
||||
System.out.printf("Values differ. CPU: %f, GPU: %f\n", cpuVals[1], gpuVals[1]);
|
||||
}
|
||||
|
||||
}
|
||||
|
||||
print_line.printf("%d\t%d\t%f\t%f\n", N, D,
|
||||
cpuTime/nb_repetitions, gpuTime/nb_repetitions);
|
||||
|
||||
}
|
||||
}
|
||||
|
||||
print_line.close();
|
||||
|
||||
return;
|
||||
}
|
||||
|
||||
|
||||
|
||||
/**
|
||||
* Benchmark on 2D correlated Gaussian.
|
||||
*
|
||||
* Reference:
|
||||
*
|
||||
* http://math.stackexchange.com/questions/446093/generate-correlated-normal-random-variables
|
||||
*
|
||||
*/
|
||||
public static void Benchmark2(String filename) throws Exception {
|
||||
|
||||
boolean append_to_file = false;
|
||||
FileWriter write = new FileWriter(filename, append_to_file);
|
||||
PrintWriter print_line = new PrintWriter(write);
|
||||
|
||||
int[] N_vec = new int[] {500, 1000, 2000, 4000};
|
||||
double[] corr_vec = new double[] {0, 0.2, 0.4, 0.6};
|
||||
|
||||
RandomGenerator rg = new RandomGenerator();
|
||||
int nb_repetitions = 5;
|
||||
|
||||
for (int i = 0; i < N_vec.length; i++) {
|
||||
for (int j = 0; j < corr_vec.length; j++) {
|
||||
|
||||
int N = N_vec[i];
|
||||
double corr = corr_vec[j];
|
||||
|
||||
double[] r1 = rg.generateNormalData(N, 0, 1);
|
||||
double[] r2 = rg.generateNormalData(N, 0, 1);
|
||||
double[][] src = new double[N][1];
|
||||
double[][] tgt = new double[N][1];
|
||||
for (int t = 0; t < N; t++) {
|
||||
src[t][0] = r1[t];
|
||||
tgt[t][0] = corr*src[t][0] + Math.sqrt(1 - corr*corr)*r2[t];
|
||||
}
|
||||
|
||||
double cpuTime = 0, gpuTime = 0;
|
||||
for (int r = 0; r < nb_repetitions; r++) {
|
||||
// CPU calculation
|
||||
double[] cpuVals = runEstimator(src, tgt, false);
|
||||
cpuTime += cpuVals[0];
|
||||
|
||||
// GPU calculation
|
||||
double[] gpuVals = runEstimator(src, tgt, true);
|
||||
gpuTime += gpuVals[0];
|
||||
|
||||
if (Math.abs(cpuVals[1] - gpuVals[1]) > 1e-4) {
|
||||
System.out.printf("Values differ. CPU: %f, GPU: %f\n", cpuVals[1], gpuVals[1]);
|
||||
}
|
||||
}
|
||||
|
||||
// Save time in milliseconds
|
||||
print_line.printf("%d\t%f\t%f\t%f\n", N, corr,
|
||||
cpuTime/nb_repetitions, gpuTime/nb_repetitions);
|
||||
|
||||
}
|
||||
}
|
||||
|
||||
print_line.close();
|
||||
|
||||
return;
|
||||
}
|
||||
|
||||
/**
|
||||
* Run several of the benchmarks above and write the results to a
|
||||
* .txt file.
|
||||
*
|
||||
* @param args List of file names to save the results of the benchmark.
|
||||
*/
|
||||
public static void main(String[] args) throws Exception {
|
||||
// Benchmark0();
|
||||
Benchmark1(args[0]);
|
||||
Benchmark2(args[1]);
|
||||
return;
|
||||
}
|
||||
|
||||
}
|
||||
|
|
@ -0,0 +1,47 @@
|
|||
#!/usr/bin/python
|
||||
|
||||
"""
|
||||
plotExample10BenchmarkResults.py: Show 3D plots with the result of benchmarking
|
||||
the CPU and GPU implementations of the KSG algorithm for mutual information.
|
||||
"""
|
||||
|
||||
__author__ = "Pedro AM Mediano"
|
||||
__email__ = "pmediano@imperial.ac.uk"
|
||||
__license__ = "GPL"
|
||||
|
||||
import sys
|
||||
import numpy as np
|
||||
import matplotlib.pyplot as pl
|
||||
from mpl_toolkits.mplot3d import Axes3D
|
||||
|
||||
## Load files
|
||||
filename1 = sys.argv[1]
|
||||
filename2 = sys.argv[2]
|
||||
res1 = np.loadtxt(filename1)
|
||||
res2 = np.loadtxt(filename2)
|
||||
|
||||
## There is a bug with the legends of 3D scatterplots, so we have to save
|
||||
# the plot style for the legend separately
|
||||
img1 = pl.Line2D([0],[0], linestyle="none", c='r', marker='o')
|
||||
img2 = pl.Line2D([0],[0], linestyle="none", c='b', marker='^')
|
||||
|
||||
## Plot times for first benchmark
|
||||
fig1 = pl.figure()
|
||||
ax1 = fig1.add_subplot(111, projection='3d')
|
||||
ax1.scatter(np.log10(res1[:,0]), res1[:,1], res1[:,2], c='r', marker='o')
|
||||
ax1.scatter(np.log10(res1[:,0]), res1[:,1], res1[:,3], c='b', marker='^')
|
||||
ax1.legend([img1, img2], ['CPU', 'GPU'], numpoints = 1)
|
||||
ax1.set_xlabel(r'$\log_{10} N$'); ax1.set_ylabel('D'); ax1.set_zlabel('Time [ms]');
|
||||
ax1.set_title(filename1.partition(".")[0])
|
||||
|
||||
## Plot times for second benchmark
|
||||
fig2 = pl.figure()
|
||||
ax2 = fig2.add_subplot(111, projection='3d')
|
||||
ax2.scatter(np.log10(res2[:,0]), res2[:,1], res2[:,2], c='r', marker='o')
|
||||
ax2.scatter(np.log10(res2[:,0]), res2[:,1], res2[:,3], c='b', marker='^')
|
||||
ax2.legend([img1, img2], ['CPU', 'GPU'], numpoints = 1)
|
||||
ax2.set_xlabel(r'$\log_{10} N$'); ax2.set_ylabel(r'$\rho$'); ax2.set_zlabel('Time [ms]');
|
||||
ax2.set_title(filename2.partition(".")[0])
|
||||
|
||||
pl.show()
|
||||
|
||||
File diff suppressed because it is too large
Load Diff
|
|
@ -23,7 +23,9 @@ import java.util.PriorityQueue;
|
|||
|
||||
import infodynamics.measures.continuous.MutualInfoCalculatorMultiVariate;
|
||||
import infodynamics.utils.MathsUtils;
|
||||
import infodynamics.utils.MatrixUtils;
|
||||
import infodynamics.utils.NeighbourNodeData;
|
||||
import infodynamics.utils.EuclideanUtils;
|
||||
|
||||
/**
|
||||
* <p>Computes the differential mutual information of two given multivariate sets of
|
||||
|
|
@ -80,7 +82,7 @@ public class MutualInfoCalculatorMultiVariateKraskov1
|
|||
// First element in the PQ is the kth NN,
|
||||
// and epsilon = kthNnData.distance
|
||||
NeighbourNodeData kthNnData = nnPQ.poll();
|
||||
|
||||
|
||||
// Count the number of points whose x distance is less
|
||||
// than eps, and whose y distance is less than
|
||||
// epsilon = kthNnData.distance
|
||||
|
|
@ -88,7 +90,7 @@ public class MutualInfoCalculatorMultiVariateKraskov1
|
|||
t, kthNnData.distance, dynCorrExclTime);
|
||||
int n_y = nnSearcherDest.countPointsStrictlyWithinR(
|
||||
t, kthNnData.distance, dynCorrExclTime);
|
||||
|
||||
|
||||
sumNx += n_x;
|
||||
sumNy += n_y;
|
||||
// And take the digammas:
|
||||
|
|
@ -103,7 +105,7 @@ public class MutualInfoCalculatorMultiVariateKraskov1
|
|||
localMi[t-startTimePoint] = digammaK - digammaNxPlusOne - digammaNyPlusOne + digammaN;
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
if (debug) {
|
||||
Calendar rightNow2 = Calendar.getInstance();
|
||||
long endTime = rightNow2.getTimeInMillis();
|
||||
|
|
@ -119,4 +121,5 @@ public class MutualInfoCalculatorMultiVariateKraskov1
|
|||
return new double[] {sumDiGammas, sumNx, sumNy};
|
||||
}
|
||||
}
|
||||
|
||||
}
|
||||
|
|
|
|||
|
|
@ -129,4 +129,5 @@ public class MutualInfoCalculatorMultiVariateKraskov2
|
|||
return new double[] {sumDiGammas, sumNx, sumNy};
|
||||
}
|
||||
}
|
||||
|
||||
}
|
||||
|
|
|
|||
|
|
@ -0,0 +1,39 @@
|
|||
|
||||
## Set these if the CUDA Toolkit is not in your PATH
|
||||
# Location of the CUDA Toolkit binaries and libraries
|
||||
# CUDA_PATH ?= /usr/local/cuda
|
||||
# CUDA_INC_PATH ?= $(CUDA_PATH)/include
|
||||
# CUDA_LIB_PATH ?= $(CUDA_PATH)/lib
|
||||
# CUDA_BIN_PATH ?= $(CUDA_PATH)/bin
|
||||
|
||||
# Common binaries and flags. TODO: the architecture flag may not be necessary
|
||||
NVCC ?= nvcc
|
||||
GCC ?= g++
|
||||
CCFLAGS ?= -O3
|
||||
# NVCCFLAGS ?= -O3 -arch=sm_30
|
||||
NVCCFLAGS ?= -O3
|
||||
|
||||
ifdef DEBUG
|
||||
CCFLAGS += -g
|
||||
NVCCFLAGS += -g -G
|
||||
endif
|
||||
|
||||
# Common includes and paths for CUDA. This assumes the CUDA toolkit is in PATH
|
||||
INCLUDES := -I. -I/usr/lib/jvm/java-7-openjdk-amd64/include
|
||||
LDFLAGS += -L. -lcuda -lcudart
|
||||
|
||||
# Target rules
|
||||
all: libKraskov.so
|
||||
|
||||
gpuKnnLibrary.o: gpuKnnLibrary.cu helperfunctions.cu gpuKnnBF_kernel.cu
|
||||
${NVCC} ${NVCCFLAGS} ${INCLUDES} -Xcompiler -fPIC -c gpuKnnLibrary.cu
|
||||
|
||||
libgpuKnnLibrary.a: gpuKnnLibrary.o
|
||||
${AR} -r libgpuKnnLibrary.a gpuKnnLibrary.o
|
||||
|
||||
libKraskov.so: libgpuKnnLibrary.a kraskovCuda.c
|
||||
${NVCC} ${NVCCFLAGS} ${INCLUDES} -Xcompiler -fPIC -shared -o libKraskov.so kraskovCuda.c ${LDFLAGS} -lgpuKnnLibrary
|
||||
|
||||
clean:
|
||||
rm -f *.o *.so *.a
|
||||
|
||||
|
|
@ -0,0 +1,120 @@
|
|||
/*************************************
|
||||
An ANSI-C implementation of the digamma-function for real arguments based
|
||||
on the Chebyshev expansion proposed in appendix E of
|
||||
http://arXiv.org/abs/math.CA/0403344 . This is identical to the implementation
|
||||
by Jet Wimp, Math. Comp. vol 15 no 74 (1961) pp 174 (see Table 1).
|
||||
For other implementations see
|
||||
the GSL implementation for Psi(Digamma) in
|
||||
http://www.gnu.org/software/gsl/manual/html_node/Psi-_0028Digamma_0029-Function.html
|
||||
|
||||
Richard J. Mathar, 2005-11-24
|
||||
***************************************/
|
||||
#include <math.h>
|
||||
|
||||
#ifndef M_PIl
|
||||
/** The constant Pi in high precision */
|
||||
#define M_PIl 3.1415926535897932384626433832795029L
|
||||
#endif
|
||||
#ifndef M_GAMMAl
|
||||
/** Euler's constant in high precision */
|
||||
#define M_GAMMAl 0.5772156649015328606065120900824024L
|
||||
#endif
|
||||
#ifndef M_LN2l
|
||||
/** the natural logarithm of 2 in high precision */
|
||||
#define M_LN2l 0.6931471805599453094172321214581766L
|
||||
#endif
|
||||
|
||||
/** The digamma function in long double precision.
|
||||
* @param x the real value of the argument
|
||||
* @return the value of the digamma (psi) function at that point
|
||||
* @author Richard J. Mathar
|
||||
* @since 2005-11-24
|
||||
*/
|
||||
#ifdef __cplusplus
|
||||
extern "C" {
|
||||
#endif
|
||||
long double cpuDigamma(long double x)
|
||||
{
|
||||
int n;
|
||||
/* force into the interval 1..3 */
|
||||
if( x < 0.0L )
|
||||
return cpuDigamma(1.0L-x)+M_PIl/tanl(M_PIl*(1.0L-x)) ; /* reflection formula */
|
||||
else if( x < 1.0L )
|
||||
return cpuDigamma(1.0L+x)-1.0L/x ;
|
||||
else if ( x == 1.0L)
|
||||
return -M_GAMMAl ;
|
||||
else if ( x == 2.0L)
|
||||
return 1.0L-M_GAMMAl ;
|
||||
else if ( x == 3.0L)
|
||||
return 1.5L-M_GAMMAl ;
|
||||
else if ( x > 3.0L)
|
||||
/* duplication formula */
|
||||
return 0.5L*(cpuDigamma(x/2.0L)+cpuDigamma((x+1.0L)/2.0L))+M_LN2l ;
|
||||
else {
|
||||
/* Just for your information, the following lines contain
|
||||
* the Maple source code to re-generate the table that is
|
||||
* eventually becoming the Kncoe[] array below
|
||||
* interface(prettyprint=0) :
|
||||
* Digits := 63 :
|
||||
* r := 0 :
|
||||
*
|
||||
* for l from 1 to 60 do
|
||||
* d := binomial(-1/2,l) :
|
||||
* r := r+d*(-1)^l*(Zeta(2*l+1) -1) ;
|
||||
* evalf(r) ;
|
||||
* print(%,evalf(1+Psi(1)-r)) ;
|
||||
*o d :
|
||||
*
|
||||
* for N from 1 to 28 do
|
||||
* r := 0 :
|
||||
* n := N-1 :
|
||||
*
|
||||
* for l from iquo(n+3,2) to 70 do
|
||||
* d := 0 :
|
||||
* for s from 0 to n+1 do
|
||||
* d := d+(-1)^s*binomial(n+1,s)*binomial((s-1)/2,l) :
|
||||
* od :
|
||||
* if 2*l-n > 1 then
|
||||
* r := r+d*(-1)^l*(Zeta(2*l-n) -1) :
|
||||
* fi :
|
||||
* od :
|
||||
* print(evalf((-1)^n*2*r)) ;
|
||||
*od :
|
||||
*quit :
|
||||
*/
|
||||
static long double Kncoe[] = { .30459198558715155634315638246624251L,
|
||||
.72037977439182833573548891941219706L, -.12454959243861367729528855995001087L,
|
||||
.27769457331927827002810119567456810e-1L, -.67762371439822456447373550186163070e-2L,
|
||||
.17238755142247705209823876688592170e-2L, -.44817699064252933515310345718960928e-3L,
|
||||
.11793660000155572716272710617753373e-3L, -.31253894280980134452125172274246963e-4L,
|
||||
.83173997012173283398932708991137488e-5L, -.22191427643780045431149221890172210e-5L,
|
||||
.59302266729329346291029599913617915e-6L, -.15863051191470655433559920279603632e-6L,
|
||||
.42459203983193603241777510648681429e-7L, -.11369129616951114238848106591780146e-7L,
|
||||
.304502217295931698401459168423403510e-8L, -.81568455080753152802915013641723686e-9L,
|
||||
.21852324749975455125936715817306383e-9L, -.58546491441689515680751900276454407e-10L,
|
||||
.15686348450871204869813586459513648e-10L, -.42029496273143231373796179302482033e-11L,
|
||||
.11261435719264907097227520956710754e-11L, -.30174353636860279765375177200637590e-12L,
|
||||
.80850955256389526647406571868193768e-13L, -.21663779809421233144009565199997351e-13L,
|
||||
.58047634271339391495076374966835526e-14L, -.15553767189204733561108869588173845e-14L,
|
||||
.41676108598040807753707828039353330e-15L, -.11167065064221317094734023242188463e-15L } ;
|
||||
|
||||
register long double Tn_1 = 1.0L ; /* T_{n-1}(x), started at n=1 */
|
||||
register long double Tn = x-2.0L ; /* T_{n}(x) , started at n=1 */
|
||||
register long double resul = Kncoe[0] + Kncoe[1]*Tn ;
|
||||
|
||||
x -= 2.0L ;
|
||||
|
||||
for(n = 2 ; n < sizeof(Kncoe)/sizeof(long double) ;n++)
|
||||
{
|
||||
const long double Tn1 = 2.0L * x * Tn - Tn_1 ; /* Chebyshev recursion, Eq. 22.7.4 Abramowitz-Stegun */
|
||||
resul += Kncoe[n]*Tn1 ;
|
||||
Tn_1 = Tn ;
|
||||
Tn = Tn1 ;
|
||||
}
|
||||
return resul ;
|
||||
}
|
||||
}
|
||||
#ifdef __cplusplus
|
||||
}
|
||||
#endif
|
||||
|
||||
|
|
@ -0,0 +1,455 @@
|
|||
#ifndef _TEMPLATE_KERNEL_H_
|
||||
#define _TEMPLATE_KERNEL_H_
|
||||
|
||||
//#include <cpugpuKnn_common.h>
|
||||
#ifndef INFINITY
|
||||
#define INFINITY 0x7F800000
|
||||
#endif
|
||||
|
||||
/**
|
||||
* Typedef function pointer for arbitrary norm functions.
|
||||
*/
|
||||
typedef float (*normFunction_t)(const float* g_uquery, const float* g_vpoint,
|
||||
int pointdim, int signallength);
|
||||
|
||||
/**
|
||||
* Calculate max norm (L-inf) between two points.
|
||||
*/
|
||||
__device__ __host__ float
|
||||
maxMetricPoints(const float* g_uquery, const float* g_vpoint, int pointdim, int signallength){
|
||||
float r_u1;
|
||||
float r_v1;
|
||||
float r_d1,r_dim=0;
|
||||
|
||||
r_dim=0;
|
||||
for(int d=0; d<pointdim; d++){
|
||||
r_u1 = *(g_uquery+d*signallength);
|
||||
r_v1 = *(g_vpoint+d*signallength);
|
||||
r_d1 = r_v1 - r_u1;
|
||||
r_d1 = r_d1 < 0? -r_d1: r_d1; //abs
|
||||
r_dim= r_dim < r_d1? r_d1: r_dim;
|
||||
}
|
||||
return r_dim;
|
||||
}
|
||||
|
||||
/**
|
||||
* Calculate squared Euclidean norm (L2) between two points (note we keep the
|
||||
* distance squared to avoid taking sqrt all the time).
|
||||
*/
|
||||
__device__ __host__ float
|
||||
squareMetricPoints(const float* g_uquery, const float* g_vpoint, int pointdim, int signallength){
|
||||
float r_u1;
|
||||
float r_v1;
|
||||
float r_d1, r_dim = 0.0f;
|
||||
|
||||
for (int d = 0; d < pointdim; d++) {
|
||||
r_u1 = *(g_uquery+d*signallength);
|
||||
r_v1 = *(g_vpoint+d*signallength);
|
||||
r_d1 = r_v1 - r_u1;
|
||||
r_dim += r_d1 * r_d1;
|
||||
}
|
||||
return r_dim;
|
||||
}
|
||||
|
||||
__device__ normFunction_t pMaxNorm = maxMetricPoints;
|
||||
__device__ normFunction_t pSquareNorm = squareMetricPoints;
|
||||
|
||||
/**
|
||||
* Insert point in current list of nearest neighbours.
|
||||
*/
|
||||
__device__ float
|
||||
insertPointKlist(int kth, float distance, int indexv,float* kdistances, int* kindexes){
|
||||
int k=0;
|
||||
while( (distance>*(kdistances+k)) && (k<kth-1)){k++;}
|
||||
//Move value to the next
|
||||
for(int k2=kth-1;k2>k;k2--){
|
||||
*(kdistances+k2)=*(kdistances+k2-1);
|
||||
*(kindexes+k2)=*(kindexes+k2-1);
|
||||
}
|
||||
//Replace
|
||||
*(kdistances+k)=distance;
|
||||
*(kindexes+k)=indexv;
|
||||
|
||||
//printf("\n -> Modificacion pila: %.f %.f. New max distance: %.f", *kdistances, *(kdistances+1), *(kdistances+kth-1));
|
||||
return *(kdistances+kth-1);
|
||||
}
|
||||
|
||||
|
||||
/*
|
||||
* Main KNN kernel. Find the nearest k neighbours to each point according
|
||||
* to the supplied norm function.
|
||||
*/
|
||||
__global__ void
|
||||
kernelKNNshared(const float* g_uquery, const float* g_vpointset,
|
||||
int *g_indexes, float* g_distances, const int pointdim,
|
||||
const int triallength, const int signallength, const int kth,
|
||||
const int exclude, normFunction_t *normFunction) {
|
||||
|
||||
// Shared memory
|
||||
extern __shared__ char array[];
|
||||
float *kdistances;
|
||||
int *kindexes;
|
||||
kdistances = (float*)array;
|
||||
kindexes = (int*)array+kth*blockDim.x;
|
||||
|
||||
const unsigned int tid = threadIdx.x + blockDim.x*blockIdx.x;
|
||||
const unsigned int itrial = tid / triallength; // indextrial
|
||||
|
||||
if(tid<signallength){
|
||||
|
||||
for(int k=0;k<kth;k++){
|
||||
kdistances[threadIdx.x*kth+k] = INFINITY;
|
||||
}
|
||||
|
||||
__syncthreads();
|
||||
|
||||
float r_kdist=INFINITY;
|
||||
unsigned int indexi = tid-triallength*itrial;
|
||||
|
||||
for(int t=0; t<triallength; t++){
|
||||
int indexu = tid;
|
||||
int indexv = (t + itrial*triallength);
|
||||
int condition1=indexi-exclude;
|
||||
int condition2=indexi+exclude;
|
||||
if((t<condition1)||(t>condition2)){
|
||||
float temp_dist = normFunction[0](g_uquery+indexu, g_vpointset+indexv,pointdim, signallength);
|
||||
if(temp_dist <= r_kdist){
|
||||
r_kdist = insertPointKlist(kth,temp_dist,t,kdistances+threadIdx.x*kth,kindexes+threadIdx.x*kth);
|
||||
}
|
||||
}
|
||||
//printf("tid:%d indexes: %d, %d distances: %.f %.f\n",tid, *kindexes, *(kindexes+1), *kdistances, *(kdistances+1));
|
||||
}
|
||||
|
||||
__syncthreads();
|
||||
// COPY TO GLOBAL MEMORY
|
||||
for (int k = 0; k < kth; k++) {
|
||||
g_indexes[tid+k*signallength] = kindexes[threadIdx.x*kth+k];
|
||||
g_distances[tid+k*signallength] = kdistances[threadIdx.x*kth+k];//*(kdistances+k);
|
||||
}
|
||||
}
|
||||
|
||||
}
|
||||
|
||||
|
||||
/*
|
||||
* Range search for one data point using bruteforce.
|
||||
*/
|
||||
__global__ void
|
||||
kernelBFRSshared(const float* g_uquery, const float* g_vpointset, int *g_npoints, int pointdim, int triallength, int signallength, int exclude, float radius)
|
||||
{
|
||||
|
||||
// Shared memory
|
||||
extern __shared__ char array[];
|
||||
int *s_npointsrange;
|
||||
s_npointsrange = (int*)array;
|
||||
|
||||
const unsigned int tid = threadIdx.x + blockDim.x*blockIdx.x;
|
||||
const unsigned int itrial = tid / triallength; // indextrial
|
||||
|
||||
if(tid<signallength){
|
||||
|
||||
s_npointsrange[threadIdx.x] = 0;
|
||||
__syncthreads();
|
||||
|
||||
|
||||
unsigned int indexi = tid-triallength*itrial;
|
||||
for(int t=0; t<triallength; t++){
|
||||
int indexu = tid;
|
||||
int indexv = (t + itrial*triallength);
|
||||
int condition1=indexi-exclude;
|
||||
int condition2=indexi+exclude;
|
||||
if((t<condition1)||(t>condition2)){
|
||||
float temp_dist = maxMetricPoints(g_uquery+indexu, g_vpointset+indexv,pointdim, signallength);
|
||||
if(temp_dist <= radius){
|
||||
s_npointsrange[threadIdx.x]++;
|
||||
}
|
||||
}
|
||||
|
||||
}
|
||||
|
||||
__syncthreads();
|
||||
|
||||
// COPY TO GLOBAL MEMORY
|
||||
g_npoints[tid] = s_npointsrange[threadIdx.x];
|
||||
|
||||
}
|
||||
}
|
||||
|
||||
/*
|
||||
* Range search using bruteforce in multiple GPUs
|
||||
*/
|
||||
__global__ void
|
||||
kernelBFRSMultishared(const float* g_uquery, const float* g_vpointset, int *g_npoints, int pointdim, int triallength, int signallength, int exclude, const float* vecradius)
|
||||
{
|
||||
|
||||
// shared memory
|
||||
extern __shared__ char array[];
|
||||
int *s_npointsrange;
|
||||
s_npointsrange = (int*)array;
|
||||
float radius=0;
|
||||
const unsigned int tid = threadIdx.x + blockDim.x*blockIdx.x;
|
||||
const unsigned int itrial = tid / triallength; // indextrial
|
||||
|
||||
if(tid<signallength){
|
||||
|
||||
s_npointsrange[threadIdx.x] = 0;
|
||||
__syncthreads();
|
||||
|
||||
radius = *(vecradius+itrial);
|
||||
unsigned int indexi = tid-triallength*itrial;
|
||||
for(int t=0; t<triallength; t++){
|
||||
int indexu = tid;
|
||||
int indexv = (t + itrial*triallength);
|
||||
int condition1=indexi-exclude;
|
||||
int condition2=indexi+exclude;
|
||||
if((t<condition1)||(t>condition2)){
|
||||
float temp_dist = maxMetricPoints(g_uquery+indexu, g_vpointset+indexv,pointdim, signallength);
|
||||
if(temp_dist <= radius){
|
||||
s_npointsrange[threadIdx.x]++;
|
||||
}
|
||||
}
|
||||
|
||||
}
|
||||
|
||||
__syncthreads();
|
||||
//printf("\ntid:%d npoints: %d\n",tid, s_npointsrange[threadIdx.x]);
|
||||
//COPY TO GLOBAL MEMORY
|
||||
g_npoints[tid] = s_npointsrange[threadIdx.x];
|
||||
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
/*
|
||||
* Range search for all data points using bruteforce.
|
||||
*/
|
||||
__global__ void
|
||||
kernelBFRSAllshared(const float* g_uquery, const float* g_vpointset,
|
||||
int *g_npoints, int pointdim, int triallength, int signallength,
|
||||
int exclude, const float* vecradius, normFunction_t *normFunction) {
|
||||
|
||||
// Shared memory
|
||||
extern __shared__ char array[];
|
||||
int *s_npointsrange;
|
||||
s_npointsrange = (int *) array;
|
||||
float radius = 0;
|
||||
const unsigned int tid = threadIdx.x + blockDim.x*blockIdx.x;
|
||||
const unsigned int itrial = tid / triallength; // indextrial
|
||||
|
||||
if(tid<signallength){
|
||||
|
||||
s_npointsrange[threadIdx.x] = 0;
|
||||
__syncthreads();
|
||||
|
||||
radius = *(vecradius+tid);
|
||||
unsigned int indexi = tid-triallength*itrial;
|
||||
for (int t=0; t<triallength; t++){
|
||||
int indexu = tid;
|
||||
int indexv = (t + itrial*triallength);
|
||||
int condition1=indexi-exclude;
|
||||
int condition2=indexi+exclude;
|
||||
if((t<condition1)||(t>condition2)){
|
||||
float temp_dist = normFunction[0](g_uquery+indexu, g_vpointset+indexv,pointdim, signallength);
|
||||
// PEDRO: For KSG algorithm 1 this should be strictly less than R, and in TRENTOOL code it's less or equal. It's a float comparison, so I don't think it matters anyway.
|
||||
if(temp_dist < radius){
|
||||
s_npointsrange[threadIdx.x]++;
|
||||
}
|
||||
}
|
||||
|
||||
}
|
||||
|
||||
__syncthreads();
|
||||
|
||||
//COPY TO GLOBAL MEMORY
|
||||
g_npoints[tid] = s_npointsrange[threadIdx.x];
|
||||
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
|
||||
__device__ void digammaXp1(double *y);
|
||||
|
||||
/**
|
||||
* Optimized reduction kernel, obtained from Mark Harris' lecture on
|
||||
* GPU optimization:
|
||||
* https://docs.nvidia.com/cuda/samples/6_Advanced/reduction/doc/reduction.pdf
|
||||
*/
|
||||
template <unsigned int blockSize>
|
||||
// __device__ void warpReduce(volatile int *sdata, unsigned int tid) {
|
||||
__device__ void warpReduce(volatile float *sdata, unsigned int tid) {
|
||||
if (blockSize >= 64) sdata[tid] += sdata[tid + 32];
|
||||
if (blockSize >= 32) sdata[tid] += sdata[tid + 16];
|
||||
if (blockSize >= 16) sdata[tid] += sdata[tid + 8];
|
||||
if (blockSize >= 8) sdata[tid] += sdata[tid + 4];
|
||||
if (blockSize >= 4) sdata[tid] += sdata[tid + 2];
|
||||
if (blockSize >= 2) sdata[tid] += sdata[tid + 1];
|
||||
}
|
||||
|
||||
// template <unsigned int blockSize>
|
||||
// __global__ void reduce6(int *g_idata, int *g_odata, unsigned int n) {
|
||||
__global__ void reduce6(int *g_nx, int *g_ny, float *g_odata, unsigned int n) {
|
||||
// extern __shared__ int sdata[];
|
||||
extern __shared__ float sdata[];
|
||||
unsigned int tid = threadIdx.x;
|
||||
const unsigned int blockSize = 512;
|
||||
unsigned int i = blockIdx.x*(blockSize*2) + tid;
|
||||
unsigned int gridSize = blockSize*2*gridDim.x;
|
||||
sdata[tid] = 0;
|
||||
|
||||
while (i<n) {
|
||||
double dgX1 = (double) g_nx[i];
|
||||
double dgY1 = (double) g_ny[i];
|
||||
double dgX2 = (double) g_nx[i+blockSize];
|
||||
double dgY2 = (double) g_ny[i+blockSize];
|
||||
|
||||
digammaXp1(&dgX1);
|
||||
digammaXp1(&dgY1);
|
||||
digammaXp1(&dgX2);
|
||||
digammaXp1(&dgY2);
|
||||
|
||||
sdata[tid] = dgX1 + dgY1 + dgX2 + dgY2;
|
||||
i += gridSize;
|
||||
}
|
||||
__syncthreads();
|
||||
|
||||
if (blockSize >= 512) { if (tid < 256) { sdata[tid] += sdata[tid + 256]; } __syncthreads(); }
|
||||
if (blockSize >= 256) { if (tid < 128) { sdata[tid] += sdata[tid + 128]; } __syncthreads(); }
|
||||
if (blockSize >= 128) { if (tid < 64) { sdata[tid] += sdata[tid + 64]; } __syncthreads(); }
|
||||
if (tid < 32) warpReduce<blockSize>(sdata, tid);
|
||||
if (tid == 0) g_odata[blockIdx.x] = sdata[0];
|
||||
}
|
||||
|
||||
|
||||
__global__ void reduce1(int *g_nx, int *g_ny, float *g_odata, unsigned int N) {
|
||||
extern __shared__ float sdata2[];
|
||||
// each thread loads one element from global to shared mem
|
||||
unsigned int tid = threadIdx.x;
|
||||
unsigned int i = blockIdx.x*blockDim.x + threadIdx.x;
|
||||
|
||||
if (i < N) {
|
||||
double dgX = g_nx[i];
|
||||
double dgY = g_ny[i];
|
||||
digammaXp1(&dgX);
|
||||
digammaXp1(&dgY);
|
||||
sdata2[tid] = dgX + dgY;
|
||||
} else {
|
||||
sdata2[tid] = 0;
|
||||
}
|
||||
__syncthreads();
|
||||
|
||||
// do reduction in shared mem
|
||||
for (unsigned int s=1; s < blockDim.x; s *= 2) {
|
||||
if (tid % (2*s) == 0) {
|
||||
sdata2[tid] += sdata2[tid + s];
|
||||
}
|
||||
__syncthreads();
|
||||
}
|
||||
// write result for this block to global mem
|
||||
if (tid == 0) g_odata[blockIdx.x] = sdata2[0];
|
||||
}
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
/**
|
||||
* Taken from NVIDIA dev forum:
|
||||
* https://devtalk.nvidia.com/default/topic/516516/kernel-launch-failure-in-matlab/
|
||||
*/
|
||||
__device__ void digammaXp1(double *y) {
|
||||
|
||||
// double x = *y;
|
||||
double x = *y + 1;
|
||||
double neginf = -INFINITY;
|
||||
const double c = 12,
|
||||
digamma1 = -0.57721566490153286,
|
||||
trigamma1 = 1.6449340668482264365, /* pi^2/6 */
|
||||
s = 1e-6,
|
||||
s3 = 1./12,
|
||||
s4 = 1./120,
|
||||
s5 = 1./252,
|
||||
s6 = 1./240,
|
||||
s7 = 1./132;
|
||||
// s8 = 691./32760,
|
||||
// s9 = 1./12,
|
||||
// s10 = 3617./8160;
|
||||
double result;
|
||||
|
||||
/* Illegal arguments */
|
||||
if((x == neginf) || isnan(x)) {
|
||||
*y = NAN;
|
||||
return;
|
||||
}
|
||||
|
||||
/* Singularities */
|
||||
if((x <= 0) && (floor(x) == x)) {
|
||||
*y = neginf;
|
||||
return;
|
||||
}
|
||||
|
||||
/* Negative values */
|
||||
|
||||
/* Use the reflection formula (Jeffrey 11.1.6):
|
||||
* digamma(-x) = digamma(x+1) + pi*cot(pi*x)
|
||||
*
|
||||
* This is related to the identity
|
||||
* digamma(-x) = digamma(x+1) - digamma(z) + digamma(1-z)
|
||||
* where z is the fractional part of x
|
||||
* For example:
|
||||
|
||||
* digamma(-3.1) = 1/3.1 + 1/2.1 + 1/1.1 + 1/0.1 + digamma(1-0.1)
|
||||
* = digamma(4.1) - digamma(0.1) + digamma(1-0.1)
|
||||
* Then we use
|
||||
* digamma(1-z) - digamma(z) = pi*cot(pi*z)
|
||||
*
|
||||
|
||||
if(x < 0) {
|
||||
*p = digamma(p,1-x) + M_PI/tan(-M_PI*x);
|
||||
return;
|
||||
}
|
||||
*/
|
||||
|
||||
/* Use Taylor series if argument <= S */
|
||||
if(x <= s) {
|
||||
*y = digamma1 - 1/x + trigamma1*x;
|
||||
return;
|
||||
}
|
||||
|
||||
/* Reduce to digamma(X + N) where (X + N) >= C */
|
||||
result = 0;
|
||||
while(x < c) {
|
||||
result -= 1/x;
|
||||
x++;
|
||||
}
|
||||
|
||||
/* Use de Moivre's expansion if argument >= C */
|
||||
/* This expansion can be computed in Maple via asympt(Psi(x),x) */
|
||||
if(x >= c) {
|
||||
double r = 1/x, t;
|
||||
result += log(x) - 0.5*r;
|
||||
r *= r;
|
||||
#if 0
|
||||
result -= r * (s3 - r * (s4 - r * (s5 - r * (s6 - r * s7))));
|
||||
#else
|
||||
/* this version for lame compilers */
|
||||
t = (s5 - r * (s6 - r * s7));
|
||||
result -= r * (s3 - r * (s4 - r * t));
|
||||
#endif
|
||||
}
|
||||
|
||||
/* assign the result to the pointer*/
|
||||
*y = result;
|
||||
return;
|
||||
|
||||
}
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
#endif // #ifndef _TEMPLATE_KERNEL_H_
|
||||
|
||||
|
||||
|
|
@ -0,0 +1,278 @@
|
|||
#include "gpuKnnBF_kernel.cu"
|
||||
#include <stdio.h>
|
||||
#include "helperfunctions.cu"
|
||||
|
||||
#ifdef __cplusplus
|
||||
extern "C" {
|
||||
#endif
|
||||
int cudaFindKnn(int* h_bf_indexes, float* h_bf_distances, float* h_pointset,
|
||||
float* h_query, int kth, int thelier, int nchunks, int pointdim,
|
||||
int signallength, unsigned int useMaxNorm) {
|
||||
float *d_bf_pointset, *d_bf_query;
|
||||
int *d_bf_indexes;
|
||||
float *d_bf_distances;
|
||||
|
||||
unsigned int meminputsignalquerypointset= pointdim * signallength * sizeof(float);
|
||||
unsigned int mem_bfcl_outputsignaldistances= kth * signallength * sizeof(float);
|
||||
unsigned int mem_bfcl_outputsignalindexes = kth * signallength * sizeof(int);
|
||||
|
||||
checkCudaErrors( cudaMalloc( (void**) &(d_bf_query), meminputsignalquerypointset));
|
||||
checkCudaErrors( cudaMalloc( (void**) &(d_bf_pointset), meminputsignalquerypointset));
|
||||
//GPU output
|
||||
checkCudaErrors( cudaMalloc( (void**) &(d_bf_distances), mem_bfcl_outputsignaldistances ));
|
||||
checkCudaErrors( cudaMalloc( (void**) &(d_bf_indexes), mem_bfcl_outputsignalindexes ));
|
||||
|
||||
cudaError_t error = cudaGetLastError();
|
||||
if(error!=cudaSuccess){
|
||||
fprintf(stderr,"%s",cudaGetErrorString(error));
|
||||
return 0;
|
||||
}
|
||||
|
||||
//Upload input data
|
||||
checkCudaErrors( cudaMemcpy(d_bf_query, h_query, meminputsignalquerypointset, cudaMemcpyHostToDevice ));
|
||||
checkCudaErrors( cudaMemcpy(d_bf_pointset, h_pointset, meminputsignalquerypointset, cudaMemcpyHostToDevice ));
|
||||
error = cudaGetLastError();
|
||||
if(error!=cudaSuccess){
|
||||
fprintf(stderr,"%s",cudaGetErrorString(error));
|
||||
return 0;
|
||||
}
|
||||
// Kernel parameters
|
||||
dim3 threads(1,1,1);
|
||||
dim3 grid(1,1,1);
|
||||
threads.x = 512;
|
||||
grid.x = (signallength-1)/threads.x + 1;
|
||||
int memkernel = kth*sizeof(float)*threads.x+\
|
||||
kth*sizeof(int)*threads.x;
|
||||
int triallength = signallength / nchunks;
|
||||
|
||||
// Pointer to the function used to calculate norms
|
||||
normFunction_t *normFunction;
|
||||
cudaMalloc( (void **) &normFunction, sizeof(normFunction_t) );
|
||||
if (useMaxNorm) {
|
||||
cudaMemcpyFromSymbol(normFunction, pMaxNorm, sizeof(normFunction_t));
|
||||
} else {
|
||||
cudaMemcpyFromSymbol(normFunction, pSquareNorm, sizeof(normFunction_t));
|
||||
}
|
||||
|
||||
// Launch kernel
|
||||
kernelKNNshared<<<grid.x, threads.x, memkernel>>>(d_bf_query, d_bf_pointset,
|
||||
d_bf_indexes, d_bf_distances, pointdim, triallength, signallength, kth,
|
||||
thelier, normFunction);
|
||||
|
||||
checkCudaErrors( cudaDeviceSynchronize() );
|
||||
|
||||
//Download result
|
||||
checkCudaErrors( cudaMemcpy( h_bf_distances, d_bf_distances, mem_bfcl_outputsignaldistances, cudaMemcpyDeviceToHost) );
|
||||
checkCudaErrors( cudaMemcpy( h_bf_indexes, d_bf_indexes, mem_bfcl_outputsignalindexes, cudaMemcpyDeviceToHost) );
|
||||
error = cudaGetLastError();
|
||||
if(error!=cudaSuccess){
|
||||
fprintf(stderr,"%s",cudaGetErrorString(error));
|
||||
return 0;
|
||||
}
|
||||
|
||||
//Free resources
|
||||
checkCudaErrors(cudaFree(d_bf_query));
|
||||
checkCudaErrors(cudaFree(d_bf_pointset));
|
||||
checkCudaErrors(cudaFree(d_bf_distances));
|
||||
checkCudaErrors(cudaFree(d_bf_indexes));
|
||||
cudaFree(normFunction);
|
||||
// cudaDeviceReset();
|
||||
if(error!=cudaSuccess){
|
||||
fprintf(stderr,"%s",cudaGetErrorString(error));
|
||||
return 0;
|
||||
}
|
||||
|
||||
return 1;
|
||||
}
|
||||
#ifdef __cplusplus
|
||||
}
|
||||
#endif
|
||||
|
||||
#ifdef __cplusplus
|
||||
extern "C" {
|
||||
#endif
|
||||
int cudaFindKnnSetGPU(int* h_bf_indexes, float* h_bf_distances,
|
||||
float* h_pointset, float* h_query, int kth, int thelier, int nchunks,
|
||||
int pointdim, int signallength, unsigned int useMaxNorm, int deviceid) {
|
||||
cudaSetDevice(deviceid);
|
||||
return cudaFindKnn(h_bf_indexes, h_bf_distances, h_pointset, h_query,
|
||||
kth, thelier, nchunks, pointdim, signallength, useMaxNorm);
|
||||
}
|
||||
#ifdef __cplusplus
|
||||
}
|
||||
#endif
|
||||
|
||||
/*
|
||||
* Range search being radius a vector of length number points in queryset/pointset
|
||||
*/
|
||||
#ifdef __cplusplus
|
||||
extern "C" {
|
||||
#endif
|
||||
int cudaFindRSAll(int* h_bf_npointsrange, float* h_pointset, float* h_query,
|
||||
float* h_vecradius, int thelier, int nchunks, int pointdim,
|
||||
int signallength, unsigned int useMaxNorm) {
|
||||
|
||||
float *d_bf_pointset, *d_bf_query, *d_bf_vecradius;
|
||||
int *d_bf_npointsrange;
|
||||
|
||||
unsigned int meminputsignalquerypointset= pointdim * signallength * sizeof(float);
|
||||
unsigned int mem_bfcl_outputsignalnpointsrange= signallength * sizeof(int);
|
||||
unsigned int mem_bfcl_inputvecradius = signallength * sizeof(float);
|
||||
|
||||
checkCudaErrors( cudaMalloc( (void**) &(d_bf_query), meminputsignalquerypointset));
|
||||
checkCudaErrors( cudaMalloc( (void**) &(d_bf_pointset), meminputsignalquerypointset));
|
||||
checkCudaErrors( cudaMalloc( (void**) &(d_bf_npointsrange), mem_bfcl_outputsignalnpointsrange ));
|
||||
checkCudaErrors( cudaMalloc( (void**) &(d_bf_vecradius), mem_bfcl_inputvecradius ));
|
||||
|
||||
cudaError_t error = cudaGetLastError();
|
||||
if(error!=cudaSuccess){
|
||||
fprintf(stderr,"%s",cudaGetErrorString(error));
|
||||
return 0;
|
||||
}
|
||||
//Upload input data
|
||||
checkCudaErrors( cudaMemcpy(d_bf_query, h_query, meminputsignalquerypointset, cudaMemcpyHostToDevice ));
|
||||
checkCudaErrors( cudaMemcpy(d_bf_pointset, h_pointset, meminputsignalquerypointset, cudaMemcpyHostToDevice ));
|
||||
checkCudaErrors( cudaMemcpy(d_bf_vecradius, h_vecradius, mem_bfcl_inputvecradius, cudaMemcpyHostToDevice ));
|
||||
|
||||
error = cudaGetLastError();
|
||||
if(error!=cudaSuccess){
|
||||
fprintf(stderr,"%s",cudaGetErrorString(error));
|
||||
return 0;
|
||||
}
|
||||
|
||||
|
||||
// Kernel parameters
|
||||
dim3 threads(1,1,1);
|
||||
dim3 grid(1,1,1);
|
||||
threads.x = 512;
|
||||
grid.x = (signallength-1)/threads.x + 1;
|
||||
int memkernel = sizeof(int)*threads.x;
|
||||
int triallength = signallength / nchunks;
|
||||
|
||||
// Pointer to the function used to calculate norms
|
||||
normFunction_t *normFunction;
|
||||
cudaMalloc( (void **) &normFunction, sizeof(normFunction_t) );
|
||||
if (useMaxNorm) {
|
||||
cudaMemcpyFromSymbol(normFunction, pMaxNorm, sizeof(normFunction_t));
|
||||
} else {
|
||||
cudaMemcpyFromSymbol(normFunction, pSquareNorm, sizeof(normFunction_t));
|
||||
}
|
||||
|
||||
// Launch kernel
|
||||
kernelBFRSAllshared<<< grid.x, threads.x, memkernel>>>(
|
||||
d_bf_query, d_bf_pointset, d_bf_npointsrange, pointdim,
|
||||
triallength, signallength, thelier, d_bf_vecradius, normFunction);
|
||||
|
||||
checkCudaErrors(cudaDeviceSynchronize());
|
||||
|
||||
checkCudaErrors( cudaMemcpy( h_bf_npointsrange, d_bf_npointsrange,mem_bfcl_outputsignalnpointsrange, cudaMemcpyDeviceToHost) );
|
||||
|
||||
|
||||
if(error!=cudaSuccess){
|
||||
fprintf(stderr,"%s",cudaGetErrorString(error));
|
||||
return 0;
|
||||
}
|
||||
|
||||
// Free resources
|
||||
checkCudaErrors(cudaFree(d_bf_query));
|
||||
checkCudaErrors(cudaFree(d_bf_pointset));
|
||||
checkCudaErrors(cudaFree(d_bf_npointsrange));
|
||||
checkCudaErrors(cudaFree(d_bf_vecradius));
|
||||
checkCudaErrors(cudaFree(normFunction));
|
||||
cudaDeviceReset();
|
||||
if(error!=cudaSuccess){
|
||||
fprintf(stderr,"%s",cudaGetErrorString(error));
|
||||
return 0;
|
||||
}
|
||||
|
||||
return 1;
|
||||
}
|
||||
#ifdef __cplusplus
|
||||
}
|
||||
#endif
|
||||
|
||||
#ifdef __cplusplus
|
||||
extern "C" {
|
||||
#endif
|
||||
int cudaFindRSAllSetGPU(int* h_bf_npointsrange, float* h_pointset,
|
||||
float* h_query, float* h_vecradius, int thelier, int nchunks,
|
||||
int pointdim, int signallength, unsigned int useMaxNorm, int deviceid) {
|
||||
cudaSetDevice(deviceid);
|
||||
return cudaFindRSAll(h_bf_npointsrange, h_pointset, h_query, h_vecradius,
|
||||
thelier, nchunks, pointdim, signallength, useMaxNorm);
|
||||
}
|
||||
#ifdef __cplusplus
|
||||
}
|
||||
#endif
|
||||
|
||||
|
||||
#ifdef __cplusplus
|
||||
extern "C" {
|
||||
#endif
|
||||
int findRadiiAlgorithm2(float *radii, const float *data, const int *indexes,
|
||||
unsigned int k, unsigned int dim, unsigned int N) {
|
||||
|
||||
unsigned int i, j;
|
||||
|
||||
for (j = 0; j < N; j++) {
|
||||
radii[j] = 0.0f;
|
||||
for (i = 0; i < k; i++) {
|
||||
float d = maxMetricPoints(data + j, data + indexes[j + i*N], dim, N);
|
||||
if (d > radii[j]) {
|
||||
radii[j] = d;
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
return 1;
|
||||
|
||||
}
|
||||
#ifdef __cplusplus
|
||||
}
|
||||
#endif
|
||||
|
||||
|
||||
|
||||
#ifdef __cplusplus
|
||||
extern "C" {
|
||||
#endif
|
||||
int computeSumDigammas(float *sumDiGammas, int *nx, int *ny, unsigned int N) {
|
||||
|
||||
int *d_nx, *d_ny;
|
||||
float *d_sumDiGammas;
|
||||
|
||||
checkCudaErrors( cudaMalloc((void **) &d_nx, N * sizeof(int)) );
|
||||
checkCudaErrors( cudaMalloc((void **) &d_ny, N * sizeof(int)) );
|
||||
checkCudaErrors( cudaMalloc((void **) &d_sumDiGammas, sizeof(float)) );
|
||||
|
||||
checkCudaErrors( cudaMemcpy(d_nx, nx, N*sizeof(int), cudaMemcpyHostToDevice) );
|
||||
checkCudaErrors( cudaMemcpy(d_ny, ny, N*sizeof(int), cudaMemcpyHostToDevice) );
|
||||
|
||||
unsigned int threads_per_block = 512;
|
||||
dim3 n_blocks, n_threads;
|
||||
n_blocks.x = ((N % threads_per_block) != 0) ? (N / threads_per_block + 1) : (N / threads_per_block);
|
||||
n_threads.x = threads_per_block;
|
||||
|
||||
printf("blocks = %i, threads = %i\n", n_blocks.x, n_threads.x);
|
||||
|
||||
// reduce6<<<n_blocks, n_threads>>>(d_nx, d_ny, d_sumDiGammas, N);
|
||||
reduce1<<<n_blocks, n_threads>>>(d_nx, d_ny, d_sumDiGammas, N);
|
||||
|
||||
checkCudaErrors( cudaDeviceSynchronize() );
|
||||
|
||||
checkCudaErrors( cudaMemcpy(sumDiGammas, d_sumDiGammas, sizeof(float), cudaMemcpyDeviceToHost) );
|
||||
|
||||
cudaFree(d_nx);
|
||||
cudaFree(d_ny);
|
||||
cudaFree(d_sumDiGammas);
|
||||
|
||||
*sumDiGammas = 0.0;
|
||||
|
||||
return 1;
|
||||
|
||||
}
|
||||
#ifdef __cplusplus
|
||||
}
|
||||
#endif
|
||||
|
||||
|
||||
|
|
@ -0,0 +1,26 @@
|
|||
#ifndef GPUKNNLIBRARY_H
|
||||
#define GPUKNNLIBRARY_H
|
||||
|
||||
int cudaFindKnn(int* h_bf_indexes, float* h_bf_distances, float* h_pointset,
|
||||
float* h_query, int kth, int thelier, int nchunks, int pointdim,
|
||||
int signallength, unsigned int useMaxNorm);
|
||||
|
||||
int cudaFindKnnSetGPU(int* h_bf_indexes, float* h_bf_distances,
|
||||
float* h_pointset, float* h_query, int kth, int thelier, int nchunks,
|
||||
int pointdim, int signallength, unsigned int useMaxNorm, int deviceid);
|
||||
|
||||
int cudaFindRSAll(int* h_bf_npointsrange, float* h_pointset,
|
||||
float* h_query, float* h_vecradius, int thelier, int nchunks,
|
||||
int pointdim, int signallength, unsigned int useMaxNorm);
|
||||
|
||||
int cudaFindRSAllSetGPU(int* h_bf_npointsrange, float* h_pointset,
|
||||
float* h_query, float* h_vecradius, int thelier, int nchunks,
|
||||
int pointdim, int signallength, unsigned int useMaxNorm, int deviceid);
|
||||
|
||||
int findRadiiAlgorithm2(float *radii, const float *data, const int *indexes,
|
||||
unsigned int k, unsigned int dim, unsigned int N);
|
||||
|
||||
int computeSumDigammas(float *sumDiGammas, int *nx, int *ny, unsigned int N);
|
||||
|
||||
#endif
|
||||
|
||||
|
|
@ -0,0 +1,230 @@
|
|||
////////////////////////////////////////////////////////////////////////////////
|
||||
// These are CUDA Helper functions
|
||||
|
||||
#include <stdlib.h>
|
||||
#include <stdio.h>
|
||||
|
||||
// This will output the proper CUDA error strings in the event that a CUDA host call returns an error
|
||||
#define checkCudaErrors(err) __checkCudaErrors (err, __FILE__, __LINE__)
|
||||
|
||||
inline void __checkCudaErrors(cudaError err, const char *file, const int line)
|
||||
{
|
||||
if (cudaSuccess != err)
|
||||
{
|
||||
fprintf(stderr, "%s(%i) : CUDA Runtime API error %d: %s.\n",file, line, (int)err, cudaGetErrorString(err));
|
||||
exit(-1);
|
||||
}
|
||||
}
|
||||
|
||||
// This will output the proper error string when calling cudaGetLastError
|
||||
#define getLastCudaError(msg) __getLastCudaError (msg, __FILE__, __LINE__)
|
||||
|
||||
inline void __getLastCudaError(const char *errorMessage, const char *file, const int line)
|
||||
{
|
||||
cudaError_t err = cudaGetLastError();
|
||||
|
||||
if (cudaSuccess != err)
|
||||
{
|
||||
fprintf(stderr, "%s(%i) : getLastCudaError() CUDA error : %s : (%d) %s.\n",
|
||||
file, line, errorMessage, (int)err, cudaGetErrorString(err));
|
||||
exit(-1);
|
||||
}
|
||||
}
|
||||
|
||||
// General GPU Device CUDA Initialization
|
||||
int gpuDeviceInit(int devID)
|
||||
{
|
||||
int deviceCount;
|
||||
checkCudaErrors(cudaGetDeviceCount(&deviceCount));
|
||||
|
||||
if (deviceCount == 0)
|
||||
{
|
||||
fprintf(stderr, "gpuDeviceInit() CUDA error: no devices supporting CUDA.\n");
|
||||
exit(-1);
|
||||
}
|
||||
|
||||
if (devID < 0)
|
||||
{
|
||||
devID = 0;
|
||||
}
|
||||
|
||||
if (devID > deviceCount-1)
|
||||
{
|
||||
fprintf(stderr, "\n");
|
||||
fprintf(stderr, ">> %d CUDA capable GPU device(s) detected. <<\n", deviceCount);
|
||||
fprintf(stderr, ">> gpuDeviceInit (-device=%d) is not a valid GPU device. <<\n", devID);
|
||||
fprintf(stderr, "\n");
|
||||
return -devID;
|
||||
}
|
||||
|
||||
cudaDeviceProp deviceProp;
|
||||
checkCudaErrors(cudaGetDeviceProperties(&deviceProp, devID));
|
||||
|
||||
if(!deviceProp.deviceOverlap){
|
||||
printf("\n Device will not handle overlaps, so no speed up from streams \n");
|
||||
}
|
||||
if (deviceProp.major < 1)
|
||||
{
|
||||
fprintf(stderr, "gpuDeviceInit(): GPU device does not support CUDA.\n");
|
||||
exit(-1);
|
||||
}
|
||||
|
||||
checkCudaErrors(cudaSetDevice(devID));
|
||||
printf("gpuDeviceInit() CUDA Device [%d]: \"%s\n", devID, deviceProp.name);
|
||||
|
||||
|
||||
return devID;
|
||||
}
|
||||
|
||||
#ifndef MAX
|
||||
#define MAX(a,b) (a > b ? a : b)
|
||||
#endif
|
||||
|
||||
// Beginning of GPU Architecture definitions
|
||||
inline int _ConvertSMVer2Cores(int major, int minor)
|
||||
{
|
||||
// Defines for GPU Architecture types (using the SM version to determine the # of cores per SM
|
||||
typedef struct
|
||||
{
|
||||
int SM; // 0xMm (hexidecimal notation), M = SM Major version, and m = SM minor version
|
||||
int Cores;
|
||||
} sSMtoCores;
|
||||
|
||||
sSMtoCores nGpuArchCoresPerSM[] =
|
||||
{
|
||||
{ 0x10, 8 }, // Tesla Generation (SM 1.0) G80 class
|
||||
{ 0x11, 8 }, // Tesla Generation (SM 1.1) G8x class
|
||||
{ 0x12, 8 }, // Tesla Generation (SM 1.2) G9x class
|
||||
{ 0x13, 8 }, // Tesla Generation (SM 1.3) GT200 class
|
||||
{ 0x20, 32 }, // Fermi Generation (SM 2.0) GF100 class
|
||||
{ 0x21, 48 }, // Fermi Generation (SM 2.1) GF10x class
|
||||
{ 0x30, 192}, // Fermi Generation (SM 3.0) GK10x class
|
||||
{ -1, -1 }
|
||||
};
|
||||
|
||||
int index = 0;
|
||||
|
||||
while (nGpuArchCoresPerSM[index].SM != -1)
|
||||
{
|
||||
if (nGpuArchCoresPerSM[index].SM == ((major << 4) + minor))
|
||||
{
|
||||
return nGpuArchCoresPerSM[index].Cores;
|
||||
}
|
||||
|
||||
index++;
|
||||
}
|
||||
|
||||
printf("MapSMtoCores undefined SM %d.%d is undefined (please update to the latest SDK)!\n", major, minor);
|
||||
return -1;
|
||||
}
|
||||
// end of GPU Architecture definitions
|
||||
|
||||
|
||||
// This function returns the best GPU (with maximum GFLOPS)
|
||||
int gpuGetMaxGflopsDeviceId()
|
||||
{
|
||||
int current_device = 0, sm_per_multiproc = 0;
|
||||
int max_compute_perf = 0, max_perf_device = 0;
|
||||
int device_count = 0, best_SM_arch = 0;
|
||||
cudaDeviceProp deviceProp;
|
||||
cudaGetDeviceCount(&device_count);
|
||||
|
||||
// Find the best major SM Architecture GPU device
|
||||
while (current_device < device_count)
|
||||
{
|
||||
cudaGetDeviceProperties(&deviceProp, current_device);
|
||||
|
||||
if (deviceProp.major > 0 && deviceProp.major < 9999)
|
||||
{
|
||||
best_SM_arch = MAX(best_SM_arch, deviceProp.major);
|
||||
}
|
||||
|
||||
current_device++;
|
||||
}
|
||||
|
||||
// Find the best CUDA capable GPU device
|
||||
current_device = 0;
|
||||
|
||||
while (current_device < device_count)
|
||||
{
|
||||
cudaGetDeviceProperties(&deviceProp, current_device);
|
||||
|
||||
if (deviceProp.major == 9999 && deviceProp.minor == 9999)
|
||||
{
|
||||
sm_per_multiproc = 1;
|
||||
}
|
||||
else
|
||||
{
|
||||
sm_per_multiproc = _ConvertSMVer2Cores(deviceProp.major, deviceProp.minor);
|
||||
}
|
||||
|
||||
int compute_perf = deviceProp.multiProcessorCount * sm_per_multiproc * deviceProp.clockRate;
|
||||
|
||||
if (compute_perf > max_compute_perf)
|
||||
{
|
||||
// If we find GPU with SM major > 2, search only these
|
||||
if (best_SM_arch > 2)
|
||||
{
|
||||
// If our device==dest_SM_arch, choose this, or else pass
|
||||
if (deviceProp.major == best_SM_arch)
|
||||
{
|
||||
max_compute_perf = compute_perf;
|
||||
max_perf_device = current_device;
|
||||
}
|
||||
}
|
||||
else
|
||||
{
|
||||
max_compute_perf = compute_perf;
|
||||
max_perf_device = current_device;
|
||||
}
|
||||
}
|
||||
|
||||
++current_device;
|
||||
}
|
||||
|
||||
return max_perf_device;
|
||||
}
|
||||
|
||||
|
||||
// Initialization code to find the best CUDA Device
|
||||
int findCudaDevice(int argc, const char **argv)
|
||||
{
|
||||
cudaDeviceProp deviceProp;
|
||||
int devID = 0;
|
||||
|
||||
// If the command-line has a device number specified, use it
|
||||
/*if (checkCmdLineFlag(argc, argv, "device"))
|
||||
{
|
||||
devID = getCmdLineArgumentInt(argc, argv, "device=");
|
||||
|
||||
if (devID < 0)
|
||||
{
|
||||
printf("Invalid command line parameter\n ");
|
||||
exit(-1);
|
||||
}
|
||||
else
|
||||
{
|
||||
devID = gpuDeviceInit(devID);
|
||||
|
||||
if (devID < 0)
|
||||
{
|
||||
printf("exiting...\n");
|
||||
//shrQAFinishExit(argc, (const char **)argv, QA_FAILED);
|
||||
exit(-1);
|
||||
}
|
||||
}
|
||||
}
|
||||
else
|
||||
{*/
|
||||
// Otherwise pick the device with highest Gflops/s
|
||||
devID = gpuGetMaxGflopsDeviceId();
|
||||
checkCudaErrors(cudaSetDevice(devID));
|
||||
checkCudaErrors(cudaGetDeviceProperties(&deviceProp, devID));
|
||||
printf("GPU Device %d: \"%s\" with compute capability %d.%d\n\n", devID, deviceProp.name, deviceProp.major, deviceProp.minor);
|
||||
//}
|
||||
|
||||
return devID;
|
||||
}
|
||||
// end of CUDA Helper Functions
|
||||
|
||||
|
||||
|
|
@ -0,0 +1,316 @@
|
|||
#include <jni.h>
|
||||
#include <stdlib.h>
|
||||
#include <stdio.h>
|
||||
|
||||
#include "digamma.c"
|
||||
#include "gpuKnnLibrary.h"
|
||||
|
||||
#define check(ans) { _check((ans), __FILE__, __LINE__); }
|
||||
|
||||
#ifdef __cplusplus
|
||||
extern "C" {
|
||||
#endif
|
||||
/*
|
||||
* Class: infodynamics_measures_continuous_kraskov_MutualInfoCalculatorMultiVariateKraskov
|
||||
* Method: MIKraskov
|
||||
* Signature: (I[DI[DIIZZZ)[D
|
||||
*/
|
||||
JNIEXPORT jdoubleArray JNICALL
|
||||
Java_infodynamics_measures_continuous_kraskov_MutualInfoCalculatorMultiVariateKraskov_MIKraskov(
|
||||
JNIEnv *env, jobject thisObj, jint j_N,
|
||||
jobjectArray j_sourceArray, jint j_dimx,
|
||||
jobjectArray j_destArray, jint j_dimy,
|
||||
jint j_k, jboolean j_returnLocals, jboolean j_useMaxNorm, jboolean j_isAlgorithm1) {
|
||||
|
||||
// Declare variables
|
||||
// Variables starting with j_ are Java types
|
||||
float *source, *dest, *pointset, *radii, *distances;
|
||||
int *nx, *ny, *indexes;
|
||||
jdouble *jd_source, *jd_dest, *result;
|
||||
unsigned int N, dimx, dimy, dims, k, returnLocals, useMaxNorm, isAlgorithm1, result_size;
|
||||
float sumNx, sumNy, sumDiGammas;
|
||||
unsigned int i, j;
|
||||
long double digammaK, digammaN;
|
||||
int thelier = 0, nchunks = 1, err;
|
||||
const int verbose = 0;
|
||||
|
||||
if (verbose) printf("Verbose is on.\n");
|
||||
|
||||
// Check that incoming data has correct size
|
||||
jsize sourceLength = (*env)->GetArrayLength(env, j_sourceArray);
|
||||
jsize destLength = (*env)->GetArrayLength(env, j_destArray);
|
||||
|
||||
if (sourceLength != j_N || destLength != j_N) {
|
||||
cudaDeviceReset();
|
||||
jclass Exception = (*env)->FindClass(env, "java/lang/Exception");
|
||||
(*env)->ThrowNew(env, Exception, "Data has wrong length.");
|
||||
}
|
||||
|
||||
// 1. Load data into GPU
|
||||
// =====================
|
||||
// 1a. Copy atomic variables from Java
|
||||
N = j_N;
|
||||
k = j_k;
|
||||
dimx = j_dimx;
|
||||
dimy = j_dimy;
|
||||
dims = dimx + dimy;
|
||||
returnLocals = j_returnLocals ? 1 : 0;
|
||||
useMaxNorm = j_useMaxNorm ? 1 : 0;
|
||||
isAlgorithm1 = j_isAlgorithm1 ? 1 : 0;
|
||||
|
||||
if (verbose > 0) {
|
||||
printf("Loading data from Java to C.\n");
|
||||
printf("N = %i, dimx = %i, dimy = %i\n", N, dimx, dimy);
|
||||
}
|
||||
|
||||
// 1b. Get arrays from Java
|
||||
source = (float *) malloc(N * dimx * sizeof(float));
|
||||
dest = (float *) malloc(N * dimy * sizeof(float));
|
||||
pointset = (float *) malloc(N * dims * sizeof(float));
|
||||
|
||||
if (NULL == source || NULL == dest || NULL == pointset) {
|
||||
printf("Error allocating data.\n");
|
||||
cudaDeviceReset();
|
||||
jclass Exception = (*env)->FindClass(env, "java/lang/Exception");
|
||||
(*env)->ThrowNew(env, Exception, "Error allocating data.");
|
||||
}
|
||||
|
||||
for (i = 0; i < N; i++) {
|
||||
jdoubleArray j_sourceRow = (*env)->GetObjectArrayElement(env, j_sourceArray, i);
|
||||
jdoubleArray j_destRow = (*env)->GetObjectArrayElement(env, j_destArray, i);
|
||||
|
||||
jdouble *sourceRow = (*env)->GetDoubleArrayElements(env, j_sourceRow, NULL);
|
||||
jdouble *destRow = (*env)->GetDoubleArrayElements(env, j_destRow, NULL);
|
||||
|
||||
// Data in java are doubles, but GPUs need floats.
|
||||
// We have to cast them manually, so we can't memcopy
|
||||
|
||||
// The following two for-loops get two matrices in T-by-D indexing (i.e.
|
||||
// first dimension is time, second is variable) and return the data in
|
||||
// column-major form
|
||||
for (j = 0; j < dimx; j++) {
|
||||
register float x = (float) sourceRow[j];
|
||||
source[N*j + i] = x;
|
||||
pointset[N*j + i] = x;
|
||||
}
|
||||
|
||||
for (j = 0; j < dimy; j++) {
|
||||
register float x = (float) destRow[j];
|
||||
dest[N*j + i] = x;
|
||||
pointset[N*dimx + N*j + i] = x;
|
||||
}
|
||||
|
||||
// // The following two for-loops get the same matrices but return arrays
|
||||
// // in row-major form. Code is left here for future development purposes.
|
||||
// for (j = 0; j < dimx; j++) {
|
||||
// register float x = (float) sourceRow[j];
|
||||
// source[dimx*i + j] = x;
|
||||
// pointset[dims*i + j] = x;
|
||||
// }
|
||||
|
||||
// for (j = 0; j < dimy; j++) {
|
||||
// register float x = (float) destRow[j];
|
||||
// dest[dimy*i + j] = x
|
||||
// pointset[i*dims + dimx + j] = x;
|
||||
// }
|
||||
|
||||
|
||||
(*env)->ReleaseDoubleArrayElements(env, j_sourceRow, sourceRow, 0);
|
||||
(*env)->ReleaseDoubleArrayElements(env, j_destRow, destRow, 0);
|
||||
|
||||
(*env)->DeleteLocalRef(env, j_sourceRow);
|
||||
(*env)->DeleteLocalRef(env, j_destRow);
|
||||
|
||||
// FIXME: I'm not entirely sure I'm freeing all the memory here. I should
|
||||
// check for memory leaks more carefully.
|
||||
|
||||
}
|
||||
|
||||
if (verbose > 1) {
|
||||
printf("Flat source array: ");
|
||||
for (i = 0; i < N*dimx; i++) printf("%f, ", source[i]);
|
||||
printf("\n");
|
||||
printf("Flat dest array: ");
|
||||
for (i = 0; i < N*dimy; i++) printf("%f, ", dest[i]);
|
||||
printf("\n");
|
||||
printf("Flat pointset: ");
|
||||
for (i = 0; i < N*dims; i++) printf("%f, ", pointset[i]);
|
||||
printf("\n");
|
||||
}
|
||||
|
||||
|
||||
// 2. Get indices and distances to the k NN of each point
|
||||
// ======================
|
||||
|
||||
indexes = (int *) malloc(N * k * sizeof(int));
|
||||
distances = (float *) malloc(N * k * sizeof(float));
|
||||
|
||||
err = cudaFindKnn(indexes, distances, pointset, pointset, k,
|
||||
thelier, nchunks, dims, N, useMaxNorm);
|
||||
if (err != 1) {
|
||||
cudaDeviceReset();
|
||||
jclass Exception = (*env)->FindClass(env, "java/lang/Exception");
|
||||
(*env)->ThrowNew(env, Exception, "Cuda error finding nearest neighbours.");
|
||||
}
|
||||
|
||||
if (verbose > 1) {
|
||||
printf("Distances: \n");
|
||||
for (j = 0; j < N; j++) {
|
||||
for (i = 0; i < k; i++) {
|
||||
printf("%f\t", distances[j + i*N]);
|
||||
}
|
||||
printf("\n");
|
||||
}
|
||||
printf("Indexes: \n");
|
||||
for (j = 0; j < N; j++) {
|
||||
for (i = 0; i < k; i++) {
|
||||
printf("%i\t", indexes[j + i*N]);
|
||||
}
|
||||
printf("\n");
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
// 3. Find distance R from each point to its k-th neighbour in joint space
|
||||
// ======================
|
||||
|
||||
if (verbose > 0) {
|
||||
printf("Getting radii from calculated dists.\n");
|
||||
}
|
||||
|
||||
radii = (float *) malloc(N * sizeof(float));
|
||||
for (i = 0; i < N; i++) {
|
||||
radii[i] = distances[N*(k-1) + i];
|
||||
}
|
||||
|
||||
if (verbose > 1) {
|
||||
for (i = 0; i < N; i++) {
|
||||
printf("%f\t", radii[i]);
|
||||
}
|
||||
printf("\n");
|
||||
}
|
||||
|
||||
|
||||
// 4. Count points strictly within R in the X-space
|
||||
// ======================
|
||||
|
||||
if (verbose > 0) {
|
||||
printf("Counting points within radius.\n");
|
||||
}
|
||||
|
||||
// If we're using algorithm 2 then we need the radius in the marginal space,
|
||||
// not the joint (as calculated above)
|
||||
if (!isAlgorithm1) {
|
||||
findRadiiAlgorithm2(radii, source, indexes, k, dimx, N);
|
||||
}
|
||||
|
||||
nx = (int *) malloc(N * sizeof(int));
|
||||
err = cudaFindRSAll(nx, source, source, radii, thelier, nchunks, dimx, N, useMaxNorm);
|
||||
if (err != 1) {
|
||||
cudaDeviceReset();
|
||||
jclass Exception = (*env)->FindClass(env, "java/lang/Exception");
|
||||
(*env)->ThrowNew(env, Exception, "Cuda error during source range search.");
|
||||
}
|
||||
|
||||
if (verbose > 0) {
|
||||
printf("X count: ");
|
||||
for (i = 0; i < N; i++) printf("%i, ", nx[i]);
|
||||
printf("\n\n");
|
||||
}
|
||||
|
||||
// 4. Count points strictly within R in the X-space
|
||||
// ======================
|
||||
|
||||
if (!isAlgorithm1) {
|
||||
findRadiiAlgorithm2(radii, dest, indexes, k, dimy, N);
|
||||
}
|
||||
|
||||
ny = (int *) malloc(N * sizeof(int));
|
||||
err = cudaFindRSAll(ny, dest, dest, radii, thelier, nchunks, dimy, N, useMaxNorm);
|
||||
if (err != 1) {
|
||||
cudaDeviceReset();
|
||||
jclass Exception = (*env)->FindClass(env, "java/lang/Exception");
|
||||
(*env)->ThrowNew(env, Exception, "Cuda error during dest range search.");
|
||||
}
|
||||
|
||||
if (verbose > 0) {
|
||||
printf("Y count: ");
|
||||
for (i = 0; i < N; i++) printf("%i, ", ny[i]);
|
||||
printf("\n\n");
|
||||
}
|
||||
|
||||
|
||||
// 6. Set locals or digammas for return
|
||||
// ======================
|
||||
if (verbose > 0) {
|
||||
printf("Start calculating digammas\n");
|
||||
}
|
||||
|
||||
if (returnLocals) {
|
||||
result = (jdouble *) malloc(N * sizeof(jdouble));
|
||||
result_size = N;
|
||||
digammaK = cpuDigamma(k);
|
||||
digammaN = cpuDigamma(N);
|
||||
for (i = 0; i < N; i++) {
|
||||
result[i] = digammaK + digammaN - cpuDigamma(nx[i] + 1) - cpuDigamma(ny[i] + 1);
|
||||
}
|
||||
|
||||
} else {
|
||||
|
||||
result = (jdouble *) malloc(3 * sizeof(jdouble));
|
||||
result_size = 3;
|
||||
|
||||
sumNx = 0;
|
||||
sumNy = 0;
|
||||
sumDiGammas = 0;
|
||||
for (i = 0; i < N; i++) {
|
||||
sumNx += nx[i];
|
||||
sumNy += ny[i];
|
||||
sumDiGammas += cpuDigamma(nx[i] + 1) + cpuDigamma(ny[i] + 1);
|
||||
|
||||
}
|
||||
|
||||
result[0] = (jdouble) sumDiGammas;
|
||||
result[1] = (jdouble) sumNx;
|
||||
result[2] = (jdouble) sumNy;
|
||||
|
||||
}
|
||||
|
||||
|
||||
// 7. Free memory and return
|
||||
// =========================
|
||||
|
||||
if (verbose > 0) {
|
||||
printf("Free memory\n");
|
||||
}
|
||||
|
||||
if (source) free(source);
|
||||
if (dest) free(dest);
|
||||
if (pointset) free(pointset);
|
||||
if (radii) free(radii);
|
||||
if (distances) free(distances);
|
||||
if (indexes) free(indexes);
|
||||
if (nx) free(nx);
|
||||
if (ny) free(ny);
|
||||
|
||||
// Set Java array for return
|
||||
jdoubleArray outJNIArray = (*env)->NewDoubleArray(env, result_size); // allocate
|
||||
if (NULL == outJNIArray) {
|
||||
cudaDeviceReset();
|
||||
jclass Exception = (*env)->FindClass(env, "java/lang/Exception");
|
||||
(*env)->ThrowNew(env, Exception, "Error creating return array.");
|
||||
}
|
||||
(*env)->SetDoubleArrayRegion(env, outJNIArray, 0 , result_size, result); // copy
|
||||
|
||||
if (result) free(result);
|
||||
|
||||
return outJNIArray;
|
||||
|
||||
|
||||
} // End of function MIKraskov1
|
||||
|
||||
|
||||
#ifdef __cplusplus
|
||||
}
|
||||
#endif
|
||||
|
||||
Loading…
Reference in New Issue