Forward and Reverse Communication Interfaces#

The suite presents two separate interfaces to all the iterative solvers, a direct one, aoclsparse_itsol_d_solve() (aoclsparse_itsol_s_solve()), and a reverse communication (RCI) one aoclsparse_itsol_d_rci_solve() (aoclsparse_itsol_s_rci_solve()). While the underlying algorithms are exactly the same, the difference lies in how data is communicated to the solvers.

The direct communication interface expects to have explicit access to the coefficient matrix \(A\). On the other hand, the reverse communication interface makes no assumption on the matrix storage. Thus when the solver requires some matrix operation such as a matrix-vector product, it returns control to the user and asks the user perform the operation and provide the results by calling again the RCI solver.

Recommended Workflow#

For solving a linear system of equations, the following workflow is recommended:

• Call aoclsparse_itsol_s_init() or aoclsparse_itsol_d_init() to initialize aoclsparse_itsol_handle.

• Choose the solver and adjust its behaviour by setting optional parameters with aoclsparse_itsol_option_set(), see there all options available.

• If the reverse communication interface is desired, define the system’s input with aoclsparse_itsol_s_rci_input() (or aoclsparse_itsol_d_rci_input()).

• Solve the system with either using direct interface aoclsparse_itsol_s_solve() (or aoclsparse_itsol_d_solve()) or reverse communication interface aoclsparse_itsol_s_rci_solve() (or aoclsparse_itsol_d_rci_solve())

• Free the memory with aoclsparse_itsol_destroy().

Information Array#

The array rinfo[100] is used by the solvers (e.g. aoclsparse_itsol_s_solve() or aoclsparse_itsol_d_rci_solve()) to report back useful convergence metrics and other solver statistics. The user callback monit is also equipped with this array and can be used to view or monitor the state of the solver. The solver will populate the following entries with the most recent iteration data

References#

aoclsparse_itsol_rci_job#

enum aoclsparse_itsol_rci_job#

Values of ircomm used by the iterative solver reverse communication interface (RCI) aoclsparse_itsol_d_rci_solve and aoclsparse_itsol_s_rci_solve to communicate back to the user which operation is required.

Values:

enumerator aoclsparse_rci_interrupt#

if set by the user, signals the solver to terminate. This is never set by the solver. Terminate.

enumerator aoclsparse_rci_stop#

found a solution within specified tolerance (see options “cg rel tolerance”, “cg abs tolerance”, “gmres rel tolerance”, and “gmres abs tolerance” in Options). Terminate, vector x contains the solution.

enumerator aoclsparse_rci_start#

initial value of the ircomm flag, no action required. Call solver.

enumerator aoclsparse_rci_mv#

perform the matrix-vector product \( v = Au\). Return control to solver.

enumerator aoclsparse_rci_precond#

perform a preconditioning step on the vector \(u\) and store in \(v\). If the preconditioner \(M\) has explicit matrix form, then applying the preconditioner would result in the operations \( v=Mu \) or \(v=M^{-1}u\). The latter would be performed by solving the linear system of equations \(Mv=u\). Return control to solver.

enumerator aoclsparse_rci_stopping_criterion#

perform a monitoring step and check for custom stopping criteria. If using a positive tolerance value for the convergence options (see aoclsparse_rci_stop), then this step can be ignored and control can be returned to solver.

aoclsparse_itsol_?_init()#

aoclsparse_status aoclsparse_itsol_s_init(aoclsparse_itsol_handle *handle)#

aoclsparse_status aoclsparse_itsol_d_init(aoclsparse_itsol_handle *handle)#

Initialize a problem handle ( aoclsparse_itsol_handle) for the iterative solvers suite of the library.

aoclsparse_itsol_s_init and aoclsparse_itsol_d_init initialize a data structure referred to as problem handle. This handle is used by iterative solvers (itsol) suite to setup options, define which solver to use, etc.

Note

