CHOLMOD редкая проблема плотного умножения - PullRequest
Новая
       54

CHOLMOD редкая проблема плотного умножения

0 голосов
/ 19 июня 2019

Я тестирую cholmod для решения разреженной матричной системы.В качестве b уравнения Ax = b я хочу использовать что-то, что приведет к известному результату, поэтому я решил объявить xe вектор единиц и b = A * xe.

Я пробовал несколько разных подходовсделать разреженное матричное умножение плотных векторов без какого-либо успеха.В качестве первого теста я использую матрицу ex15 из коллекции матриц SuitSparse 

#include "cholmod.h"
#include <iostream>
#pragma warning(disable:4996)
int main()
{
    const char* matFile = "C:/www/CholeskyMatGPU/ex15.mtx";
    const char* bfile = "C:/www/CholeskyMatGPU/b.mtx";
    FILE* bp = fopen(bfile, "w");
    FILE* fp = fopen(matFile, "r");
    cholmod_sparse *A;
    cholmod_dense *x, *xe, *b;
    cholmod_factor *L;

    cholmod_common* c = (cholmod_common*)malloc(sizeof(cholmod_common));
    cholmod_l_start(c); /* start CHOLMOD */
    c->useGPU = 0;
    c->supernodal = CHOLMOD_SUPERNODAL;

    A = cholmod_l_read_sparse(fp, c); /* read in a matrix */
    cholmod_l_print_sparse(A, "A", c); /* print the matrix */
    fclose(fp);

    double alpha[2] = { 1., 0. };
    double beta[2] = { 0., 0. };
    xe = cholmod_l_ones(A->nrow, 1, A->xtype, c); 
    b = cholmod_l_allocate_dense(A->nrow, 1, A->nrow, CHOLMOD_LONG, c); //attempt1
    b = cholmod_copy_dense(xe, c); //attempt2
    b = cholmod_l_zeros(A->nrow, 1, A->xtype, c); //attempt3
    cholmod_sdmult(A, 0, alpha, beta, xe, b, c); //b = A*xe
    cholmod_l_write_dense(bp, b, "", c);
    L = cholmod_l_analyze(A, c); /* analyze */
    cholmod_l_factorize(A, L, c); /* factorize */
    x = cholmod_l_solve(CHOLMOD_A, L, b, c); /* solve Ax=b */
}

Я ожидаю, что b будет результатом операции A * xe, поэтому что-то вроде 

-4487487.55116591
228083519.121975
0.0561846979753665
18213301.5299736
...
30383411.1062435
8073.91311889887
0.00602254606925790
-7101038.89577146
378419792.843878
0.0931721845637848

(проверено в matlab), но показанный код будет печатать вектор нулей, что-то не так с тем, что я делаю с sdmult. 

A = Problem.A;
n = length(A);
xe(1:n,1) = 1;
b=A*xe;
R=chol(A);
tmp=(R'\b);
x=R\tmp;

1 Ответ

0 голосов
/ 21 июня 2019

Ошибка была невероятно глупой, в имени метода нет длинного префикса!Это не выдаст никакой ошибки, но не выполнит умножение.С cholmod_l_sdmult(A, 0, alpha, beta, xe, b, c); вместо cholmod_sdmult(A, 0, alpha, beta, xe, b, c); все работает нормально ..

...