aoclsparse_itsol_?_rci_solve() - 5.2 English - 68552

AOCL API Guide (68552)
Document ID
68552
Release Date
2025-12-29
Version
5.2 English
aoclsparse_status aoclsparse_itsol_s_rci_solve(aoclsparse_itsol_handle handle, aoclsparse_itsol_rci_job *ircomm, float **u, float **v, float *x, float rinfo[100])#
aoclsparse_status aoclsparse_itsol_d_rci_solve(aoclsparse_itsol_handle handle, aoclsparse_itsol_rci_job *ircomm, double **u, double **v, double *x, double rinfo[100])#
aoclsparse_status aoclsparse_itsol_c_rci_solve(aoclsparse_itsol_handle handle, aoclsparse_itsol_rci_job *ircomm, aoclsparse_float_complex **u, aoclsparse_float_complex **v, aoclsparse_float_complex *x, float rinfo[100])#
aoclsparse_status aoclsparse_itsol_z_rci_solve(aoclsparse_itsol_handle handle, aoclsparse_itsol_rci_job *ircomm, aoclsparse_double_complex **u, aoclsparse_double_complex **v, aoclsparse_double_complex *x, double rinfo[100])#

Reverse Communication Interface (RCI) to the iterative solvers (itsol) suite.

This function solves the linear system of equations

\[ Ax=b, \]
where the matrix of coefficients \(A\) is not required to be provided explicitly. The right hand-side is the dense vector b and the vector of unknowns is x. If \(A\) is symmetric and positive definite then set the option “iterative method” to “cg” to solve the problem using the Conjugate Gradient method, alternatively set the option to “gmres” to solve using GMRes. See the Options for a list of available options to modify the behaviour of each solver.

The reverse communication interface (RCI), also know as _matrix-free_ interface does not require the user to explicitly provide the matrix \(A\). During the solve process whenever the algorithm requires a matrix operation (matrix-vector or transposed matrix-vector products), it returns control to the user with a flag ircomm indicating what operation is requested. Once the user performs the requested task it must call this function again to resume the solve.

The expected workflow is as follows:

  1. Call aoclsparse_itsol_s_init or aoclsparse_itsol_d_init or aoclsparse_itsol_c_init or aoclsparse_itsol_z_init to initialize the problem handle ( aoclsparse_itsol_handle)

  2. Choose the solver and adjust its behaviour by setting optional parameters with aoclsparse_itsol_option_set, see also Options.

  3. Define the problem size and right-hand side vector \(b\) with aoclsparse_itsol_d_rci_input.

  4. Solve the system with either aoclsparse_itsol_s_rci_solve or aoclsparse_itsol_d_rci_solve or aoclsparse_itsol_c_rci_solve or aoclsparse_itsol_z_rci_solve.

  5. If there is another linear system of equations to solve with the same matrix but a different right-hand side \(b\), then repeat from step 3.

  6. If solver terminated successfully then vector x contains the solution.

  7. Free the memory with aoclsparse_itsol_destroy.