Once the handle is no longer needed, it can be destroyed and the memory released by calling aoclsparse_itsol_destroy.

Parameters

handle – [inout] the pointer to the problem handle data structure.

Return values

• aoclsparse_status_success – the operation completed successfully.

• aoclsparse_status_memory_error – internal memory allocation error.

• aoclsparse_status_invalid_pointer – the pointer to the problem handle is invalid.

• aoclsparse_status_internal_error – an unexpected error occurred.

aoclsparse_itsol_destroy()#

void aoclsparse_itsol_destroy(aoclsparse_itsol_handle *handle)#

Free the memory reserved in a problem handle previously initialized by aoclsparse_itsol_s_init or aoclsparse_itsol_d_init.

Once the problem handle is no longer needed, calling this function to deallocate the memory is advisable to avoid memory leaks.

Note

Passing a handle that has not been initialized by aoclsparse_itsol_s_init or aoclsparse_itsol_d_init may have unpredictable results.

Parameters

handle – [inout] pointer to a problem handle.

aoclsparse_itsol_?_solve()#

aoclsparse_status aoclsparse_itsol_s_solve(aoclsparse_itsol_handle handle, aoclsparse_int n, aoclsparse_matrix mat, const aoclsparse_mat_descr descr, const float *b, float *x, float rinfo[100], aoclsparse_int precond(aoclsparse_int flag, aoclsparse_int n, const float *u, float *v, void *udata), aoclsparse_int monit(aoclsparse_int n, const float *x, const float *r, float rinfo[100], void *udata), void *udata)#

aoclsparse_status aoclsparse_itsol_d_solve(aoclsparse_itsol_handle handle, aoclsparse_int n, aoclsparse_matrix mat, const aoclsparse_mat_descr descr, const double *b, double *x, double rinfo[100], aoclsparse_int precond(aoclsparse_int flag, aoclsparse_int n, const double *u, double *v, void *udata), aoclsparse_int monit(aoclsparse_int n, const double *x, const double *r, double rinfo[100], void *udata), void *udata)#

Forward communication interface to the iterative solvers suite of the library.

This function solves the linear system of equations

\[ Ax=b, \]

where the matrix of coefficients \(A\) is defined by mat. 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 expected workflow is as follows:

1. Call aoclsparse_itsol_s_init or aoclsparse_itsol_d_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. Solve the system by calling aoclsparse_itsol_s_solve or aoclsparse_itsol_d_solve.

4. 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.

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

6. Free the memory with aoclsparse_itsol_destroy.

This interface requires to explicitly provide the matrix \(A\) and its descriptor descr, this kind of interface is also known as _forward communication_ which contrasts with *reverse communication* in which case the matrix \(A\) and its descriptor descr need not be explicitly available. For more details on the latter, see aoclsparse_itsol_d_rci_solve or aoclsparse_itsol_s_rci_solve.

