main.cpp

#include <mkl.h>

#include <array>
#include <cmath>
#include <iomanip>
#include <iostream>

int main()
{
    constexpr MKL_INT n = 5;
    constexpr MKL_INT nnz = 9;

    // Zero-based CSR for a real, structurally nonsymmetric matrix.
    std::array<MKL_INT, n + 1> row_ptr{0, 3, 5, 7, 8, nnz};
    std::array<MKL_INT, nnz> column_indices{0, 2, 3, 1, 4, 2, 3, 3, 4};
    std::array<double, nnz> values{1.0, 2.0, 2.0, 1.0, 3.0,
                                    1.0, 2.0, 1.0, 2.0};

    const std::array<double, n> exact{1.0, 1.0, 1.0, 1.0, 1.0};
    std::array<double, n> rhs{};
    std::array<double, n> solution{};

    for (MKL_INT row = 0; row < n; ++row) {
        for (MKL_INT position = row_ptr[row];
             position < row_ptr[row + 1]; ++position) {
            rhs[row] += values[position] * exact[column_indices[position]];
        }
    }

    void* internal_data[64]{};
    MKL_INT parameters[64]{};
    MKL_INT matrix_type = 11;  // real and nonsymmetric
    pardisoinit(internal_data, &matrix_type, parameters);
    parameters[34] = 1;  // use zero-based CSR indices

    MKL_INT max_factorizations = 1;
    MKL_INT matrix_number = 1;
    MKL_INT phase = 13;  // analyze, factorize, and solve
    MKL_INT right_hand_sides = 1;
    MKL_INT message_level = 0;
    MKL_INT error = 0;

    pardiso(internal_data,
            &max_factorizations,
            &matrix_number,
            &matrix_type,
            &phase,
            &n,
            values.data(),
            row_ptr.data(),
            column_indices.data(),
            nullptr,
            &right_hand_sides,
            parameters,
            &message_level,
            rhs.data(),
            solution.data(),
            &error);

    if (error != 0) {
        std::cerr << "PARDISO solve failed with error " << error << '\n';
    }

    phase = -1;
    MKL_INT release_error = 0;
    pardiso(internal_data,
            &max_factorizations,
            &matrix_number,
            &matrix_type,
            &phase,
            &n,
            nullptr,
            row_ptr.data(),
            column_indices.data(),
            nullptr,
            &right_hand_sides,
            parameters,
            &message_level,
            nullptr,
            nullptr,
            &release_error);

    if (error != 0) return 1;
    if (release_error != 0) {
        std::cerr << "PARDISO cleanup failed with error " << release_error << '\n';
        return 1;
    }

    double error_squared = 0.0;
    double exact_squared = 0.0;
    for (MKL_INT index = 0; index < n; ++index) {
        const double difference = solution[index] - exact[index];
        error_squared += difference * difference;
        exact_squared += exact[index] * exact[index];
    }
    const double relative_error = std::sqrt(error_squared / exact_squared);

    std::cout << "Solution:";
    for (const double value : solution)
        std::cout << ' ' << std::setprecision(12) << value;
    std::cout << "\nRelative solution error: " << relative_error << '\n';

    return relative_error < 1.0e-12 ? 0 : 1;
}