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

Question

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

Возможно, Knuth относится к ОСНОВНОМУ.

В ОСНОВНОМ (более старый, я не знаю о VB), ДЛЯ циклов были реализованы по-другому.

, Например, циклично выполняться десять раз

ДЛЯ N=1 К 10

...

NEXT N

----ИЛИ----для нахождения суммы четных чисел ниже 100

ДЛЯ N=0 К 100 шагам 2

...

NEXT N

В C у Вас может быть любое условие в "для" проверить на выход из цикла. Более гибкий, я думаю.

9

c++ boost linear-algebra lapack umfpack

задан denfromufa 12 October 2013 в 10:28

5 ответов

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

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

3

ответ дан 4 December 2019 в 06:57

Предполагая, что ваши огромные матрицы являются разреженными, а я надеюсь, что они имеют такой размер, взгляните на проект PARDISO , который является разреженным линейным решателем. Это то, что вам понадобится, если вы хотите обрабатывать матрицы такого размера, как вы сказали. Позволяет эффективно хранить только ненулевые значения и намного быстрее, чем решение той же системы плотных матриц.

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

6

ответ дан 4 December 2019 в 06:57

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

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

3

ответ дан 4 December 2019 в 06:57

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

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

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

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

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

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

увеличение количества ядер должно уменьшить время, близкое к линейному. Вдобавок ко всему идет ввод-вывод.

6

ответ дан 4 December 2019 в 06:57

Другие вопросы по тегам:

c++ boost linear-algebra lapack umfpack

Похожие вопросы:

score 14 · Accepted Answer

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

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

Ниже приведен образец код (на основе 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 как 500000 на моей машине выполнение всей программы занимает около секунды. Valgrind сообщает, что было выделено 369 239 649 байтов (немногим более 352 МБ).

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