python — Сингулярное разложение по значению с LAPACK: проблемы с большими матрицами

Я использую интерфейс C LAPACK для вычисления разложения по сингулярным числам (SVD) матрицы. Для этого я использую рутину dgesvd_,

Я создал простой скрипт C ++, который создает случайную матрицу (с M строки и N столбцы), и это вычисляет его SVD. Код этого скрипта следующий:

#include <iostream>
#include <random>
#include "lapacke.h"
#define M 150
#define N 3
#define MEAN 0
#define STD 0.25int main(int argc, char *argv[])
{

std::default_random_engine generator;
std::normal_distribution<double> distribution(MEAN, STD);

int m = M, n = N, lda = M, ldu = M, ldvt = N, info, lwork;
double wkopt;
double* work;
double s[n], u[ldu * m], vt[ldvt * n], a[lda * n];

for(unsigned int k = 0; k < (lda * n); k++){
a[k] = distribution(generator);
}

lwork = -1;
dgesvd_("All", "All", &m, &n, a, &lda, s, u, &ldu, vt, &ldvt, &wkopt, &lwork, &info);
lwork = (int)wkopt;
work = (double*) malloc(lwork * sizeof(double));
dgesvd_("All", "All", &m, &n, a, &lda, s, u, &ldu, vt, &ldvt, work, &lwork, &info);

if(info > 0){
std::cerr<< "The algorithm computing SVD failed to converge\n"     ;
exit(1);
}

free(work);

std::cout<<"Singular values:\n";
for(unsigned int k = 0; k < n; k++){
std::cout<<s[k]<<" ";
}
std::cout<<"\n";
return(0);
}

Матрица из 150 строк и 3 столбцов, кажется, правильно вычисляет SVD. Однако, когда матрица имеет большее количество строк (например, матрица с 1500 строками и 3 столбцами), выполнение скомпилированного скрипта вызывает такую ​​ошибку: Ошибка сегментации (ядро сброшено).

Я пытался сделать что-то подобное в скрипте Python:

import numpy as np
M = 1500
N = 3
MEAN = np.zeros(3)
STD = 0.25
result = np.linalg.svd(points)

Когда я использую Python, значения единственного числа, кажется, вычисляются правильно без какой-либо ошибки, хотя метод NumPy, который я использовал, выполняется с помощью процедуры LAPACK dgesdd_ (Я пробовал это и та же ошибка возникает).

Кто-нибудь знает, почему происходит эта ошибка? Любая помощь для решения проблемы будет принята с благодарностью. Благодарю.

Pd: я использую версию 3.6.0 lapacke и Ubuntu 16.04 LTS.

0

Решение

Распределение матриц по куче таким образом решает проблему:

double* s = new double[n];
double* u = new double[ldu * m];
double* vt = new double[ldvt * n];
double* a = new double[lda * n];
0

Другие решения

Других решений пока нет …

По вопросам рекламы [email protected]