Intel MKL Sparse QR Solve in C++ returns not initialized error

user_11595345 0
Viewed 377

When attempting to use mkl_sparse_s_qr_solve, I receive a result of all 0's and the error status of SPARSE_STATUS_NOT_INITIALIZED which means that the handle/matrix is empty.

I have tried reading through the documentation thoroughly and printing all arrays that are used to instantiate the CSR sparse matrix required for the solve, and all the arrays contain the correct values.

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;
}

Therefore, for some reason I cannot solve with A because it either thinks A is empty or something else is uninitialized. I do not think A is empty (I wanted to check but there was no convenient way to print A) as the initialization of the matrix A returns as a successful operation (the only thing I am slightly doubtful of is the row beginning i_b and row end i_e indices).

Can anyone please offer some guidance?

intel-mkl sparse-matrix c++
1 Answers

This is not how you are supposed to use mkl_sparse_?_qr_solve. Sparse systems are solved in 3 steps (phases):

  1. Reorder.
  2. Factorize.
  3. Solve.

First, you have to call mkl_sparse_qr_reorder, then mkl_sparse_?_qr_factorize, and only then mkl_sparse_?_qr_solve:

Try to insert the following code before 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) { ... }

Or just use mkl_sparse_?_qr that will do all 3 steps for you in a single call. Separation of the process into three steps gives you more freedom. For example, if you want to solve several systems with the same A, you can save time by calling mkl_sparse_qr_reorder and mkl_sparse_?_qr_factorize only once.

Not directly related, but don't use int instead of MKL_INT. When MKL_ILP64 is defined, MKL_INT is not int, but long long int.

user_16251870
Related