4.5 Задача Дирихле для уравнения Пуассона

Для демонстрации решения задач дифференциальных уравнений с частными производными разберем простое одномерное уравнение Пуассона:

∂2u
∂x2-= f,x ∈ [0,l].
Оно с помощью использования конечных разностей в половинных узлах легко принимает следующий вид СЛАУ:
ui+1---2ui +-ui-1-= fi,i = 0,N-- 1,
      h2
или в матричном виде:
Ax  = b,
где A – трех-диагональная матрица со значениями на средней диагонали 2, на соседних -1.

Выберем метод сопряженных градиентов для решения данной системы, который гарантированно сходится за N итераций. Алгоритм:

1.
Вычислить r(0) = p(0) = b - Ax(0) для начального приближения x(0).
2.
Цикл do.
3.
α =   i-1 i-1
((rpi-1,,Arpi-)1).
4.
xi = xi-1 + αpi-1.
5.
xi = xi-1 + αpi-1.
6.
temp = ri-1 + αApi-1.
7.
ϵ = ∘ ----------
  (ri- 1,ri-1).
8.
β = ((termi-p1,t,reim-p1)).
9.
pi = temp + β temp.
10.
ri = temp.
11.
while ϵ > ε

Итак, на основе метода, посмотрим какие операции надо распараллелить: сложение векторов, скалярное умножение векторов, умножение скаляра на вектор и умножение матрицы на вектор. Идея состоит в том, чтобы ускорить вычисление каждой итерации путем ускорения наиболее тяжелой операции. В случае если распараллелить все операции скорость возрастет не значительно, на больших сетках и даже совсем на небольшое значение будет хуже на малых сетках. Так как Вы из предыдущих глав уже научились распараллеливать сложение и скалярное умножение векторов, распараллелить на основе данного кода Вам не составит труда.

Так выглядит операция умножения матрицы на вектор в последовательном варианте:

for (int i = 0; i < NSIZE; i++) { 
    for (int j = 0; j < NSIZE; j++) { 
        Y[i] += A[i,j] * X[j]; 
    } 
}

Так выглядит данная операция на Cuda:

// вычисление = AX 
__global__ void matdotvec(float *A, float *X, float *Y) { 
    __shared__ float Xds[threads]; // используем разделяемую память 
    intbx = blockIdx.x;    //индекс блока; 
    inttx = threadIdx.x;   //индекс нити в блоке; 
    // вычисление индекса 
    int Row = bx * blockDim.x + tx; 
    floatPvalue = 0; 
    for (unsigned int m = 0; m < (NSIZE-1)/threads+1; m++) { 
        if (m*threads + tx< NSIZE) 
        //переназначение вектора на нужный диапазон 
             Xds[tx] = X[m*threads + tx]; 
        else 
             Xds[tx] = 0; 
    __syncthreads();// точка синхронизации 
    for (unsigned int k = 0; k < threads; k++) 
        if((Row<NSIZE)&&(m*threads + k < NSIZE)) 
            //собирание умножения матрицы на вектор внутри блока 
            Pvalue += A[m*threads+Row*NSIZE+k] * Xds[k]; 
    __syncthreads();// точка синхронизации 
    } 
    if(Row < NSIZE) { 
        Y[Row] = Pvalue; // запись в результирующий вектор 
    } 
    __syncthreads();// точка синхронизации 
}

Итак, первые строки нам уже хорошо известны. Тут производится назначение индексов. Также в первых строках находиться выделение разделяемой памяти она используется внутри блока. Именно в неё будет записываться актуальный индекс для умножения куска строки матрицы на кусок вектора. Данное умножение записывается в переменную Pvalue, которая действительна для данного умножения. Остается лишь записать в нужную позицию результирующего вектора полученное значение Pvalue.

Полный код программы:

#include<iostream> 
#include <stdio.h> 
#include <math.h> 
#include <stdlib.h> 
 
#define imin(a,b) (a<b?a:b) // выбирает минимальное из a и b; 
 
const int NSIZE = 2000; // число неизвестных 
const int threads = 512; // кол-во нитей в блоке 
const int blocks = imin(32, ((NSIZE+threads - 1) / threads)); 
// кол-во блоков 
 
float *matrA; 
float *vecX, *vecP, *vecR, *vecB, *vecTemp, *vecAP; 
float alpha, beta, normaR; 
float *dev_vec1, *dev_vec2, *dev_matrA; 
// объявление переменных 
 
void initvector() // инициализация векторов 
{ 
    vecTemp = (float *)malloc(NSIZE * sizeof(float)); 
    vecX = (float *) malloc(NSIZE * sizeof(float)); 
    vecP = (float *) malloc(NSIZE * sizeof(float)); 
    vecR = (float *) malloc(NSIZE * sizeof(float)); 
    vecB = (float *) malloc(NSIZE * sizeof(float)); 
    vecAP = (float *) malloc(NSIZE * sizeof(float)); 
    for (int i = 0; i < NSIZE; i++) {// задаем значения 
        vecX[i] = vecP[i] = vecR[i] = vecTemp[i] = 0.0; 
        vecB[i] = 1.0 / (NSIZE * NSIZE); 
    } 
    // выделение памяти для векторов Cuda 
    cudaMalloc( (void**)&dev_vec1, NSIZE*sizeof(float) ) ; 
    cudaMalloc( (void**)&dev_vec2, NSIZE*sizeof(float) ) ; 
} 
 
