Intel MKL Sparse QR Solve в C ++ возвращает неинициализированную ошибку - PullRequest
3 голосов
/ 14 июня 2019

При попытке использовать mkl_sparse_s_qr_solve я получаю результат со всеми 0 и статус ошибки SPARSE_STATUS_NOT_INITIALIZED, что означает, что ручка / матрица пуста.

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

int main() {
    std::vector<int> i_b(22);
    std::vector<int> i_e(22);
    // j, v, and b vectors are just examples.
    // Regardless of the resulting numbers of the solve, I just
    // cannot get the SPARSE_STATUS_NOT_INITIALIZED error to stop
    // for the qr solver.
    std::vector<int> j(44, 0);
    std::vector<float> v(44, 1.0);
    std::vector<float> b(22, 1.0);

    {
        struct IncGenerator {
            int current_;
            IncGenerator(int start) : current_(start) {}
            int operator() () {
                current_ += 2;
                return current_;
            }
        };
        // Fill i_b with {0, 2, 4, ..., 42}
        IncGenerator g(-2);
        std::generate(i_b.begin(), i_b.end(), g);

        // Fill i_e with {2, 4, 6, ..., 44}
        IncGenerator f(0);
        std::generate(i_e.begin(), i_e.end(), f);
    }

    // ...

    // j, v, and b arrays are all the correct values
    // confirmed. The sparse A matrix should have 2 values
    // per row, with 22 rows, and 15 columns.

    int out;
    sparse_matrix_t A;
    out = mkl_sparse_s_create_csr(&A, SPARSE_INDEX_BASE_ZERO, 22, 15, &i_b[0], &i_e[0], &j[0], &v[0]);

    switch (out) {
    case SPARSE_STATUS_SUCCESS:
        std::cout << "Successfully created matrix!" << std::endl;
        break;
    case SPARSE_STATUS_NOT_INITIALIZED:
        std::cout << "Not initialized." << std::endl;
        break;
    case SPARSE_STATUS_ALLOC_FAILED:
        std::cout << "Internal memory allocation failed." << std::endl;
        break;
    default:
        std::cout << "Unknown." << std::endl;
        break;
    }

    std::vector<float> X(22 * 15);
    out = mkl_sparse_s_qr_solve(SPARSE_OPERATION_NON_TRANSPOSE, A, NULL, SPARSE_LAYOUT_COLUMN_MAJOR, 1, &X[0], 15, &asv[0], 22);
    switch (out) {
    case SPARSE_STATUS_SUCCESS:
        std::cout << "Successfully solved!" << std::endl;
        break;
    case SPARSE_STATUS_NOT_INITIALIZED:
        std::cout << "Not initialized." << std::endl;
        break;
    case SPARSE_STATUS_ALLOC_FAILED:
        std::cout << "Internal memory allocation failed." << std::endl;
        break;
    default:
        std::cout << "Unknown." << std::endl;
        break;
    }

    return 0;
}

Поэтому по какой-то причине я не могу решить с помощью A, потому что он либо думает, что A пуст, либо что-то еще неинициализировано.Я не думаю, что A пусто (я хотел проверить, но не было удобного способа печати A), поскольку инициализация матрицы A возвращается как успешная операция (единственное, в чем я немного сомневаюсь, этоначало строки i_b и конец строки i_e индексы).

Может кто-нибудь предложить какие-нибудь рекомендации?

1 Ответ

2 голосов
/ 14 июня 2019

Это не то, как вы должны использовать mkl_sparse_?_qr_solve.Разреженные системы решаются в 3 этапа (фазы):

  1. Переупорядочить.
  2. Разложить.
  3. Решить.

Сначала вынужно позвонить mkl_sparse_qr_reorder, затем mkl_sparse_?_qr_factorize, и только потом mkl_sparse_?_qr_solve:

Попробуйте вставить следующий код перед mkl_sparse_?_qr_solve:

struct matrix_descr descr;
descr.type = SPARSE_MATRIX_TYPE_GENERAL;
out = mkl_sparse_qr_reorder(A, descr);
switch (out) { ... }

out = mkl_sparse_?_qr_factorize(A, NULL);
switch (out) { ... }

Или просто использовать mkl_sparse_?_qr это сделает все 3 шага для вас за один звонок.Разделение процесса на три этапа дает вам больше свободы.Например, если вы хотите решить несколько систем с одним и тем же A, вы можете сэкономить время, вызывая mkl_sparse_qr_reorder и mkl_sparse_?_qr_factorize только один раз.

Не имеет прямого отношения, но не используйте int вместо MKL_INT.Когда определено MKL_ILP64, MKL_INT - это не int, а long long int.

...