These reverse communication interfaces complement the _forward communication_ interfaces aoclsparse_itsol_d_rci_solve, aoclsparse_itsol_s_rci_solve, aoclsparse_itsol_c_rci_solve and aoclsparse_itsol_z_rci_solve.

  1/* ************************************************************************
  2 * Copyright (c) 2022-2025 Advanced Micro Devices, Inc.
  3 *
  4 * Permission is hereby granted, free of charge, to any person obtaining a copy
  5 * of this software and associated documentation files (the "Software"), to deal
  6 * in the Software without restriction, including without limitation the rights
  7 * to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
  8 * copies of the Software, and to permit persons to whom the Software is
  9 * furnished to do so, subject to the following conditions:
 10 *
 11 * The above copyright notice and this permission notice shall be included in
 12 * all copies or substantial portions of the Software.
 13 *
 14 * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
 15 * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
 16 * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
 17 * AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
 18 * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
 19 * OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
 20 * THE SOFTWARE.
 21 *
 22 * ************************************************************************ */
 23
 24#include "aoclsparse.h"
 25
 26#include <iomanip>
 27#include <iostream>
 28#include <math.h>
 29#include <vector>
 30
 31// Define custom log printer
 32void printer(double rinfo[100], bool header)
 33{
 34    std::ios_base::fmtflags fmt = std::cout.flags();
 35    fmt |= std::cout.scientific | std::cout.right;
 36    if(header)
 37        std::cout << std::setw(5) << std::right << " iter"
 38                  << " " << std::setw(16) << std::right << "optim"
 39                  << "  " << std::setw(3) << std::endl;
 40    std::cout << std::setw(5) << std::right << (int)rinfo[30] << " " << std::setw(16) << std::right
 41              << std::scientific << rinfo[0] << "  " << std::endl;
 42    std::cout << std::resetiosflags(fmt);
 43}
 44
 45int main()
 46{
 47    // CSR symmetric matrix. Only the lower triangle is stored
 48    std::vector<aoclsparse_int> icrow, icol;
 49    std::vector<double>         aval;
 50    aoclsparse_int              n = 8, nnz = 18;
 51    icrow.assign({0, 1, 2, 5, 6, 8, 11, 15, 18});
 52    icol.assign({0, 1, 0, 1, 2, 3, 1, 4, 0, 4, 5, 0, 3, 4, 6, 2, 5, 7});
 53    aval.assign({19, 10, 1, 8, 11, 13, 2, 11, 2, 1, 9, 7, 9, 5, 12, 5, 5, 9});
 54
 55    // create matrix and its descriptor
 56    aoclsparse_matrix     A;
 57    aoclsparse_index_base base = aoclsparse_index_base_zero;
 58    aoclsparse_mat_descr  descr_a;
 59    aoclsparse_operation  trans = aoclsparse_operation_none;
 60    aoclsparse_create_dcsr(&A, base, n, n, nnz, icrow.data(), icol.data(), aval.data());
 61    aoclsparse_create_mat_descr(&descr_a);
 62    aoclsparse_set_mat_type(descr_a, aoclsparse_matrix_type_symmetric);
 63    aoclsparse_set_mat_fill_mode(descr_a, aoclsparse_fill_mode_lower);
 64    aoclsparse_set_mat_index_base(descr_a, base);
 65    aoclsparse_set_sv_hint(A, trans, descr_a, 100);
 66    aoclsparse_optimize(A);
 67
 68    // Initialize initial point x0 and right hand side b
 69    std::vector<double> x, b, expected_sol, y;
 70    x.assign({1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0});
 71    b.assign({0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0});
 72    expected_sol.assign({1.E0, 0.E0, 1.E0, 0.E0, 1.E0, 0.E0, 1.E0, 0.E0});
 73    double alpha = 1.0, beta = 0.;
 74    y.resize(n);
 75    aoclsparse_dmv(trans, &alpha, A, descr_a, expected_sol.data(), &beta, b.data());
 76
 77    // create CG handle
 78    aoclsparse_itsol_handle handle = nullptr;
 79    aoclsparse_itsol_d_init(&handle);
 80
 81    // Change options, set user defined preconditioner
 82    if(aoclsparse_itsol_option_set(handle, "CG Iteration Limit", "50") != aoclsparse_status_success
 83       || aoclsparse_itsol_option_set(handle, "CG preconditioner", "user")
 84              != aoclsparse_status_success)
 85        std::cout << "Warning an option could not be set" << std::endl;
 86
 87    // initialize size and rhs inside the handle
 88    aoclsparse_itsol_d_rci_input(handle, n, b.data());
 89
 90    // Call CG solver
 91    aoclsparse_status        status;
 92    aoclsparse_itsol_rci_job ircomm = aoclsparse_rci_start;
 93    double                  *u      = nullptr;
 94    double                  *v      = nullptr;
 95    double                   rinfo[100];
 96    double                   tol = 1.0e-5;
 97    bool                     hdr;
 98    std::cout << std::endl;
 99    while(ircomm != aoclsparse_rci_stop)
100    {
101        status = aoclsparse_itsol_d_rci_solve(handle, &ircomm, &u, &v, x.data(), rinfo);
102        if(status != aoclsparse_status_success)
103            break;
104        switch(ircomm)
105        {
106        case aoclsparse_rci_mv:
107            // Compute v = Au
108            beta   = 0.0;
109            alpha  = 1.0;
110            status = aoclsparse_dmv(trans, &alpha, A, descr_a, u, &beta, v);
111            if(status != aoclsparse_status_success)
112                ircomm = aoclsparse_rci_stop;
113            break;
114
115        case aoclsparse_rci_precond:
116            // apply Symmetric Gauss-Seidel preconditioner step
117            status = aoclsparse_dtrsv(aoclsparse_operation_none, alpha, A, descr_a, u, y.data());
118            if(status != aoclsparse_status_success)
119            {
120                ircomm = aoclsparse_rci_stop;
121                break;
122            }
123            for(aoclsparse_int i = 0; i < n; i++)
124                y[i] *= aval[icrow[i + 1] - 1];
125            status
126                = aoclsparse_dtrsv(aoclsparse_operation_transpose, alpha, A, descr_a, y.data(), v);
127            if(status != aoclsparse_status_success)
128                ircomm = aoclsparse_rci_stop;
129            break;
130
131        case aoclsparse_rci_stopping_criterion:
132            // No operations required, can be used to monitor the progress of the solve
133            // or defining a custom stopping criterion
134            // print iteration log
135            hdr = ((int)rinfo[30] % 100) == 0;
136            printer(rinfo, hdr);
137            // request solver to stop if custom criterion is met
138            if(rinfo[0] < tol)
139            {
140                std::cout << "User stop. Final residual: " << rinfo[0] << std::endl;
141                ircomm = aoclsparse_rci_interrupt;
142            }
143            break;
144
145        default:
146            break;
147        }
148    }
149    // Print the final results if the internal stopping criterion or the user defined one were met
150    switch(status)
151    {
152    case aoclsparse_status_user_stop:
153    case aoclsparse_status_success:
154        std::cout.precision(2);
155        std::cout << std::scientific;
156        aoclsparse_itsol_handle_prn_options(handle);
157        std::cout << std::endl
158                  << "Solution found: (residual = " << rinfo[0] << " in " << (int)rinfo[30]
159                  << " iterations)" << std::endl
160                  << "   Final X* = ";
161        for(int i = 0; i < n; i++)
162            std::cout << std::setw(9) << x[i] << " ";
163        std::cout << std::endl;
164        std::cout << "Expected X* = ";
165        for(int i = 0; i < n; i++)
166            std::cout << std::setw(9) << expected_sol[i] << " ";
167        std::cout << std::endl;
168        break;
169
170    case aoclsparse_status_maxit:
171        std::cout << "solve stopped after " << (int)rinfo[30] << " iterations" << std::endl
172                  << "residual = " << rinfo[0] << std::endl;
173        break;
174
175    default:
176        std::cout << "Something unexpected happened! " << status << std::endl;
177    }
178    aoclsparse_itsol_destroy(&handle);
179    aoclsparse_destroy_mat_descr(descr_a);
180    aoclsparse_destroy(&A);
181
182    return 0;
183}


  1/* ************************************************************************
  2 * Copyright (c) 2022-2023 Advanced Micro Devices, Inc.
  3 *
  4 * Permission is hereby granted, free of charge, to any person obtaining a copy
  5 * of this software and associated documentation files (the "Software"), to deal
  6 * in the Software without restriction, including without limitation the rights
  7 * to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
  8 * copies of the Software, and to permit persons to whom the Software is
  9 * furnished to do so, subject to the following conditions:
 10 *
 11 * The above copyright notice and this permission notice shall be included in
 12 * all copies or substantial portions of the Software.
 13 *
 14 * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
 15 * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
 16 * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
 17 * AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
 18 * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
 19 * OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
 20 * THE SOFTWARE.
 21 *
 22 * ************************************************************************ */
 23
 24#include "aoclsparse.h"
 25
 26#include <cfloat>
 27#include <iomanip>
 28#include <iostream>
 29#include <math.h>
 30#include <vector>
 31
 32#define PREMATURE_STOP_TOLERANCE 1e-4
 33const char gmres_sol[] = "gmres_direct_d";
 34
 35double calculate_l2Norm_solvers(const double *xSol, const double *x, aoclsparse_int n)
 36{
 37    double norm = 0.0f;
 38    //residual = xSol - x
 39    //xSol obtained by initial spmv
 40    //x obtained by ILU-GMRES iterations
 41    for(aoclsparse_int i = 0; i < n; i++)
 42    {
 43        double a = xSol[i] - x[i];
 44        norm += a * a;
 45    }
 46    norm = sqrt(norm);
 47    return norm;
 48}
 49
 50aoclsparse_int monit(aoclsparse_int n,
 51                     const double  *x,
 52                     const double  *r __attribute__((unused)),
 53                     double        *rinfo,
 54                     void          *udata __attribute__((unused)))
 55{
 56    int    it  = (int)rinfo[30];
 57    double tol = PREMATURE_STOP_TOLERANCE;
 58
 59    std::ios oldState(nullptr);
 60    oldState.copyfmt(std::cout);
 61
 62    std::ios_base::fmtflags fmt = std::cout.flags();
 63    fmt |= std::ios_base::scientific | std::ios_base::right | std::ios_base::showpos;
 64
 65    if(!(it % 10))
 66    {
 67        std::cout << std::setw(5) << std::right << " iter"
 68                  << " " << std::setw(16) << std::right << "optim";
 69        for(int i = 0; i < n; i++)
 70            std::cout << std::setw(8) << std::right << "x[" << i << "]";
 71        std::cout << std::endl;
 72    }
 73    std::cout << std::setw(5) << std::right << (int)rinfo[30] << " " << std::setw(16) << std::right
 74              << std::scientific << std::setprecision(8) << rinfo[0];
 75    std::cout << std::setprecision(2) << std::showpos;
 76    for(int i = 0; i < n; i++)
 77        std::cout << " " << x[i];
 78    std::cout << std::endl;
 79    std::cout << std::resetiosflags(fmt);
 80
 81    //reset std::cout state
 82    std::cout.copyfmt(oldState);
 83
 84    if(rinfo[0] < tol) // check for premature stop
 85    {
 86        return 1; // request to interrupt
 87    }
 88    return 0;
 89}
 90
 91int main()
 92{
 93    std::vector<aoclsparse_int> csr_row_ptr;
 94    std::vector<aoclsparse_int> csr_col_ind;
 95    std::vector<double>         csr_val;
 96
 97    int               m, n, nnz;
 98    aoclsparse_status exit_status = aoclsparse_status_success;
 99
100    std::string filename = "cage4.mtx";
101    //symmetry = 0; //real unsymmetric, symmetry = 0%, spd = no
102    //https://www.cise.ufl.edu/research/sparse/MM/vanHeukelum/cage4.tar.gz
103    n   = 9;
104    m   = 9;
105    nnz = 49;
106    csr_row_ptr.assign({0, 5, 10, 15, 20, 26, 32, 38, 44, 49});
107    csr_col_ind.assign({0, 1, 3, 4, 7, 0, 1, 2, 4, 5, 1, 2, 3, 5, 6, 0, 2, 3, 6, 7, 0, 1, 4, 5, 6,
108                        8, 1, 2, 4, 5, 7, 8, 2, 3, 4, 6, 7, 8, 0, 3, 5, 6, 7, 8, 4, 5, 6, 7, 8});
109    csr_val.assign({0.75, 0.14, 0.11, 0.14, 0.11, 0.08, 0.69, 0.11, 0.08, 0.11, 0.09, 0.67, 0.08,
110                    0.09, 0.08, 0.09, 0.14, 0.73, 0.14, 0.09, 0.04, 0.04, 0.54, 0.14, 0.11, 0.25,
111                    0.05, 0.05, 0.08, 0.45, 0.08, 0.15, 0.04, 0.04, 0.09, 0.47, 0.09, 0.18, 0.05,
112                    0.05, 0.14, 0.11, 0.55, 0.25, 0.08, 0.08, 0.09, 0.08, 0.17});
113
114    // Create aocl sparse matrix
115    aoclsparse_matrix     A;
116    aoclsparse_index_base base = aoclsparse_index_base_zero;
117    aoclsparse_mat_descr  descr_a;
118    aoclsparse_operation  trans = aoclsparse_operation_none;
119    aoclsparse_create_dcsr(&A,
120                           base,
121                           (aoclsparse_int)n,
122                           (aoclsparse_int)n,
123                           (aoclsparse_int)nnz,
124                           csr_row_ptr.data(),
125                           csr_col_ind.data(),
126                           csr_val.data());
127    aoclsparse_create_mat_descr(&descr_a);
128    aoclsparse_set_mat_type(descr_a, aoclsparse_matrix_type_symmetric);
129    aoclsparse_set_mat_fill_mode(descr_a, aoclsparse_fill_mode_lower);
130    aoclsparse_set_mat_index_base(descr_a, base);
131    aoclsparse_set_sv_hint(A, trans, descr_a, 100);
132    aoclsparse_optimize(A);
133
134    // Initialize initial point x0 and right hand side b
135    double *expected_sol = NULL;
136    double *x            = NULL;
137    double *b            = NULL;
138    double  norm         = 0.0;
139    double  rinfo[100];
140    double  alpha = 1.0, beta = 0.;
141    int     rs_iters = 7;
142    char    rs_iters_string[16];
143
144    expected_sol = new double[n];
145    x            = new double[n];
146    b            = new double[n];
147
148    double init_x = 1.0, ref_x = 0.5;
149    for(int i = 0; i < n; i++)
150    {
151        expected_sol[i] = ref_x;
152        x[i]            = init_x;
153    }
154
155    aoclsparse_dmv(trans, &alpha, A, descr_a, expected_sol, &beta, b);
156
157    double                  tol    = PREMATURE_STOP_TOLERANCE;
158    aoclsparse_itsol_handle handle = NULL;
159    // create GMRES handle
160    aoclsparse_itsol_d_init(&handle);
161
162    exit_status = aoclsparse_itsol_option_set(handle, "iterative method", "GMRES");
163    if(exit_status != aoclsparse_status_success)
164        printf("Warning an iterative method option could not be set, exit status = %d\n",
165               exit_status);
166
167    exit_status = aoclsparse_itsol_option_set(handle, "gmres preconditioner", "ILU0");
168    if(exit_status != aoclsparse_status_success)
169        printf("Warning gmres preconditioner option could not be set, exit status = %d\n",
170               exit_status);
171
172    exit_status = aoclsparse_itsol_option_set(handle, "gmres iteration limit", "50");
173    if(exit_status != aoclsparse_status_success)
174        printf("Warning gmres iteration limit option could not be set, exit status = %d\n",
175               exit_status);
176
177    sprintf(rs_iters_string, "%d", (int)rs_iters);
178    exit_status = aoclsparse_itsol_option_set(handle, "gmres restart iterations", rs_iters_string);
179    if(exit_status != aoclsparse_status_success)
180        printf("Warning gmres restart iterations option could not be set, exit status = %d\n",
181               exit_status);
182
183    aoclsparse_itsol_handle_prn_options(handle);
184
185    // Call GMRES solver
186    aoclsparse_status status;
187
188    status = aoclsparse_itsol_d_solve(handle, n, A, descr_a, b, x, rinfo, NULL, monit, &n);
189    if(status == aoclsparse_status_success || status == aoclsparse_status_user_stop
190       || aoclsparse_status_maxit == status)
191    {
192        norm = calculate_l2Norm_solvers(expected_sol, x, n);
193        printf("input = %s\n", filename.c_str());
194        printf("solver = %s\n", gmres_sol);
195        printf("no of rows = %d\n", (int)m);
196        printf("no of cols = %d\n", (int)n);
197        printf("no of nnz = %d\n", (int)nnz);
198        printf("monitoring tolerance = %e\n", tol);
199        printf("restart iterations = %d\n", (int)rs_iters);
200        printf("residual achieved = %e\n", rinfo[0]);
201        printf("total iterations = %d\n", (int)rinfo[30]);
202        printf("l2 Norm = %e\n", norm);
203    }
204    else
205    {
206        std::cout << "Something unexpected happened!, Status = " << status << std::endl;
207    }
208
209    delete[] expected_sol;
210    delete[] x;
211    delete[] b;
212    aoclsparse_itsol_destroy(&handle);
213    aoclsparse_destroy_mat_descr(descr_a);
214    aoclsparse_destroy(&A);
215    printf("\n");
216    fflush(stdout);
217
218    return 0;
219}


  1/* ************************************************************************
  2 * Copyright (c) 2022-2025 Advanced Micro Devices, Inc.
  3 *
  4 * Permission is hereby granted, free of charge, to any person obtaining a copy
  5 * of this software and associated documentation files (the "Software"), to deal
  6 * in the Software without restriction, including without limitation the rights
  7 * to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
  8 * copies of the Software, and to permit persons to whom the Software is
  9 * furnished to do so, subject to the following conditions:
 10 *
 11 * The above copyright notice and this permission notice shall be included in
 12 * all copies or substantial portions of the Software.
 13 *
 14 * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
 15 * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
 16 * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
 17 * AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
 18 * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
 19 * OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
 20 * THE SOFTWARE.
 21 *
 22 * ************************************************************************ */
 23
 24#include "aoclsparse.h"
 25
 26#include <iomanip>
 27#include <iostream>
 28#include <math.h>
 29#include <vector>
 30
 31// Define custom log printer
 32void printer(float rinfo[100], bool header)
 33{
 34    std::ios_base::fmtflags fmt = std::cout.flags();
 35    fmt |= std::cout.scientific | std::cout.right;
 36    if(header)
 37        std::cout << std::setw(5) << std::right << " iter"
 38                  << " " << std::setw(16) << std::right << "optim"
 39                  << "  " << std::setw(3) << std::endl;
 40    std::cout << std::setw(5) << std::right << (int)rinfo[30] << " " << std::setw(16) << std::right
 41              << std::scientific << rinfo[0] << "  " << std::endl;
 42    std::cout << std::resetiosflags(fmt);
 43}
 44
 45int main()
 46{
 47    // CSR symmetric matrix. Only the lower triangle is stored
 48    std::vector<aoclsparse_int> icrow, icol;
 49    std::vector<float>          aval;
 50    aoclsparse_int              n = 8, nnz = 18;
 51    icrow.assign({0, 1, 2, 5, 6, 8, 11, 15, 18});
 52    icol.assign({0, 1, 0, 1, 2, 3, 1, 4, 0, 4, 5, 0, 3, 4, 6, 2, 5, 7});
 53    aval.assign({19, 10, 1, 8, 11, 13, 2, 11, 2, 1, 9, 7, 9, 5, 12, 5, 5, 9});
 54
 55    // Create aocl sparse matrix and its descriptor
 56    aoclsparse_matrix     A;
 57    aoclsparse_index_base base = aoclsparse_index_base_zero;
 58    aoclsparse_mat_descr  descr_a;
 59    aoclsparse_operation  trans = aoclsparse_operation_none;
 60    aoclsparse_create_scsr(&A, base, n, n, nnz, icrow.data(), icol.data(), aval.data());
 61    aoclsparse_create_mat_descr(&descr_a);
 62    aoclsparse_set_mat_type(descr_a, aoclsparse_matrix_type_symmetric);
 63    aoclsparse_set_mat_fill_mode(descr_a, aoclsparse_fill_mode_lower);
 64    aoclsparse_set_mat_index_base(descr_a, base);
 65    aoclsparse_set_sv_hint(A, trans, descr_a, 100);
 66    aoclsparse_optimize(A);
 67
 68    // create CG handle
 69    aoclsparse_itsol_handle handle = nullptr;
 70    aoclsparse_itsol_s_init(&handle);
 71
 72    // Change options (update to use )
 73    if(aoclsparse_itsol_option_set(handle, "CG Rel Tolerance", "1.0e-06")
 74           != aoclsparse_status_success
 75       || aoclsparse_itsol_option_set(handle, "CG preconditioner", "user")
 76              != aoclsparse_status_success)
 77        std::cout << "Warning an option could not be set" << std::endl;
 78
 79    // Initialize initial point x0 and right hand side b
 80    std::vector<float> x, b, expected_sol, y;
 81    x.assign({1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0});
 82    b.assign({0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0});
 83    expected_sol.assign({1.E0, 0.E0, 1.E0, 0.E0, 1.E0, 0.E0, 1.E0, 0.E0});
 84    float alpha = 1.0, beta = 0.;
 85    y.resize(n);
 86    aoclsparse_smv(trans, &alpha, A, descr_a, expected_sol.data(), &beta, b.data());
 87
 88    // initialize size and rhs inside the handle
 89    aoclsparse_itsol_s_rci_input(handle, n, b.data());
 90
 91    // Call CG solver
 92    aoclsparse_itsol_rci_job ircomm = aoclsparse_rci_start;
 93    aoclsparse_status        status;
 94    float                   *u = nullptr;
 95    float                   *v = nullptr;
 96    float                    rinfo[100];
 97    float                    tol = 1.0e-5;
 98    bool                     hdr;
 99    std::cout << std::endl;
100    while(ircomm != aoclsparse_rci_stop)
101    {
102        status = aoclsparse_itsol_s_rci_solve(handle, &ircomm, &u, &v, x.data(), rinfo);
103        if(status != aoclsparse_status_success)
104            break;
105        switch(ircomm)
106        {
107        case aoclsparse_rci_mv:
108            // Compute v = Au
109            beta   = 0.0;
110            alpha  = 1.0;
111            status = aoclsparse_smv(trans, &alpha, A, descr_a, u, &beta, v);
112            if(status != aoclsparse_status_success)
113                ircomm = aoclsparse_rci_stop;
114            break;
115
116        case aoclsparse_rci_precond:
117            // apply Symmetric Gauss-Seidel preconditioner step
118            status = aoclsparse_strsv(aoclsparse_operation_none, alpha, A, descr_a, u, y.data());
119            if(status != aoclsparse_status_success)
120            {
121                ircomm = aoclsparse_rci_stop;
122                break;
123            }
124            for(aoclsparse_int i = 0; i < n; i++)
125                y[i] *= aval[icrow[i + 1] - 1];
126            status
127                = aoclsparse_strsv(aoclsparse_operation_transpose, alpha, A, descr_a, y.data(), v);
128            if(status != aoclsparse_status_success)
129                ircomm = aoclsparse_rci_stop;
130            break;
131
132        case aoclsparse_rci_stopping_criterion:
133            // No operations required, can be used to monitor the progress of the solve
134            // or defining a custom stopping criterion
135            // print iteration log
136            hdr = ((int)rinfo[30] % 100) == 0;
137            printer(rinfo, hdr);
138            // request solver to stop if custom criterion is met
139            if(rinfo[0] < tol)
140            {
141                std::cout << "User stop. Final residual: " << rinfo[0] << std::endl;
142                ircomm = aoclsparse_rci_interrupt;
143            }
144            break;
145
146        default:
147            break;
148        }
149    }
150    // Print the final results if the internal stopping criterion or the user defined one were met
151    switch(status)
152    {
153    case aoclsparse_status_user_stop:
154    case aoclsparse_status_success:
155        std::cout.precision(2);
156        std::cout << std::scientific;
157        aoclsparse_itsol_handle_prn_options(handle);
158        std::cout << std::endl
159                  << "Solution found: (residual = " << rinfo[0] << " in " << (int)rinfo[30]
160                  << " iterations)" << std::endl
161                  << "   Final X* = ";
162        for(int i = 0; i < n; i++)
163            std::cout << std::setw(9) << x[i] << " ";
164        std::cout << std::endl;
165        std::cout << "Expected X* = ";
166        for(int i = 0; i < n; i++)
167            std::cout << std::setw(9) << expected_sol[i] << " ";
168        std::cout << std::endl;
169        break;
170
171    case aoclsparse_status_maxit:
172        std::cout << "solve stopped after " << (int)rinfo[30] << " iterations" << std::endl
173                  << "residual = " << rinfo[0] << std::endl;
174        break;
175
176    default:
177        std::cout << "Something unexpected happened!" << std::endl;
178    }
179    aoclsparse_itsol_destroy(&handle);
180    aoclsparse_destroy_mat_descr(descr_a);
181    aoclsparse_destroy(&A);
182
183    return 0;
184}


  1/* ************************************************************************
  2 * Copyright (c) 2022-2023 Advanced Micro Devices, Inc.
  3 *
  4 * Permission is hereby granted, free of charge, to any person obtaining a copy
  5 * of this software and associated documentation files (the "Software"), to deal
  6 * in the Software without restriction, including without limitation the rights
  7 * to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
  8 * copies of the Software, and to permit persons to whom the Software is
  9 * furnished to do so, subject to the following conditions:
 10 *
 11 * The above copyright notice and this permission notice shall be included in
 12 * all copies or substantial portions of the Software.
 13 *
 14 * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
 15 * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
 16 * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
 17 * AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
 18 * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
 19 * OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
 20 * THE SOFTWARE.
 21 *
 22 * ************************************************************************ */
 23
 24#include "aoclsparse.h"
 25
 26#include <cfloat>
 27#include <iomanip>
 28#include <iostream>
 29#include <math.h>
 30#include <vector>
 31
 32#define PREMATURE_STOP_TOLERANCE 1e-4
 33const char gmres_sol[] = "gmres_direct_s";
 34
 35float calculate_l2Norm_solvers(const float *xSol, const float *x, aoclsparse_int n)
 36{
 37    float norm = 0.0f;
 38    //residual = xSol - x
 39    //xSol obtained by initial spmv
 40    //x obtained by ILU-GMRES iterations
 41    for(aoclsparse_int i = 0; i < n; i++)
 42    {
 43        float a = xSol[i] - x[i];
 44        norm += a * a;
 45    }
 46    norm = sqrt(norm);
 47    return norm;
 48}
 49
 50aoclsparse_int monit(aoclsparse_int n,
 51                     const float   *x,
 52                     const float   *r __attribute__((unused)),
 53                     float         *rinfo,
 54                     void          *udata __attribute__((unused)))
 55{
 56    int   it  = (int)rinfo[30];
 57    float tol = PREMATURE_STOP_TOLERANCE;
 58
 59    std::ios oldState(nullptr);
 60    oldState.copyfmt(std::cout);
 61
 62    std::ios_base::fmtflags fmt = std::cout.flags();
 63    fmt |= std::ios_base::scientific | std::ios_base::right | std::ios_base::showpos;
 64
 65    if(!(it % 10))
 66    {
 67        std::cout << std::setw(5) << std::right << " iter"
 68                  << " " << std::setw(16) << std::right << "optim";
 69        for(int i = 0; i < n; i++)
 70            std::cout << std::setw(8) << std::right << "x[" << i << "]";
 71        std::cout << std::endl;
 72    }
 73    std::cout << std::setw(5) << std::right << (int)rinfo[30] << " " << std::setw(16) << std::right
 74              << std::scientific << std::setprecision(8) << rinfo[0];
 75    std::cout << std::setprecision(2) << std::showpos;
 76    for(int i = 0; i < n; i++)
 77        std::cout << " " << x[i];
 78    std::cout << std::endl;
 79    std::cout << std::resetiosflags(fmt);
 80
 81    //reset std::cout state
 82    std::cout.copyfmt(oldState);
 83
 84    if(rinfo[0] < tol) // check for premature stop
 85    {
 86        return 1; // request to interrupt
 87    }
 88    return 0;
 89}
 90
 91int main()
 92{
 93    // CSR symmetric matrix. Only the lower triangle is stored
 94    std::vector<aoclsparse_int> csr_row_ptr;
 95    std::vector<aoclsparse_int> csr_col_ind;
 96    std::vector<float>          csr_val;
 97
 98    int               m, n, nnz;
 99    aoclsparse_status exit_status = aoclsparse_status_success;
100
101    std::string filename = "cage4.mtx";
102    //symmetry = 0; //real unsymmetric, symmetry = 0%, spd = no
103    //https://www.cise.ufl.edu/research/sparse/MM/vanHeukelum/cage4.tar.gz
104    n   = 9;
105    m   = 9;
106    nnz = 49;
107    csr_row_ptr.assign({0, 5, 10, 15, 20, 26, 32, 38, 44, 49});
108    csr_col_ind.assign({0, 1, 3, 4, 7, 0, 1, 2, 4, 5, 1, 2, 3, 5, 6, 0, 2, 3, 6, 7, 0, 1, 4, 5, 6,
109                        8, 1, 2, 4, 5, 7, 8, 2, 3, 4, 6, 7, 8, 0, 3, 5, 6, 7, 8, 4, 5, 6, 7, 8});
110    csr_val.assign({0.75, 0.14, 0.11, 0.14, 0.11, 0.08, 0.69, 0.11, 0.08, 0.11, 0.09, 0.67, 0.08,
111                    0.09, 0.08, 0.09, 0.14, 0.73, 0.14, 0.09, 0.04, 0.04, 0.54, 0.14, 0.11, 0.25,
112                    0.05, 0.05, 0.08, 0.45, 0.08, 0.15, 0.04, 0.04, 0.09, 0.47, 0.09, 0.18, 0.05,
113                    0.05, 0.14, 0.11, 0.55, 0.25, 0.08, 0.08, 0.09, 0.08, 0.17});
114
115    // Create aocl sparse matrix
116    aoclsparse_matrix     A;
117    aoclsparse_index_base base = aoclsparse_index_base_zero;
118    aoclsparse_mat_descr  descr_a;
119    aoclsparse_operation  trans = aoclsparse_operation_none;
120    aoclsparse_create_scsr(&A,
121                           base,
122                           (aoclsparse_int)n,
123                           (aoclsparse_int)n,
124                           (aoclsparse_int)nnz,
125                           csr_row_ptr.data(),
126                           csr_col_ind.data(),
127                           csr_val.data());
128    aoclsparse_create_mat_descr(&descr_a);
129    aoclsparse_set_mat_type(descr_a, aoclsparse_matrix_type_symmetric);
130    aoclsparse_set_mat_fill_mode(descr_a, aoclsparse_fill_mode_lower);
131    aoclsparse_set_mat_index_base(descr_a, base);
132    aoclsparse_set_sv_hint(A, trans, descr_a, 100);
133    aoclsparse_optimize(A);
134
135    // Initialize initial point x0 and right hand side b
136    float *expected_sol = NULL;
137    float *x            = NULL;
138    float *b            = NULL;
139    float  norm         = 0.0;
140    float  rinfo[100];
141    float  alpha = 1.0, beta = 0.;
142    int    rs_iters = 7;
143    char   rs_iters_string[16];
144
145    expected_sol = new float[n];
146    x            = new float[n];
147    b            = new float[n];
148
149    float init_x = 1.0, ref_x = 0.5;
150    for(int i = 0; i < n; i++)
151    {
152        expected_sol[i] = ref_x;
153        x[i]            = init_x;
154    }
155
156    // generate RHS
157    // Reference SPMV CSR implementation
158    aoclsparse_smv(trans, &alpha, A, descr_a, expected_sol, &beta, b);
159
160    float                   tol    = PREMATURE_STOP_TOLERANCE;
161    aoclsparse_itsol_handle handle = nullptr;
162    // create GMRES handle
163    aoclsparse_itsol_s_init(&handle);
164
165    exit_status = aoclsparse_itsol_option_set(handle, "iterative method", "GMRES");
166    if(exit_status != aoclsparse_status_success)
167        printf("Warning an iterative method option could not be set, exit status = %d\n",
168               exit_status);
169
170    exit_status = aoclsparse_itsol_option_set(handle, "gmres preconditioner", "ILU0");
171    if(exit_status != aoclsparse_status_success)
172        printf("Warning gmres preconditioner option could not be set, exit status = %d\n",
173               exit_status);
174
175    exit_status = aoclsparse_itsol_option_set(handle, "gmres iteration limit", "50");
176    if(exit_status != aoclsparse_status_success)
177        printf("Warning gmres iteration limit option could not be set, exit status = %d\n",
178               exit_status);
179
180    sprintf(rs_iters_string, "%d", (int)rs_iters);
181    exit_status = aoclsparse_itsol_option_set(handle, "gmres restart iterations", rs_iters_string);
182    if(exit_status != aoclsparse_status_success)
183        printf("Warning gmres restart iterations option could not be set, exit status = %d\n",
184               exit_status);
185
186    aoclsparse_itsol_handle_prn_options(handle);
187
188    // Call GMRES solver
189    aoclsparse_status status;
190
191    status = aoclsparse_itsol_s_solve(handle, n, A, descr_a, b, x, rinfo, nullptr, monit, &n);
192    if(status == aoclsparse_status_success || status == aoclsparse_status_user_stop
193       || aoclsparse_status_maxit == status)
194    {
195        norm = calculate_l2Norm_solvers(expected_sol, x, n);
196
197        printf("input = %s\n", filename.c_str());
198        printf("solver = %s\n", gmres_sol);
199        printf("no of rows = %d\n", (int)m);
200        printf("no of cols = %d\n", (int)n);
201        printf("no of nnz = %d\n", (int)nnz);
202        printf("monitoring tolerance = %e\n", tol);
203        printf("restart iterations = %d\n", (int)rs_iters);
204        printf("residual achieved = %e\n", rinfo[0]);
205        printf("total iterations = %d\n", (int)rinfo[30]);
206        printf("l2 Norm = %e\n", norm);
207    }
208    else
209    {
210        std::cout << "Something unexpected happened!, Status = " << status << std::endl;
211    }
212
213    delete[] expected_sol;
214    delete[] x;
215    delete[] b;
216    aoclsparse_destroy_mat_descr(descr_a);
217    aoclsparse_destroy(&A);
218    aoclsparse_itsol_destroy(&handle);
219    printf("\n");
220    fflush(stdout);
221
222    return 0;
223}