Example - CG / floating point double precision (tests/examples/sample_itsol_d_cg.cpp)

  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 <iomanip>
 27#include <iostream>
 28#include <math.h>
 29#include <vector>
 30
 31aoclsparse_int monit(aoclsparse_int n,
 32                     const double  *x,
 33                     const double  *r __attribute__((unused)),
 34                     double        *rinfo,
 35                     void          *udata __attribute((unused)))
 36{
 37    int                     it  = (int)rinfo[30];
 38    std::ios_base::fmtflags fmt = std::cout.flags();
 39    fmt |= std::ios_base::scientific | std::ios_base::right | std::ios_base::showpos;
 40    if(!(it % 10))
 41    {
 42        std::cout << std::setw(5) << std::right << " iter"
 43                  << " " << std::setw(16) << std::right << "optim";
 44        for(int i = 0; i < n; i++)
 45            std::cout << std::setw(8) << std::right << "x[" << i << "]";
 46        std::cout << std::endl;
 47    }
 48    std::cout << std::setw(5) << std::right << (int)rinfo[30] << " " << std::setw(16) << std::right
 49              << std::scientific << std::setprecision(8) << rinfo[0];
 50    std::cout << std::setprecision(2) << std::showpos;
 51    for(int i = 0; i < n; i++)
 52        std::cout << " " << x[i];
 53    std::cout << std::endl;
 54    std::cout << std::resetiosflags(fmt);
 55    if(rinfo[0] < 1.0e-12) // check for premature stop
 56        return 1; // request to interrupt
 57    return 0;
 58}
 59
 60int main()
 61{
 62    // CSR symmetric matrix. Only the lower triangle is stored
 63    std::vector<aoclsparse_int> icrow, icol;
 64    std::vector<double>         aval;
 65    aoclsparse_int              n = 8, nnz = 18;
 66    icrow.assign({0, 1, 2, 5, 6, 8, 11, 15, 18});
 67    icol.assign({0, 1, 0, 1, 2, 3, 1, 4, 0, 4, 5, 0, 3, 4, 6, 2, 5, 7});
 68    aval.assign({19, 10, 1, 8, 11, 13, 2, 11, 2, 1, 9, 7, 9, 5, 12, 5, 5, 9});
 69
 70    // Create aocl sparse matrix
 71    aoclsparse_matrix     A;
 72    aoclsparse_index_base base = aoclsparse_index_base_zero;
 73    aoclsparse_mat_descr  descr_a;
 74    aoclsparse_operation  trans = aoclsparse_operation_none;
 75    aoclsparse_create_dcsr(&A, base, n, n, nnz, icrow.data(), icol.data(), aval.data());
 76    aoclsparse_create_mat_descr(&descr_a);
 77    aoclsparse_set_mat_type(descr_a, aoclsparse_matrix_type_symmetric);
 78    aoclsparse_set_mat_fill_mode(descr_a, aoclsparse_fill_mode_lower);
 79    aoclsparse_set_mat_index_base(descr_a, base);
 80    aoclsparse_set_sv_hint(A, trans, descr_a, 100);
 81    aoclsparse_optimize(A);
 82
 83    // Initialize initial point x0 and right hand side b
 84    std::vector<double> x, b, expected_sol, y;
 85    x.assign({1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0});
 86    b.assign({0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0});
 87    expected_sol.assign({1.E0, 0.E0, 1.E0, 0.E0, 1.E0, 0.E0, 1.E0, 0.E0});
 88    double alpha = 1.0, beta = 0.;
 89    y.resize(n);
 90    aoclsparse_dmv(trans, &alpha, A, descr_a, expected_sol.data(), &beta, b.data());
 91
 92    // create CG handle
 93    aoclsparse_itsol_handle handle = nullptr;
 94    aoclsparse_itsol_d_init(&handle);
 95
 96    if(aoclsparse_itsol_option_set(handle, "CG Abs Tolerance", "5.0e-6")
 97           != aoclsparse_status_success
 98       || aoclsparse_itsol_option_set(handle, "CG Preconditioner", "SGS")
 99              != aoclsparse_status_success)
100        std::cout << "Warning an option could not be set" << std::endl;
101
102    // Call CG solver
103    aoclsparse_status status;
104    double            rinfo[100];
105    status = aoclsparse_itsol_d_solve(
106        handle, n, A, descr_a, b.data(), x.data(), rinfo, nullptr, monit, &n);
107    if(status == aoclsparse_status_success)
108    {
109        std::cout.precision(2);
110        std::cout << std::scientific;
111        aoclsparse_itsol_handle_prn_options(handle);
112        std::cout << std::endl
113                  << "Solution found: (residual = " << rinfo[0] << " in " << (int)rinfo[30]
114                  << " iterations)" << std::endl
115                  << "   Final X* = ";
116        for(int i = 0; i < n; i++)
117            std::cout << std::setw(9) << x[i] << " ";
118        std::cout << std::endl;
119        std::cout << "Expected X* = ";
120        for(int i = 0; i < n; i++)
121            std::cout << std::setw(9) << expected_sol[i] << " ";
122        std::cout << std::endl;
123    }
124    else if(status == aoclsparse_status_user_stop)
125    {
126        std::cout << "User requested to terminate at interation " << (aoclsparse_int)rinfo[30]
127                  << std::endl;
128    }
129    else
130        std::cout << "Something unexpected happened. Status = " << status << std::endl;
131
132    aoclsparse_itsol_destroy(&handle);
133    aoclsparse_destroy_mat_descr(descr_a);
134    aoclsparse_destroy(&A);
135
136    return 0;
137}

