aoclsparse_ilu_?smoother() - 5.2 English - 68552

  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}