Матрица инверсии в OpenCL

Question

Я пытаюсь ускорить некоторые вычисления, используя OpenCL, и часть алгоритма состоит из инвертирования матрицы. Существует ли какая-либо библиотека с открытым исходным кодом или свободно доступный код для вычисления факторизации lu (lapack dgetrf и dgetri) матрицы или общей инверсии, написанной на OpenCL или CUDA? Матрица является реальной и квадратной, но не имеет никаких других специальных свойств, кроме этого. До сих пор мне удалось найти только основные реализации матрично-векторных операций blas в gpu.

Матрица довольно маленькая, всего около 60-100 строк и столбцов, поэтому она может быть вычислена быстрее на процессоре, но она используется как бы посередине алгоритма, поэтому мне придется перенести ее на хост, вычислить инвертировать, а затем передать результат обратно на устройство, где он затем используется в гораздо больших вычислениях.

Answer 1

Посмотрите на ВенуCL: http://viennacl.sourceforge.net/

Answer 2

У меня нет реализации в Open CL, но у обоих "Числовые рецепты" и у Джила Странга "Into to Applied Math" есть замечательные объяснения, которые было бы легко кодировать , «NR» имеет C-код, который вы можете адаптировать.

рассчитать обратное

Это неверно. Вы не вычисляете обратное с разложением LU, вы разлагаете матрицу. Если вам нужно обратное, вам придется выполнить прямую замену назад на серию векторов единиц. Это небольшая, но важная разница.

Answer 3

Я делаю исчисление до 2k x 2k в CPU, используя многопоточность, используя eigen lib, поэтому теперь это в 3,5-3,65 раза быстрее (зависит от размера матрицы), чем при использовании одного потока Я использовал процессор Intel Xeon 3,5 ГГц E5-1620 v3 и 16 ГБ оперативной памяти. (К сожалению, я удалил старую версию, чтобы добавить точные значения, но если у меня есть какой-то приоритет, я мог бы переписать sw)

Это мой матричный обратный алгоритм, с которым я сравнивал. (Это правильно, как я уже говорил, делая много тестов снова превосходный результат):

/*
Uses 2D arrays.
Main routines are:
init_2Dvector() that initializes any 2d vector (can be uchar, char, int, float or double)
multiply_2Dvector()
inverse()
*/

#include<iostream>
#include <vector>
#include <stdlib.h>
#include <time.h>

using namespace std;

/*
void print_2Dvector(vector<vector<double> >& vec)
{
    size_t xmax, ymax;
    ymax = vec.size();
    xmax = vec[0].size(); 

    int x, y;
    for (y = 0; y < ymax; y++)
    {
        for (x = 0; x < xmax; x++)
            cout << vec[y][x] << " \t";
        cout << endl;
    }
}*/

void print_2Dvector(vector<vector<double> >& vec,char *format="%lg \t")
{
    size_t xmax, ymax;
    ymax = vec.size();
    xmax = vec[0].size();

    int x, y;
    for (y = 0; y < ymax; y++)
    {
        {
            for (x = 0; x < xmax; x++)
                printf(format, vec[y][x]);
        }
        cout << endl;
    }
}

//Resizes to y_dim,x_dim any kind of 2d array:
template<typename T>
void init_2Dvector(vector<vector<T> >& vec, size_t y_dim, size_t x_dim)
{
    vec.resize(y_dim);
    for (size_t i = 0; i < vec.size(); i++)
        vec[i].resize(x_dim);
}
//Returns vec1*vec2. vec1 and 2 are not touch
vector< vector<double> > multiply_2Dvector(vector< vector<double> > vec1, vector< vector<double> > vec2)
{
    size_t xmax, ymax;
    ymax = vec1.size();   
    xmax = vec1[0].size();
    vector< vector<double> > vec3;
    if ((ymax != vec2[0].size()) || (xmax != vec2.size()))
    {
        cout << "ERROR on dim2_multiply() dimensions of vec2 not corresponding with vec1 ones" << endl; getchar(); return{};//returns a null
    }
    init_2Dvector(vec3, ymax, ymax);
    cout << "dimensions of vec3=" << vec3.size() << " x " << vec3[0].size() << endl;
    double xx;
    for (int y = 0; y < ymax; y++)
        for (int x = 0; x < ymax; x++)
        {
            xx = 0.0;
            for (int t = 0; t < xmax; t++)
                xx += vec1[y][t] * vec2[t][x];
            vec3[y][x] = xx;
        }
    return vec3;//ok
}