Example - GMRES / floating point double precision (tests/examples/sample_itsol_d_gmres.cpp)

  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}

Example - CG / floating point single precision (tests/examples/sample_itsol_s_cg.cpp)

  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 <iomanip>
 27#include <iostream>
 28#include <math.h>
 29#include <vector>
 30
 31aoclsparse_int monit(aoclsparse_int n,
 32                     const float   *x,
 33                     const float   *r __attribute__((unused)),
 34                     float         *rinfo,
 35                     void          *udata __attribute__((unused)))
 36{
 37    int                     it  = (int)rinfo[30];
 38    std::ios_base::fmtflags fmt = std::cout.flags();
 39    fmt |= std::ios_base::scientific | std::ios_base::right | std::ios_base::showpos;
 40    if(!(it % 10))
 41    {
 42        std::cout << std::setw(5) << std::right << " iter"
 43                  << " " << std::setw(16) << std::right << "residual norm2";
 44        for(int i = 0; i < n; i++)
 45            std::cout << std::setw(8) << std::right << "x[" << i << "]";
 46        std::cout << std::endl;
 47    }
 48    std::cout << std::setw(5) << std::right << (int)rinfo[30] << " " << std::setw(16) << std::right
 49              << std::scientific << std::setprecision(8) << rinfo[0];
 50    std::cout << std::setprecision(2) << std::showpos;
 51    for(int i = 0; i < n; i++)
 52        std::cout << " " << x[i];
 53    std::cout << std::endl;
 54    std::cout << std::resetiosflags(fmt);
 55    if(rinfo[0] < 1.0e-12) // check for premature stop
 56        return 1; // request to interrupt
 57    return 0;
 58}
 59
 60int main()
 61{
 62    // CSR symmetric matrix. Only the lower triangle is stored
 63    std::vector<aoclsparse_int> icrow, icol;
 64    std::vector<float>          aval;
 65    aoclsparse_int              n = 8, nnz = 18;
 66    icrow.assign({0, 1, 2, 5, 6, 8, 11, 15, 18});
 67    icol.assign({0, 1, 0, 1, 2, 3, 1, 4, 0, 4, 5, 0, 3, 4, 6, 2, 5, 7});
 68    aval.assign({19, 10, 1, 8, 11, 13, 2, 11, 2, 1, 9, 7, 9, 5, 12, 5, 5, 9});
 69
 70    // Create aocl sparse matrix
 71    aoclsparse_matrix     A;
 72    aoclsparse_index_base base = aoclsparse_index_base_zero;
 73    aoclsparse_mat_descr  descr_a;
 74    aoclsparse_operation  trans = aoclsparse_operation_none;
 75    aoclsparse_create_scsr(&A, base, n, n, nnz, icrow.data(), icol.data(), aval.data());
 76    aoclsparse_create_mat_descr(&descr_a);
 77    aoclsparse_set_mat_type(descr_a, aoclsparse_matrix_type_symmetric);
 78    aoclsparse_set_mat_fill_mode(descr_a, aoclsparse_fill_mode_lower);
 79    aoclsparse_set_mat_index_base(descr_a, base);
 80    aoclsparse_set_sv_hint(A, trans, descr_a, 100);
 81    aoclsparse_optimize(A);
 82
 83    // Initialize initial point x0 and right hand side b
 84    std::vector<float> x, b, expected_sol, y;
 85    x.assign({1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0});
 86    b.assign({0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0});
 87    expected_sol.assign({1.E0, 0.E0, 1.E0, 0.E0, 1.E0, 0.E0, 1.E0, 0.E0});
 88    float alpha = 1.0, beta = 0.;
 89    y.resize(n);
 90    aoclsparse_smv(trans, &alpha, A, descr_a, expected_sol.data(), &beta, b.data());
 91
 92    // create CG handle
 93    aoclsparse_itsol_handle handle = nullptr;
 94    aoclsparse_itsol_s_init(&handle);
 95
 96    if(aoclsparse_itsol_option_set(handle, "CG Rel Tolerance", "1.0e-06")
 97           != aoclsparse_status_success
 98       || aoclsparse_itsol_option_set(handle, "CG Preconditioner", "SGS")
 99              != aoclsparse_status_success)
100        std::cout << "Warning an option could not be set" << std::endl;
101
102    // Call CG solver
103    aoclsparse_status status;
104    float             rinfo[100];
105    status = aoclsparse_itsol_s_solve(
106        handle, n, A, descr_a, b.data(), x.data(), rinfo, nullptr, monit, &n);
107    if(status == aoclsparse_status_success)
108    {
109        std::cout.precision(2);
110        std::cout << std::scientific;
111        aoclsparse_itsol_handle_prn_options(handle);
112        std::cout << std::endl
113                  << "Solution found: (residual = " << rinfo[0] << " in " << (int)rinfo[30]
114                  << " iterations)" << std::endl
115                  << "   Final X* = ";
116        for(int i = 0; i < n; i++)
117            std::cout << std::setw(9) << x[i] << " ";
118        std::cout << std::endl;
119        std::cout << "Expected X* = ";
120        for(int i = 0; i < n; i++)
121            std::cout << std::setw(9) << expected_sol[i] << " ";
122        std::cout << std::endl;
123    }
124    else
125        std::cout << "Something unexpected happened " << status << std::endl;
126
127    aoclsparse_itsol_destroy(&handle);
128    aoclsparse_destroy_mat_descr(descr_a);
129    aoclsparse_destroy(&A);
130
131    return 0;
132}