Note

This function returns control back to the user under certain circumstances. The table in aoclsparse_itsol_rci_job indicates what actions are required to be performed by the user.

Parameters:

  • handle – [inout] problem handle. Needs to be previously initialized by aoclsparse_itsol_s_init or aoclsparse_itsol_d_init or aoclsparse_itsol_c_init or aoclsparse_itsol_z_init and then populated using either aoclsparse_itsol_s_rci_input or aoclsparse_itsol_d_rci_input, or aoclsparse_itsol_c_rci_input or aoclsparse_itsol_z_rci_input as appropriate.

  • ircomm – [inout] pointer to the reverse communication instruction flag and defined in aoclsparse_itsol_rci_job.

  • u – [inout] pointer to a generic vector of data. The solver will point to the data on which the operation defined by ircomm needs to be applied.

  • v – [inout] pointer to a generic vector of data. The solver will ask that the result of the operation defined by ircomm be stored in v.

  • x – [inout] dense vector of unknowns. On input, it should contain the initial guess from which to start the iterative process. If there is no good initial estimate guess then any arbitrary but finite values can be used. On output, it contains an estimate to the solution of the linear system of equations up to the requested tolerance, e.g. see “cg rel tolerance” or “cg abs tolerance” in Options.

  • rinfo – [out] vector containing information and stats related to the iterative solve, see Information Array. This parameter can be used to monitor progress and define a custom stopping criterion when the solver returns control to user with ircomm = aoclsparse_rci_stopping_criterion.