Эффективное решение для использования памяти C ++ для системы линейной алгебры Ax = b

Question

Я использую привязки числовой библиотеки для Boost UBLAS для решения простой линейной системы.Следующее работает нормально, за исключением того, что оно ограничено обработкой матриц A (m x m) для относительно небольших 'm'.

На практике у меня есть гораздо большая матрица с размерностью m = 10 ^ 6 (до 10 ^ 7).
Существует ли существующий подход C ++ для решения Ax = b, который эффективно использует память?

#include<boost/numeric/ublas/matrix.hpp>
#include<boost/numeric/ublas/io.hpp>
#include<boost/numeric/bindings/traits/ublas_matrix.hpp>
#include<boost/numeric/bindings/lapack/gesv.hpp>
#include <boost/numeric/bindings/traits/ublas_vector2.hpp>

// compileable with this command


//g++ -I/home/foolb/.boost/include/boost-1_38 -I/home/foolb/.boostnumbind/include/boost-numeric-bindings solve_Axb_byhand.cc -o solve_Axb_byhand -llapack


namespace ublas = boost::numeric::ublas;
namespace lapack= boost::numeric::bindings::lapack;


int main()
{
    ublas::matrix<float,ublas::column_major> A(3,3);
    ublas::vector<float> b(3);


    for(unsigned i=0;i < A.size1();i++)
        for(unsigned j =0;j < A.size2();j++)
        {
            std::cout << "enter element "<<i << j << std::endl;
            std::cin >> A(i,j);
        }

    std::cout << A << std::endl;

    b(0) = 21; b(1) = 1; b(2) = 17;

    lapack::gesv(A,b);

    std::cout << b << std::endl;


    return 0;
}

Answer 1

Краткий ответ:Не используйте Boost'ы LAPACK привязки, они были разработаны для плотных матриц, не для разреженных матриц, используйте UMFPACK вместо этого.

Длинный ответ: UMFPACK это одна из лучших библиотек для решения Ax = b, когда A большое и разреженное.

• http://www.cise.ufl.edu/research/sparse/umfpack/
• http://www.cise.ufl.edu/research/sparse/umfpack/UMFPACK/Doc/QuickStart.pdf

Ниже приведен пример кода (основанный на umfpack_simple.c), который генерирует простой A и b и решает Ax = b.

#include <stdlib.h>
#include <stdio.h>
#include "umfpack.h"

int    *Ap; 
int    *Ai;
double *Ax; 
double *b; 
double *x; 

/* Generates a sparse matrix problem: 
   A is n x n tridiagonal matrix
   A(i,i-1) = -1;
   A(i,i) = 3; 
   A(i,i+1) = -1; 
*/
void generate_sparse_matrix_problem(int n){
  int i;  /* row index */ 
  int nz; /* nonzero index */
  int nnz = 2 + 3*(n-2) + 2; /* number of nonzeros*/
  int *Ti; /* row indices */ 
  int *Tj; /* col indices */ 
  double *Tx; /* values */ 

  /* Allocate memory for triplet form */
  Ti = malloc(sizeof(int)*nnz);
  Tj = malloc(sizeof(int)*nnz);
  Tx = malloc(sizeof(double)*nnz);

  /* Allocate memory for compressed sparse column form */
  Ap = malloc(sizeof(int)*(n+1));
  Ai = malloc(sizeof(int)*nnz);
  Ax = malloc(sizeof(double)*nnz);

  /* Allocate memory for rhs and solution vector */
  x = malloc(sizeof(double)*n);
  b = malloc(sizeof(double)*n);

  /* Construct the matrix A*/
  nz = 0;
  for (i = 0; i < n; i++){
    if (i > 0){
      Ti[nz] = i;
      Tj[nz] = i-1;
      Tx[nz] = -1;
      nz++;
    }

    Ti[nz] = i;
    Tj[nz] = i;
    Tx[nz] = 3;
    nz++;

    if (i < n-1){
      Ti[nz] = i;
      Tj[nz] = i+1;
      Tx[nz] = -1;
      nz++;
    }
    b[i] = 0;
  }
  b[0] = 21; b[1] = 1; b[2] = 17;
  /* Convert Triplet to Compressed Sparse Column format */
  (void) umfpack_di_triplet_to_col(n,n,nnz,Ti,Tj,Tx,Ap,Ai,Ax,NULL);

  /* free triplet format */ 
  free(Ti); free(Tj); free(Tx);
}


int main (void)
{
    double *null = (double *) NULL ;
    int i, n;
    void *Symbolic, *Numeric ;
    n = 500000;
    generate_sparse_matrix_problem(n);
    (void) umfpack_di_symbolic (n, n, Ap, Ai, Ax, &Symbolic, null, null);
    (void) umfpack_di_numeric (Ap, Ai, Ax, Symbolic, &Numeric, null, null);
    umfpack_di_free_symbolic (&Symbolic);
    (void) umfpack_di_solve (UMFPACK_A, Ap, Ai, Ax, x, b, Numeric, null, null);
    umfpack_di_free_numeric (&Numeric);
    for (i = 0 ; i < 10 ; i++) printf ("x [%d] = %g\n", i, x [i]);
    free(b); free(x); free(Ax); free(Ai); free(Ap);
    return (0);
}