void initmatrix() { // задание трехдиагональной матрицы 
    matrA = (float *)malloc(NSIZE * NSIZE * sizeof(float)); 
    for (int i = 0; i < NSIZE; i++) 
        for (int j = 0; j < NSIZE; j++) { 
            int k = i * NSIZE + j; 
            matrA[k] = 0; 
        } 
    for (int i = 0; i < NSIZE; i++) { 
        matrA[i * NSIZE + i] = 2.0; 
        if (i > 0) { 
            matrA[(i - 1)*NSIZE + i] = -1.0; 
        } 
        if (i < NSIZE - 1) { 
            matrA[(i + 1)*NSIZE + i] = -1.0; 
        } 
    } 
    // выделяем память и копируем для матрицы на видеокарте 
    cudaMalloc( (void**)&dev_matrA, NSIZE * NSIZE * sizeof(float)); 
    cudaMemcpy(dev_matrA, matrA,  NSIZE * NSIZE * sizeof(float), cudaMemcpyHostToDevice); 
} 
 
// скалярное умножение 
float vectorbyvector (float *vect1, float *vect2) 
{ 
    float result; 
    result = 0.0; 
    for (int i = 0; i < NSIZE; i++) { 
        result += vect1[i] * vect2[i]; 
    } 
    return result; 
} 
 
// сумма векторов 
float *vectorsum(float *vect1, float *vect2, bool flag) { 
    float *result; 
    result = (float *) malloc(NSIZE * sizeof(float)); 
    for (int i = 0; i < NSIZE; i++) { 
        result[i] = 0.0; 
    } 
    if (flag) { 
        for (int i = 0; i < NSIZE; i++) { 
            result[i] = vect1[i] + vect2[i]; 
        } 
    } 
    else { 
        for (int i = 0; i < NSIZE; i++) { 
            result[i] = vect1[i] - vect2[i]; 
        } 
    } 
    return result; 
} 
 
// умножение вектора на число 
float *vectorbynumber(float number, float *vect) { 
    float *result; 
    result = (float *) malloc(NSIZE * sizeof(float)); 
    for (int i = 0; i < NSIZE; i++) { 
        result[i] = 0.0; 
    } 
    for (int i = 0; i < NSIZE; i++) { 
        result[i] += number * vect[i]; 
    } 
    return result; 
} 
__global__ void matdotvec(float *A, float *X, float *Y) { 
    __shared__ float Xds[threads]; // используем разделяемую память 
    intbx = blockIdx.x;    //индекс блока; 
    inttx = threadIdx.x;   //индекс нити в блоке; 
    // вычисление индекса 
    int Row = bx * blockDim.x + tx; 
    floatPvalue = 0; 
    for (unsigned int m = 0; m < (NSIZE-1)/threads+1; m++) { 
        if (m*threads + tx< NSIZE) 
        //переназначение вектора на нужный диапазон 
             Xds[tx] = X[m*threads + tx]; 
        else 
             Xds[tx] = 0; 
    __syncthreads();// точка синхронизации 
    for (unsigned int k = 0; k < threads; k++) 
        if((Row<NSIZE)&&(m*threads + k < NSIZE)) 
            //собирание умножения матрицы на вектор внутри блока 
            Pvalue += A[m*threads+Row*NSIZE+k] * Xds[k]; 
    __syncthreads();// точка синхронизации 
    } 
    if(Row < NSIZE) { 
        Y[Row] = Pvalue; // запись в результирующий вектор 
    } 
    __syncthreads();// точка синхронизации 
} 
int main(void) { // главная часть программы 
    initvector();// инициализация 
    initmatrix();// инициализация 
    // копирование вектора X на видеокарту 
    cudaMemcpy(dev_vec1, vecX, NSIZE * sizeof(float), cudaMemcpyHostToDevice); 
    // Вычисление AX 
    matdotvec<<<blocks,threads>>>(dev_matrA, dev_vec1, dev_vec2); 
    // копирование обратно на хост 
    cudaMemcpy(vecAP, dev_vec2, NSIZE * sizeof(float), cudaMemcpyDeviceToHost); 
    vecR = vecP = vectorsum(vecB, vecAP, false); 
 
    intiter = 0; 
    do { // реализация итераций в методе сопряженных градиентов 
        iter++; 
        // копирование вектора X на видеокарту 
        cudaMemcpy(dev_vec1, vecP, NSIZE * sizeof(float), cudaMemcpyHostToDevice); 
        // Вычисление AP 
        matdotvec<<<blocks,threads>>>(dev_matrA, dev_vec1, dev_vec2); 
        // копирование обратно на хост 
        cudaMemcpy(vecAP, dev_vec2, NSIZE * sizeof(float), cudaMemcpyDeviceToHost); 
        alpha = vectorbyvector(vecR, vecR) / vectorbyvector (vecP, vecAP); 
        vecX = vectorsum(vecX, vectorbynumber(alpha, vecP), true); 
        vecTemp = vectorsum(vecR, vectorbynumber(alpha, vecAP), false); 
        normaR = sqrt(vectorbyvector(vecR, vecR)); 
        beta = vectorbyvector(vecTemp, vecTemp) / vectorbyvector(vecR, vecR); 
        vecP = vectorsum(vecTemp, vectorbynumber(beta, vecP), true); 
        vecR = vecTemp; 
    } while (normaR> 0.00001);// пока не сойдется по норме 
    FILE *fp; // запись в файл 
    fp = fopen("output", "w"); 
    for (int i = 0; i < NSIZE; i++) { 
        fprintf(fp, "%f \n", vecX[i]); 
    } 
    fclose(fp); 
    cudaFree(dev_vec1); // освобождение переменной cuda 
    cudaFree(dev_vec2); // освобождение переменной cuda 
    cudaFree(dev_matrA);// освобождение переменной cuda 
}