amgx_ilu_solver.cpp

#include <iostream>
#include <fstream>
#include <vector>
#include <amgx_c.h>
#include "mmio.h"

void read_matrix_market_A(const char* filename, std::vector<int>& row, std::vector<int>& col, std::vector<double>& val,
                          int &n, int &nnz) {
    FILE* f = fopen(filename, "r");
    MM_typecode matcode;
    mm_read_banner(f, &matcode);
    int m, k;
    mm_read_mtx_crd_size(f, &m, &k, &nnz);
    n = m;
    row.resize(nnz);
    col.resize(nnz);
    val.resize(nnz);
    for (int i = 0; i < nnz; ++i) {
        int r, c;
        double v;
        fscanf(f, "%d %d %lf", &r, &c, &v);
        row[i] = r - 1;
        col[i] = c - 1;
        val[i] = v;
    }
    fclose(f);
}

void read_vector_b(const char* filename, std::vector<double>& b, int &n) {
    FILE* f = fopen(filename, "r");
    MM_typecode matcode;
    mm_read_banner(f, &matcode);
    int m, cols;
    fscanf(f, "%d %d", &m, &cols);
    n = m;
    b.resize(n);
    for (int i = 0; i < n; ++i) {
        fscanf(f, "%lf", &b[i]);
    }
    fclose(f);
}

int main() {
    AMGX_config_handle cfg;
    AMGX_resources_handle rsrc;
    AMGX_matrix_handle A;
    AMGX_vector_handle x, b;

    // Init
    AMGX_initialize();
    AMGX_config_create_from_file(&cfg, "solver_config.json");
    AMGX_resources_create_simple(&rsrc, cfg);

    // Read matrix
    std::vector<int> row, col;
    std::vector<double> val;
    int n, nnz;
    read_matrix_market_A("A.mtx", row, col, val, n, nnz);

    // Read RHS
    std::vector<double> b_vec;
    int rhs_len;
    read_vector_b("b.mtx", b_vec, rhs_len);

    // Create matrix and vectors
    AMGX_vector_create(&x, rsrc, AMGX_mode_dDDI);
    AMGX_vector_create(&b, rsrc, AMGX_mode_dDDI);
    AMGX_matrix_create(&A, rsrc, AMGX_mode_dDDI);

    // Upload system
    AMGX_matrix_upload_all(A, n, nnz, 1, 1, row.data(), col.data(), val.data(), NULL);

    AMGX_vector_upload(b, n, 1, b_vec.data());
    std::vector<double> x_vec(n, 0.0);
    AMGX_vector_upload(x, n, 1, x_vec.data());

    // Create solver and solve
    AMGX_solver_handle solver;
    AMGX_solver_create(&solver, rsrc, AMGX_mode_dDDI, cfg);
    AMGX_solver_setup(solver, A);
    AMGX_solver_solve(solver, b, x);
    AMGX_SOLVE_STATUS status;
    AMGX_solver_get_status(solver, &status);

    if (status != AMGX_SOLVE_SUCCESS)
    {
        std::cerr << "❌ AMGX solve failed. Status code: " << status << std::endl;
    }
    else
    {
        std::cout << "✅ AMGX solve succeeded." << std::endl;
    }

    // Download solution
    AMGX_vector_download(x, x_vec.data());

    std::ofstream out("solution.mtx");
    out << "%%MatrixMarket matrix array real general\n";
    out << n << " 1\n";
    for (int i = 0; i < n; ++i)
        out << x_vec[i] << "\n";
    out.close();

    std::cout << "Solved Ax = b and the solution written to solution.mtx\n";

    // Cleanup
    AMGX_solver_destroy(solver);
    AMGX_vector_destroy(x);
    AMGX_vector_destroy(b);
    AMGX_matrix_destroy(A);
    AMGX_resources_destroy(rsrc);
    AMGX_config_destroy(cfg);
    AMGX_finalize();
    return 0;
}