Функция generate_sparse_matrix_problem создает матрицу A и правая сторона b.Матрица сначала строится в триплетной форме. Векторы Ti, Tj и Tx полностью описывают A.Триплетную форму легко создать, но эффективные методы разреженных матриц требуют сжатого формата разреженных столбцов.Преобразование выполняется с помощью umfpack_di_triplet_to_col.

Символическая факторизация выполняется с помощью umfpack_di_symbolic.Разреженное LU-разложение A выполняется с umfpack_di_numeric.Нижнее и верхнее треугольные решения выполняются с помощью umfpack_di_solve.

С n как и в случае с 500 000, на моем компьютере выполнение всей программы занимает около секунды.Valgrind сообщает, что было выделено 369 239 649 байт (чуть более 352 МБ).

Обратите внимание на это Страница обсуждается поддержка Boost для разреженных матриц в триплетном (координатном) и сжатом формате.Если хотите, вы можете написать процедуры для преобразования этих объектов boost в простые массивы UMFPACK требуется в качестве входных данных.

Answer 2

Предполагая, что ваши огромные матрицы разрежены, а я надеюсь, что они такого размера, взгляните на ПАРДИСО проект, представляющий собой разреженный линейный решатель, - это то, что вам понадобится, если вы хотите обрабатывать матрицы такого размера, как вы сказали.Позволяет эффективно хранить только ненулевые значения и намного быстрее, чем решение той же системы плотных матриц.

Answer 3

Я предполагаю, что ваша матрица плотная.Если он разрежен, вы можете найти множество специализированных алгоритмов, как уже упоминалось выше. ДеусАдуро и даффимо.

Если в вашем распоряжении нет (достаточно большого) кластера, вы хотите взглянуть на неосновные алгоритмы. Скалапак имеет несколько неосновных решателей в составе своей пакет прототипов, смотрите документацию здесь и Google для получения более подробной информации.Поиск в Интернете "неосновных LU / (матричных) решателей / пакетов" даст вам ссылки на множество дополнительных алгоритмов и инструментов.Я не специалист в этих вопросах.

Однако для решения этой проблемы большинство людей использовали бы кластер.Пакет, который вы найдете практически в любом кластере, опять же ScaLAPACK.Кроме того, в типичном кластере обычно имеется множество других пакетов, так что вы можете выбрать то, что подходит для вашей задачи (примеры здесь и здесь).

Прежде чем приступить к программированию, вы, вероятно, захотите быстро проверить, сколько времени потребуется для решения вашей проблемы.Типичный решатель занимает около O (3 * N ^ 3) флопов (N - размерность матрицы).Если N = 100000, то, следовательно, вы смотрите на 3000000 Гфлопс.Предполагая, что ваш решатель в памяти выполняет 10 Гфлопс / с на ядро, вы рассчитываете на 3 1/2 дня работы на одном ядре.Поскольку алгоритмы хорошо масштабируются, увеличение количества ядер должно сокращать время почти линейно.Вдобавок ко всему идет ввод-вывод.

Answer 4

Не уверен насчет реализаций C ++, но есть несколько вещей, которые вы можете сделать, если возникает проблема с памятью, в зависимости от типа матрицы, с которой вы имеете дело:

1. Если ваша матрица разреженная или полосатая, вы можете использовать решатель разреженной или пропускной способности.Они не хранят нулевые элементы за пределами диапазона.
2. Вы можете использовать решатель волнового фронта, который сохраняет матрицу на диске и вводит волновой фронт матрицы только для декомпозиции.
3. Вы можете вообще избежать решения матрицы и использовать итерационные методы.
4. Вы можете попробовать решить проблему методами Монте-Карло.

Answer 5

Взгляните на список свободно доступного программного обеспечения для решения задач линейной алгебры, составленный Джеком Донгаррой и Хатемом Лтаефом.

Я думаю, что для того размера задачи, на который вы смотрите, вам, вероятно, нужен итеративный алгоритм.Если вы не хотите хранить матрицу A в разреженном формате, вы можете использовать реализацию без матрицы.Итеративным алгоритмам обычно не требуется обращаться к отдельным элементам матрицы A, им нужно только вычислить произведения матрицы на вектор Av (а иногда и A ^ T v, произведение транспонированной матрицы на вектор).Таким образом, если библиотека хорошо спроектирована, этого должно быть достаточно, если вы передадите ей класс, который знает, как выполнять матрично-векторные произведения.

Answer 6

Как следует из принятого ответа, существует UMFPACK.Но если вы используете BOOST, вы все равно можете использовать компактные матрицы в BOOST и использовать UMFPACK для решения системы.Существует привязка, которая упрощает это:

http://mathema.tician.de/software/boost-numeric-bindings

Оно датировано примерно двумя годами, но это всего лишь переплет (наряду с несколькими другими).

смотрите связанный вопрос:UMFPACK и разреженная матрица uBLAS от BOOST