Операции со сложными матрицами с использованием jCUDA
Как лучше всего работать с комплексными числами, используя jCuda? Должен ли я использовать формат cuComplex или есть какое-то другое решение (например, массив с реальными и мнимыми частями, идущими одна за другой)? Я был бы очень признателен за примеры кода Java с такими вычислениями.
Поскольку моя цель состоит в том, чтобы решать большие системы линейных уравнений с комплексными числами с использованием графического процессора, я не хотел бы присоединяться только к jCuda. Каковы альтернативные способы проведения таких расчетов с помощью графического процессора?
1 ответ
Во-первых, касательно вашего вопроса о вычислениях с Java на GPU в целом, я написал несколько слов об этом здесь.
Ваш случай применения выглядит очень конкретным. Вероятно, вам следует более подробно описать, каково ваше истинное намерение, поскольку оно будет определять все проектные решения. До сих пор я могу дать только некоторые основные советы. Решение, какое решение является наиболее подходящим, зависит от вас.
Одна из главных трудностей при соединении между миром Java и миром графических процессоров - это принципиально иная обработка памяти.
Расположение памяти в C/C++
cuComplex Структура в CUDA определяется как
typedef float2 cuFloatComplex
typedef cuFloatComplex cuComplex;
где float2 в основном что-то вроде
struct float2 {
float x;
float y;
};
(с некоторыми дополнительными спецификаторами для выравнивания и т. д.)
Теперь, когда вы выделяете массив cuComplex значения в программе C/C++, тогда вы просто напишите что-то вроде
cuComplex *c = new cuComplex[100];
В этом случае гарантируется, что память обо всех этих cuComplex значения будут единым, непрерывным блоком памяти. Этот блок памяти просто состоит из всех x а также y Значения комплексных чисел, одно за другим:
_____________________________
c -> | x0 | y0 | x1 | y1 | x2 | y2 |...
|____|____|____|____|____|____|
Этот непрерывный блок памяти можно легко скопировать на устройство: один берет указатель и вызывает вызов
cudaMemcpy(device, c, sizeof(cuComplex)*n, cudaMemcpyHostToDevice);
Расположение памяти в Java
Рассмотрим случай, когда вы создаете класс Java, который структурно равен cuComplex структурировать и выделить массив из них:
class cuComplex {
public float x;
public float y;
}
cuComplex c[] = new cuComplex[100];
Тогда у вас нет ни одного непрерывного блока памяти float ценности. Вместо этого у вас есть массив ссылок на cuComplex объекты и соответствующие x а также y значения разбросаны по всему:
____________________
c -> | c0 | c1 | c2 |...
|______|______|______|
| | |
v v v
[x,y] [x,y] [x,y]
Ключевым моментом здесь является:
Вы не можете скопировать массив (Java) cuComplex объекты на устройство!,
Это имеет несколько последствий. В комментариях вы уже ссылались на cublasSetVector метод, который принимает cuComplex массив в качестве аргумента, и где я попытался подчеркнуть, что это не самое эффективное решение, но только для удобства. На самом деле, этот метод работает путем внутреннего создания нового ByteBuffer чтобы иметь непрерывный блок памяти, заполняя это ByteBuffer со значениями из cuComplex[] массив, а затем скопировать это ByteBuffer к устройству.
Конечно, это накладывает накладные расходы, которых вы, скорее всего, хотели бы избежать в приложениях, критичных к производительности.
Есть несколько вариантов решения этой проблемы. К счастью, для комплексных чисел решение сравнительно просто:
Не используйте cuComplex структура для представления массивов комплексных чисел
Вместо этого вы должны представить свой массив комплексных чисел в виде единого непрерывного блока памяти, где чередуются действительные и мнимые части комплексных чисел, каждая из которых представляет собой единое целое. float или же double стоимость соответственно. Это позволит обеспечить наибольшую совместимость между различными бэкэндами (исключая некоторые детали, такие как требования к выравниванию).
К сожалению, это может вызвать некоторые неудобства и вызвать некоторые вопросы, и для этого не существует универсального решения.
Если попытаться обобщить это, чтобы ссылаться не только на комплексные числа, но и на "структуры" в целом, то можно применить "шаблон": можно создать интерфейсы для структур и создать коллекцию этих структур. это список экземпляров классов, реализующих этот интерфейс, которые все поддерживаются непрерывным блоком памяти. Это может быть подходящим для определенных случаев. Но для комплексных чисел издержки памяти, связанные с наличием Java-объекта для каждого комплексного числа, могут быть чрезмерно большими.
Другая крайность, только обработка сырья float[] или же double[] Массивы также могут быть не лучшим решением. Например: если у вас есть массив float Значения, которые представляют комплексные числа, тогда как вы могли бы умножить одно из этих комплексных чисел на другое?
Одно "промежуточное" решение могло бы создать интерфейс, позволяющий получать доступ к действительным и мнимым частям комплексных чисел. В реализации эти комплексные числа хранятся в одном массиве, как описано выше.
Я набросал такую реализацию здесь.
НОТА:
Это предназначено только в качестве примера, чтобы показать основную идею и показать, как она может работать вместе с чем-то вроде JCublas. Для вас может быть более подходящей другая стратегия, в зависимости от ваших реальных целей: какие еще бэкэнды (кроме JCuda) должны быть? Насколько "удобно" обрабатывать комплексные числа на стороне Java? Как должны выглядеть структуры (классы / интерфейсы) для обработки комплексных чисел на стороне Java?
Если коротко: у вас должно быть четкое представление о том, что должно делать ваше приложение / библиотека, прежде чем приступить к реализации.
import static jcuda.jcublas.JCublas2.*;
import static jcuda.jcublas.cublasOperation.CUBLAS_OP_N;
import static jcuda.runtime.JCuda.*;
import java.util.Random;
import jcuda.*;
import jcuda.jcublas.cublasHandle;
import jcuda.runtime.cudaMemcpyKind;
// An interface describing an array of complex numbers, residing
// on the host, with methods for accessing the real and imaginary
// parts of the complex numbers, as well as methods for copying
// the underlying data from and to the device
interface cuComplexHostArray
{
int size();
float getReal(int i);
float getImag(int i);
void setReal(int i, float real);
void setImag(int i, float imag);
void set(int i, cuComplex c);
void set(int i, float real, float imag);
cuComplex get(int i, cuComplex c);
void copyToDevice(Pointer devicePointer);
void copyFromDevice(Pointer devicePointer);
}
// A default implementation of a cuComplexHostArray, backed
// by a single float[] array
class DefaultCuComplexHostArray implements cuComplexHostArray
{
private final int size;
private final float data[];
DefaultCuComplexHostArray(int size)
{
this.size = size;
this.data = new float[size * 2];
}
@Override
public int size()
{
return size;
}
@Override
public float getReal(int i)
{
return data[i+i];
}
@Override
public float getImag(int i)
{
return data[i+i+1];
}
@Override
public void setReal(int i, float real)
{
data[i+i] = real;
}
@Override
public void setImag(int i, float imag)
{
data[i+i+1] = imag;
}
@Override
public void set(int i, cuComplex c)
{
data[i+i+0] = c.x;
data[i+i+1] = c.y;
}
@Override
public void set(int i, float real, float imag)
{
data[i+i+0] = real;
data[i+i+1] = imag;
}
@Override
public cuComplex get(int i, cuComplex c)
{
float real = getReal(i);
float imag = getImag(i);
if (c != null)
{
c.x = real;
c.y = imag;
return c;
}
return cuComplex.cuCmplx(real, imag);
}
@Override
public void copyToDevice(Pointer devicePointer)
{
cudaMemcpy(devicePointer, Pointer.to(data),
size * Sizeof.FLOAT * 2,
cudaMemcpyKind.cudaMemcpyHostToDevice);
}
@Override
public void copyFromDevice(Pointer devicePointer)
{
cudaMemcpy(Pointer.to(data), devicePointer,
size * Sizeof.FLOAT * 2,
cudaMemcpyKind.cudaMemcpyDeviceToHost);
}
}
// An example that performs a "gemm" with complex numbers, once
// in Java and once in JCublas2, and verifies the result
public class JCublas2ComplexSample
{
public static void main(String args[])
{
testCgemm(500);
}
public static void testCgemm(int n)
{
cuComplex alpha = cuComplex.cuCmplx(0.3f, 0.2f);
cuComplex beta = cuComplex.cuCmplx(0.1f, 0.7f);
int nn = n * n;
System.out.println("Creating input data...");
Random random = new Random(0);
cuComplex[] rhA = createRandomComplexRawArray(nn, random);
cuComplex[] rhB = createRandomComplexRawArray(nn, random);
cuComplex[] rhC = createRandomComplexRawArray(nn, random);
random = new Random(0);
cuComplexHostArray hA = createRandomComplexHostArray(nn, random);
cuComplexHostArray hB = createRandomComplexHostArray(nn, random);
cuComplexHostArray hC = createRandomComplexHostArray(nn, random);
System.out.println("Performing Cgemm with Java...");
cgemmJava(n, alpha, rhA, rhB, beta, rhC);
System.out.println("Performing Cgemm with JCublas...");
cgemmJCublas(n, alpha, hA, hB, beta, hC);
boolean passed = isCorrectResult(hC, rhC);
System.out.println("testCgemm "+(passed?"PASSED":"FAILED"));
}
private static void cgemmJCublas(
int n,
cuComplex alpha,
cuComplexHostArray A,
cuComplexHostArray B,
cuComplex beta,
cuComplexHostArray C)
{
int nn = n * n;
// Create a CUBLAS handle
cublasHandle handle = new cublasHandle();
cublasCreate(handle);
// Allocate memory on the device
Pointer dA = new Pointer();
Pointer dB = new Pointer();
Pointer dC = new Pointer();
cudaMalloc(dA, nn * Sizeof.FLOAT * 2);
cudaMalloc(dB, nn * Sizeof.FLOAT * 2);
cudaMalloc(dC, nn * Sizeof.FLOAT * 2);
// Copy the memory from the host to the device
A.copyToDevice(dA);
B.copyToDevice(dB);
C.copyToDevice(dC);
// Execute cgemm
Pointer pAlpha = Pointer.to(new float[]{alpha.x, alpha.y});
Pointer pBeta = Pointer.to(new float[]{beta.x, beta.y});
cublasCgemm(handle, CUBLAS_OP_N, CUBLAS_OP_N, n, n, n,
pAlpha, dA, n, dB, n, pBeta, dC, n);
// Copy the result from the device to the host
C.copyFromDevice(dC);
// Clean up
cudaFree(dA);
cudaFree(dB);
cudaFree(dC);
cublasDestroy(handle);
}
private static void cgemmJava(
int n,
cuComplex alpha,
cuComplex A[],
cuComplex B[],
cuComplex beta,
cuComplex C[])
{
for (int i = 0; i < n; ++i)
{
for (int j = 0; j < n; ++j)
{
cuComplex prod = cuComplex.cuCmplx(0, 0);
for (int k = 0; k < n; ++k)
{
cuComplex ab =
cuComplex.cuCmul(A[k * n + i], B[j * n + k]);
prod = cuComplex.cuCadd(prod, ab);
}
cuComplex ap = cuComplex.cuCmul(alpha, prod);
cuComplex bc = cuComplex.cuCmul(beta, C[j * n + i]);
C[j * n + i] = cuComplex.cuCadd(ap, bc);
}
}
}
private static cuComplex[] createRandomComplexRawArray(
int n, Random random)
{
cuComplex c[] = new cuComplex[n];
for (int i = 0; i < n; i++)
{
float real = random.nextFloat();
float imag = random.nextFloat();
c[i] = cuComplex.cuCmplx(real, imag);
}
return c;
}
private static cuComplexHostArray createRandomComplexHostArray(
int n, Random random)
{
cuComplexHostArray c = new DefaultCuComplexHostArray(n);
for (int i = 0; i < n; i++)
{
float real = random.nextFloat();
float imag = random.nextFloat();
c.setReal(i, real);
c.setImag(i, imag);
}
return c;
}
private static boolean isCorrectResult(
cuComplexHostArray result, cuComplex reference[])
{
float errorNormX = 0;
float errorNormY = 0;
float refNormX = 0;
float refNormY = 0;
for (int i = 0; i < result.size(); i++)
{
float diffX = reference[i].x - result.getReal(i);
float diffY = reference[i].y - result.getImag(i);
errorNormX += diffX * diffX;
errorNormY += diffY * diffY;
refNormX += reference[i].x * result.getReal(i);
refNormY += reference[i].y * result.getImag(i);
}
errorNormX = (float) Math.sqrt(errorNormX);
errorNormY = (float) Math.sqrt(errorNormY);
refNormX = (float) Math.sqrt(refNormX);
refNormY = (float) Math.sqrt(refNormY);
if (Math.abs(refNormX) < 1e-6)
{
return false;
}
if (Math.abs(refNormY) < 1e-6)
{
return false;
}
return
(errorNormX / refNormX < 1e-6f) &&
(errorNormY / refNormY < 1e-6f);
}
}
(Кстати: я, вероятно, возьму части этого ответа и расширю их, чтобы они стали образцами и / или страницами "Как..." для JCuda. Задача предоставления некоторой информации, подобной этой, уже была на моем "задании" список на некоторое время).