//returns inverse of x2, x2 is not modified
vector< vector<double> > inverse(vector< vector<double> > x)
{
    if (x.size() != x[0].size())
    {
        cout << "ERROR on inverse() not square array" << endl; getchar(); return{};//returns a null
    }

    size_t dim = x.size();
    int i, j, ord;
    vector< vector<double> > y(dim,vector<double>(dim));//output
    //init_2Dvector(y, dim, dim);
    //1. Unity array y: 
    for (i = 0; i < dim; i++)
    {
        y[i][i] = 1.0;
        for (j = i+1; j < dim; j++)
        {
            y[i][j]= y[j][i] = 0.0;
        }
    }

    double diagon, coef;
    double *ptrx, *ptry, *ptrx2, *ptry2;
    for (ord = 0; ord<dim; ord++)
    {
        //2 Hacemos diagonal de x =1
        int i2;
        if (fabs(x[ord][ord])<1e-15) //Si el elemento diagonal es 0 sumamos una columna que no sea 0 el elemento correspondiente
        {
            for (i2 = ord + 1; i2<dim; i2++)
            {
                if (fabs(x[i2][ord])>1e-15) break;
            }
            if (i2 >= dim)
                return{};//error, returns null
            for (i = 0; i<dim; i++)//sumo la linea que no es 0 el de la misma fila de ord
            {
                x[ord][i] += x[i2][i];
                y[ord][i] += y[i2][i];
            }
        }
        diagon = 1.0/x[ord][ord];
        ptry = &y[ord][0];
        ptrx = &x[ord][0];
        for (i = 0; i < dim; i++)
        {
            *ptry++ *= diagon;
            *ptrx++ *= diagon;
        }

        //Hacemos '0' la columna ord salvo elemento diagonal:
        for (i = 0; i<dim; i++)//Empezamos por primera fila
        {
            if (i == ord) continue;
            coef = x[i][ord];//elemento ha hacer 0 
            if (fabs(coef)<1e-15) continue; //si es cero se evita
            ptry = &y[i][0];
            ptry2 = &y[ord][0];
            ptrx = &x[i][0];
            ptrx2 = &x[ord][0];
            for (j = 0; j < dim; j++)
            {
                *ptry++ = *ptry - coef * (*ptry2++);//1ª matriz
                *ptrx++ = *ptrx - coef * (*ptrx2++);//2ª matriz
            }
        }
    }//end ord
    return y;
}


void test_5_inverse()
{
    vector< vector<double> > vec1 = {
        {0,-5,0,7,33,11,-1},
        {72,0,-11,7,9,33,5 },
        {-13,31,-5,15,29,30,24 },
        {-24,9,8,-23,31,-12,4 },
        {-3,-22,4,-24,-5,27,-10 },
        {-10,-21,-16,-32,-11,20,14 },
        {5,30,13,-32,29,-13,-13 }
    };
    vector< vector<double> > vec2;
    vec2 = inverse(vec1);
    vector< vector<double> > vec3;
    vec3 = multiply_2Dvector(vec1, vec2);

    cout << "initial array (must be unmodified):" << endl;
    print_2Dvector(vec1);

    cout << "Must be diagon array:" << endl;
    print_2Dvector(vec3," %8.3lf");
    cout << endl;
}


void test_6_inverse(int dim)
{
    vector< vector<double> > vec1(dim, vector<double>(dim));
    for (int i=0;i<dim;i++)
        for (int j = 0; j < dim; j++)
        {
            vec1[i][j] = (-1.0 + 2.0*rand() / RAND_MAX) * 10000;
        }

    vector< vector<double> > vec2;
    double ini, end;
    ini = (double)clock();
    vec2 = inverse(vec1);
    end = (double)clock();
    cout << "Time inverse =" << (end - ini) / CLOCKS_PER_SEC << endl;
    vector< vector<double> > vec3;
    vec3 = multiply_2Dvector(vec1, vec2);

    cout << "initial array (must be unmodified):" << endl;
    //print_2Dvector(vec1);

    cout << "Must be diagon array:" << endl;
    //print_2Dvector(vec3, " %8.3lf");
    cout << endl;
}

int main()
{
    vector< vector<double> > vec1;
    init_2Dvector(vec1, 10, 5);    //size_t ymax = vec1.size(),xmax = vec1[0].size();
    //test_2_dimension_vector();
    //test_one_dimension_vector();
    test_5_inverse();
    test_6_inverse(1300);
    cout << endl << "=== END ===" << endl; getchar(); 
    return 1;
}

Answer 4

Я знаю, что это немного поздно, но если вы пытаетесь выполнить какие-либо матричные вычисления для такой маленькой матрицы (60-100 строк), тогда вычисления будут выполняться намного быстрее на процессоре, чем на GPU из-за времени, которое требуется для копирования информации из основной памяти в память GPU. Если вы хотите сделать это, я бы посоветовал использовать параллельный язык, такой как OpenMP или MPI, так как это позволит вам распараллелить ваш код для ускорения вычислений на ЦП.

Answer 5

Проверьте CULA

http://www.culatools.com/ http://www.culatools.com/versions/basic

Answer 6

Первоначальный вопрос (которому сейчас 7 лет) фактически был решен спустя 4 года в статье , описывающей инверсию матрицы в CUDA на основе Гаусса-Джордана .Он пытается распределить вычисления по разным потокам и дает подробные показатели производительности для матриц размером до 2048.

Хотя это и не OpenCL, общие идеи будут легко переводиться с CUDA.