Как мне умножить две разреженные матрицы в C?

Question

У меня есть разреженная матрица D, и я хочу умножить D_transpose и D, чтобы получить L следующим образом:

L = D '* D;

Я использую sparseBLAS для работы с разреженными матрицами, но в документации сказано, что нечего умножать две разреженные матрицы.

Я полностью застрял и не знаю, как поступить.

Размеры D обычно составляют около 500 000 на 250 000. Я вообще не могу выделить столько памяти, так что это нужно сделать с помощью разреженных матриц.

Я сделал это с помощью MATLAB, и я не понимаю, как MATLAB делает это, если он также использует sparseBLAS под интерфейсом - или нет? Если нет, что он использует? Я тоже мог бы это использовать.

Спасибо за чтение!

РЕДАКТИРОВАТЬ: Решено. Мне нужно, чтобы матрица L умножалась на вектор. Таким образом, вместо того, чтобы сначала вычислить L, я просто сделал D'*(D*x), таким образом избегая необходимости умножения двух разреженных матриц. Я сейчас делаю только разреженные матрицы и плотные векторные умножения, которые поддерживаются sparseBLAS.

Answer 1

Библиотека Intel MKL поддерживает произведение двух разреженных матриц. Вы можете проверить режим mkl_sparse_spmmиз справки МКЛ. (Старая версия MKL использовала процедуру mkl_?csrmultcsrкоторый сейчас устарел.)

Кроме того, существует старая, но широко используемая поточно-ориентированная библиотека разреженных матриц с открытым исходным кодом. SPARSKIT(и его обновленная версия SPARSKIT2). Оба они написаны F77. есть рутина amubкоторый можно использовать для получения произведения двух разреженных матриц. Вы можете проверить, не переписал ли кто-нибудь эту библиотеку на C.

Answer 2

Это на самом деле указано в размещенной документации.

Страница 11

5.2 Использование разреженных матриц BLAS

После того, как дескриптор матрицы Sparse BLAS полностью построен (что можно проверить, проверив свойство blas_valid_handle), можно использовать дескриптор матрицы для выполнения операций. В настоящее время поддерживаются четыре операции, показанные в таблицах 3.2 и 3.3. В дополнение к выполнению операций с матрицей Sparse BLAS, можно запрашивать ее свойства через его дескриптор. В таблице 5.5 перечислены свойства, которые можно получить, вызвав процедуру get properties.

Таблица 3.3 Страница 4

USMM редкое матрично-матричное умножение

Так что поддержка вроде бы есть. Я просто не могу найти подпись для BLAS_usmm функция. Может быть, вы можете проверить в заголовках.

Изменить: Если вы получили свой sparseBLas от NIST, вы можете проверить blas_sparse_proto.h файл для BLAS_*usmm функции для подписей и параметров.

 /* Level 3 Computational Routines */

int BLAS_susmm( enum blas_order_type order, enum blas_trans_type transa,
    int nrhs, float alpha, blas_sparse_matrix A, const float *b, int ldb,
        float *c, int ldc );
int BLAS_dusmm( enum blas_order_type order, enum blas_trans_type transa,
        int nrhs, double alpha, blas_sparse_matrix A, const double *b,
        int ldb, double *c, int ldc );
int BLAS_cusmm( enum blas_order_type order, enum blas_trans_type transa,
         int nrhs, const void *alpha, blas_sparse_matrix A, const void *b, 
     int ldb, void *c, int ldc );
int BLAS_zusmm( enum blas_order_type order, enum blas_trans_type transa,
         int nrhs, const void *alpha, blas_sparse_matrix A, const void *b, 
     int ldb, void *c, int ldc );

Answer 3

Насколько я понимаю, ваша проблема в основном заключается в хранении огромной матрицы в памяти. Вы можете хранить значения в паре (строка, столбец). Например,

1 0 0
0 0 2
0 4 0

Эта матрица может быть сохранена в std::map<pair<int, int>, int> как:

map[make_pair(1, 1)] = 1
map[make_pair(2, 3)] = 2
map[make_pair(3, 2)] = 4

Теперь часть вычислений. Предполагая, что первая матрица хранится в map1вторая матрица хранится в map2 и ответ хранится в mapAns,

for each element x in map1:
    for each element y in map2:
        if x.column == y.row:
            mapAns[x.row, y.column] += x.value * y.value

Вам нужно использовать карту как пользовательскую структуру данных, если вы хотите сделать то же самое в C.