Example - GMRES / floating point single precision (tests/examples/sample_itsol_s_gmres.cpp)

  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}

Parameters

• handle – [inout] a valid problem handle, previously initialized by calling aoclsparse_itsol_s_init or aoclsparse_itsol_d_init.

• n – [in] the size of the square matrix mat.

• mat – [inout] coefficient matrix \(A\).

• descr – [inout] matrix descriptor for mat.

• b – [in] right-hand side dense vector \(b\).

• 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.

• precond – [in] (optional, can be nullptr) function pointer to a user routine that applies the preconditioning step

  \[ v = Mu \text{or } v = M^{-1}u,\]
  where \(v\) is the resulting vector of applying a preconditioning step on the vector \(u\) and \(M\) refers to the user specified preconditioner in matrix form and need not be explicitly available. The void pointer udata, is a convenience pointer that can be used by the user to point to user data and is not used by the itsol framework. If the user requests to use a predefined preconditioner already available in the suite (refer to e.g. “cg preconditioner” or “gmres preconditioner” in Options), then this parameter need not be provided.

• monit – [in] (optional, can be nullptr) function pointer to a user monitoring routine. If provided, then at each iteration, the routine is called and can be used to define a custom stopping criteria or to oversee the convergence process. In general, this function need not be provided. If provided then the solver provides n the problem size, x the current iterate, r the current residual vector ( \(r = Ax-b\)), rinfo the current solver’s stats, see Information Array, and udata a convenience pointer that can be used by the user to point to arbitrary user data and is not used by the itsol framework.

• udata – [inout] (optional, can be nullptr) user convenience pointer, it can be used by the user to pass a pointer to user data. It is not modified by the solver.