aoclsparse_?symgs() - 5.2 English - 68552

AOCL API Guide (68552)
Document ID
68552
Release Date
2025-12-29
Version
5.2 English
aoclsparse_status aoclsparse_ssymgs(aoclsparse_operation trans, aoclsparse_matrix A, const aoclsparse_mat_descr descr, const float alpha, const float *b, float *x)#
aoclsparse_status aoclsparse_dsymgs(aoclsparse_operation trans, aoclsparse_matrix A, const aoclsparse_mat_descr descr, const double alpha, const double *b, double *x)#
aoclsparse_status aoclsparse_csymgs(aoclsparse_operation trans, aoclsparse_matrix A, const aoclsparse_mat_descr descr, const aoclsparse_float_complex alpha, const aoclsparse_float_complex *b, aoclsparse_float_complex *x)#
aoclsparse_status aoclsparse_zsymgs(aoclsparse_operation trans, aoclsparse_matrix A, const aoclsparse_mat_descr descr, const aoclsparse_double_complex alpha, const aoclsparse_double_complex *b, aoclsparse_double_complex *x)#

Symmetric Gauss Seidel(SYMGS) Preconditioner for real/complex single and double data precisions.

aoclsparse_?symgs performs an iteration of Gauss Seidel preconditioning. Krylov methods such as CG (Conjugate Gradient) and GMRES (Generalized Minimal Residual) are used to solve large sparse linear systems of the form

\[ op(A)\; x = \alpha b, \]
where \(A\) is a sparse matrix of size \(m\), \(op()\) is a linear operator, \(b\) is a dense right-hand side vector and \(x\) is the unknown dense vector, while \(\alpha\) is a scalar. This Gauss Seidel(GS) relaxation is typically used either as a preconditioner for a Krylov solver directly, or as a smoother in a V –cycle of a multigrid preconditioner to accelerate the convergence rate. The Symmetric Gauss Seidel algorithm performs a forward sweep followed by a backward sweep to maintain symmetry of the matrix operation.

To solve a linear system \(Ax = b\), Gauss Seidel(GS) iteration is based on the matrix splitting

\[ A = L + D + U = -E + D - F \]
where \(-E\) or \(L\) is strictly lower triangle, \(D\) is diagonal and \(-F\) or \(D\) is strictly upper triangle. Gauss-Seidel is best derived as element-wise (refer Yousef Saad’s book Iterative Methods for Sparse Linear Systems, Second Edition, Chapter 4.1, p. 125 onwards):
\[ x_i = \frac{1}{a_{ii}} \; \left (b_i - \sum_{j=1}^{i-1} a_{ij} \; x_j - \sum_{j=i+1}^n a_{ij} \; x_j\right) \]
where the first sum is lower triangle i.e., \(-Ex\) and the second sum is upper triangle i.e., \(-Fx\). If we iterate through the rows i=1 to n and keep overwriting/reusing the new \(x_{i}\), we get forward GS, expressed in matrix form as,
\[ (D-E) \; x_{k+1} = F \; x_k + b \]
Iterating through the rows in reverse order from i=n to 1, the upper triangle keeps using the new \(x_{k+1}\) elements and we get backward GS, expressed in matrix form as,
\[ (D-F) \; x_{k+1} = E \; x_k + b \]
The above two equations can be expressed in terms of L, D and U as follows,
\[ (L + D)\;x_1 = b - U\;x_0 \]
\[ (U + D)\;x = b - L\;x_1 \]
So, Symmetric Gauss Seidel (SYMGS) can be computed using two aoclsparse_?mv and two aoclsparse_?trsv operations.

The sparse matrix \(A\) can be either a symmetric or a Hermitian matrix, whose fill is indicated by fill_mode from the matrix descriptor descr where either upper or lower triangular portion of the matrix is used. Matrix \(A\) must be of full rank, that is, the matrix must be invertible. The linear operator \(op()\) can define the transposition or conjugate transposition operations. By default, no transposition is performed. The right-hand-side vector \(b\) and the solution vector \(x\) are dense and must be of the correct size, that is \(m\). If used as fixed point iterative method, the convergence is guaranteed for strictly diagonally dominant and symmetric positive definite matrices from any starting point, x0. However, the API can be applied to wider types of input or as a preconditioning step. Refer Yousef Saad’s Iterative Methods for Sparse Linear Systems 2nd Edition, Theorem 4.9 and related literature for mathematical theory.

1. If the matrix descriptor descr specifies that the matrix \(A\) is to be regarded as having a unitary diagonal, then the main diagonal entries of matrix \(A\) are not accessed and are considered to all be ones.

2. If the matrix \(A\) is described as upper triangular, then only the upper triangular portion of the matrix is referenced. Conversely, if the matrix \(A\) is described lower triangular, then only the lower triangular portion of the matrix is used.

3. This set of APIs allocates couple of work array buffers of size \(m\) for to store intermediate results

4. If the input matrix is of triangular type, the SGS is computed using a single aoclsparse_?trsv operation and a quick return is made without going through the 3-step reference(described above)

  1/* ************************************************************************
  2 * Copyright (c) 2024 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 <cmath>
 27#include <iomanip>
 28#include <iostream>
 29#include <limits>
 30#include <vector>
 31
 32int main(void)
 33{
 34    std::cout << "------------------------------------------------" << std::endl
 35              << " Symmetric Gauss Seidel Preconditioner" << std::endl
 36              << " sample program for real data type" << std::endl
 37              << "------------------------------------------------" << std::endl
 38              << std::endl;
 39
 40    /*
 41     * This example illustrates how to use Symmetric Gauss Seidel API with an iterative solver,
 42        here we have simulated that effect using a loop. The real matrix used is
 43        A = [
 44                10   -1    2    0;
 45                -1   11   -1    3;
 46                 2   -1   10   -1;
 47                 0    3   -1    8;
 48            ]
 49     */
 50
 51    // Create a tri-diagonal matrix in CSR format
 52    const aoclsparse_int n = 4, m = 4, nnz = 14;
 53    aoclsparse_int       icrow[5] = {0, 3, 7, 11, 14};
 54    aoclsparse_int       icol[14] = {0, 1, 2, 0, 1, 2, 3, 0, 1, 2, 3, 1, 2, 3};
 55    double               aval[14] = {10, -1, 2, -1, 11, -1, 3, 2, -1, 10, -1, 3, -1, 8};
 56    aoclsparse_status    status;
 57    bool                 oki, ok;
 58
 59    const double macheps      = std::numeric_limits<double>::epsilon();
 60    const double safe_macheps = (double)2.0 * macheps;
 61    double       exp_tol      = 20.0 * sqrt(safe_macheps);
 62
 63    // Create aoclsparse matrix and its descriptor
 64    aoclsparse_matrix     A;
 65    aoclsparse_index_base base = aoclsparse_index_base_zero;
 66    aoclsparse_mat_descr  descr_a;
 67    aoclsparse_operation  trans;
 68    status = aoclsparse_create_dcsr(&A, base, m, n, nnz, icrow, icol, aval);
 69    if(status != aoclsparse_status_success)
 70    {
 71        std::cerr << "Error returned from aoclsparse_create_dcsr, status = " << status << "."
 72                  << std::endl;
 73        return 3;
 74    }
 75    aoclsparse_create_mat_descr(&descr_a);
 76
 77    double alpha = 1.0;
 78    double x[m]  = {1, 1, 1, 1};
 79
 80    double              b[m] = {14, 30, 26, 35};
 81    std::vector<double> xref;
 82    xref.assign({1.0, 2.0, 3.0, 4.0});
 83
 84    // Indicate to access the lower part of matrix A
 85    aoclsparse_set_mat_type(descr_a, aoclsparse_matrix_type_symmetric);
 86    aoclsparse_set_mat_fill_mode(descr_a, aoclsparse_fill_mode_lower);
 87    trans = aoclsparse_operation_none;
 88
 89    status = aoclsparse_set_symgs_hint(A, trans, descr_a, 1);
 90    if(status != aoclsparse_status_success)
 91    {
 92        std::cerr << "Error returned from aoclsparse_set_symgs_hint, status = " << status << "."
 93                  << std::endl;
 94        return 3;
 95    }
 96    status = aoclsparse_optimize(A);
 97    if(status != aoclsparse_status_success)
 98    {
 99        std::cerr << "Error returned from aoclsparse_optimize, status = " << status << "."
100                  << std::endl;
101        return 3;
102    }
103    // Solve till convergence
104    std::cout << std::setw(12) << " iter";
105    for(int i = 0; i < n; i++)
106        std::cout << std::setw(15) << "x[" << i << "]";
107    std::cout << std::endl;
108    for(int i = 0; i < 8; i++)
109    {
110        status = aoclsparse_dsymgs(trans, A, descr_a, alpha, b, x);
111        if(status != aoclsparse_status_success)
112        {
113            std::cerr << "Error returned from aoclsparse_dsymgs, status = " << status << "."
114                      << std::endl;
115            return 3;
116        }
117        std::cout << std::setw(12) << (int)i << " ";
118        for(int j = 0; j < m; ++j)
119        {
120            std::cout << std::setw(16) << std::scientific << x[j] << " ";
121        }
122        std::cout << std::endl;
123    }
124    // Print and check the result
125    std::cout << "Solving A x = alpha b" << std::endl;
126    std::cout << "  where x and b are dense vectors" << std::endl;
127    std::cout << "  Solution found after 8 iterations, x = " << std::endl;
128    std::cout << std::fixed;
129    std::cout.precision(6);
130    ok = true;
131    for(int i = 0; i < m; ++i)
132    {
133        oki = std::abs(x[i] - xref[i]) <= exp_tol;
134        std::cout << x[i] << (oki ? "  " : "! ");
135        ok &= oki;
136    }
137    std::cout << std::endl;
138
139    // Destroy the aoclsparse memory
140    aoclsparse_destroy_mat_descr(descr_a);
141    aoclsparse_destroy(&A);
142
143    return ok ? 0 : 5;
144}


  1/* ************************************************************************
  2 * Copyright (c) 2024 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 <cmath>
 27#include <complex>
 28#include <iomanip>
 29#include <iostream>
 30#include <limits>
 31#include <vector>
 32
 33int main(void)
 34{
 35    std::cout << "------------------------------------------------" << std::endl
 36              << " Symmetric Gauss Seidel Preconditioner" << std::endl
 37              << " sample program for complex data type" << std::endl
 38              << "------------------------------------------------" << std::endl
 39              << std::endl;
 40
 41    /*
 42     * This example illustrates how to use Symmetric Gauss Seidel API for complex data
 43        with an iterative solver, here we have simulated that effect using a loop.
 44        The complex matrix used is
 45        A = [
 46                   10+i  -1-0.1i    2+0.2i        0;
 47                -1-0.1i  11+1.1i   -1-0.1i   3+0.3i;
 48                 2+0.2i  -1-0.1i      10+i  -1-0.1i;
 49                      0   3+0.3i   -1-0.1i   8+0.8i;
 50            ]
 51     */
 52
 53    // Create a tri-diagonal matrix in CSR format
 54    const aoclsparse_int                   n = 4, m = 4, nnz = 14;
 55    aoclsparse_int                         icrow[5] = {0, 3, 7, 11, 14};
 56    aoclsparse_int                         icol[14] = {0, 1, 2, 0, 1, 2, 3, 0, 1, 2, 3, 1, 2, 3};
 57    std::vector<aoclsparse_double_complex> aval;
 58    // clang-format off
 59    aval.assign({{10,1}, {-1,-0.1}, {2,0.2}, {-1,-0.1},
 60                 {11,1.1}, {-1,-0.1}, {3,0.3}, {2,0.2},
 61                 {-1,-0.1}, {10,1}, {-1,-0.1}, {3, 0.3},
 62                 {-1,-0.1}, {8, 0.8}});
 63    // clang-format on
 64    aoclsparse_status status;
 65    bool              oki, ok;
 66
 67    const double macheps      = std::numeric_limits<double>::epsilon();
 68    const double safe_macheps = (double)2.0 * macheps;
 69    double       exp_tol      = 20.0 * sqrt(safe_macheps);
 70
 71    // Create aoclsparse matrix and its descriptor
 72    aoclsparse_matrix     A;
 73    aoclsparse_index_base base = aoclsparse_index_base_zero;
 74    aoclsparse_mat_descr  descr_a;
 75    aoclsparse_operation  trans;
 76    status = aoclsparse_create_zcsr(&A, base, m, n, nnz, icrow, icol, &aval[0]);
 77    if(status != aoclsparse_status_success)
 78    {
 79        std::cerr << "Error returned from aoclsparse_create_dcsr, status = " << status << "."
 80                  << std::endl;
 81        return 3;
 82    }
 83    aoclsparse_create_mat_descr(&descr_a);
 84
 85    aoclsparse_double_complex alpha = {1.0, 0.};
 86
 87    std::vector<aoclsparse_double_complex> x, b, xref;
 88    x.assign({{1, 0.1}, {1, 0.1}, {1, 0.1}, {1, 0.1}});
 89    b.assign({{13.859999999999999, 2.7999999999999998},
 90              {29.699999999999999, 6},
 91              {25.739999999999998, 5.2000000000000002},
 92              {34.649999999999999, 7}});
 93    xref.assign({{0.99999999992700483, 0.099999999992700442},
 94                 {2.0000000217847944, 0.20000000217847949},
 95                 {3.0000000112573724, 0.30000000112573733},
 96                 {3.9999999237056727, 0.39999999237056733}});
 97
 98    // Indicate to access the lower part of matrix A
 99    aoclsparse_set_mat_type(descr_a, aoclsparse_matrix_type_symmetric);
100    aoclsparse_set_mat_fill_mode(descr_a, aoclsparse_fill_mode_lower);
101    trans = aoclsparse_operation_none;
102
103    status = aoclsparse_set_symgs_hint(A, trans, descr_a, 1);
104    if(status != aoclsparse_status_success)
105    {
106        std::cerr << "Error returned from aoclsparse_set_symgs_hint, status = " << status << "."
107                  << std::endl;
108        return 3;
109    }
110    status = aoclsparse_optimize(A);
111    if(status != aoclsparse_status_success)
112    {
113        std::cerr << "Error returned from aoclsparse_optimize, status = " << status << "."
114                  << std::endl;
115        return 3;
116    }
117    // Solve till convergence
118    std::cout << std::setw(12) << " iter";
119    for(int i = 0; i < n; i++)
120        std::cout << std::setw(33) << "x[" << i << "]";
121    std::cout << std::endl;
122    for(int i = 0; i < 8; i++)
123    {
124        status = aoclsparse_zsymgs(trans, A, descr_a, alpha, &b[0], &x[0]);
125        if(status != aoclsparse_status_success)
126        {
127            std::cerr << "Error returned from aoclsparse_zsymgs, status = " << status << "."
128                      << std::endl;
129            return 3;
130        }
131        std::cout << std::setw(12) << (int)i << " ";
132        for(int j = 0; j < m; ++j)
133        {
134            std::cout << std::setw(8) << std::scientific << std::setprecision(6) << "(" << x[j].real
135                      << ", " << x[j].imag << ")";
136        }
137        std::cout << std::endl;
138    }
139    // Print and check the result
140    std::cout << "Solving A x = alpha b" << std::endl;
141    std::cout << "  where x and b are dense vectors" << std::endl;
142    std::cout << "  Solution found after 8 iterations, x = " << std::endl;
143    std::cout << std::fixed;
144    std::cout.precision(6);
145    ok = true;
146    for(int i = 0; i < m; ++i)
147    {
148        oki = std::abs(x[i].real - xref[i].real) <= exp_tol
149              && std::abs(x[i].imag - xref[i].imag) <= exp_tol;
150        std::cout << "(" << x[i].real << ", " << x[i].imag << "i) " << (oki ? "  " : "! ");
151        ok &= oki;
152    }
153    std::cout << std::endl;
154
155    // Destroy the aoclsparse memory
156    aoclsparse_destroy_mat_descr(descr_a);
157    aoclsparse_destroy(&A);
158
159    return ok ? 0 : 5;
160}

Note

Parameters:
Return values:

  • aoclsparse_status_success – indicates that the operation completed successfully.

  • aoclsparse_status_invalid_size – informs that either m, n or nnz is invalid. The error code also informs if the given sparse matrix \(A\) is not square.

  • aoclsparse_status_invalid_value – informs that either base, trans, matrix type descr->type or fill mode descr->fill_mode is invalid. If the sparse matrix \(A\) is not of full rank, the error code is returned to indicate that the linear system cannot be solved.

  • aoclsparse_status_invalid_pointer – informs that either descr, A, b, or x pointer is invalid.

  • aoclsparse_status_not_implemented – this error occurs when the provided matrix’s aoclsparse_fill_mode is aoclsparse_diag_type_unit or the input format is not aoclsparse_csr_mat, or when aoclsparse_matrix_type is aoclsparse_matrix_type_general and trans is aoclsparse_operation_conjugate_transpose.

aoclsparse_status aoclsparse_ssymgs_mv(aoclsparse_operation trans, aoclsparse_matrix A, const aoclsparse_mat_descr descr, const float alpha, const float *b, float *x, float *y)#
aoclsparse_status aoclsparse_dsymgs_mv(aoclsparse_operation trans, aoclsparse_matrix A, const aoclsparse_mat_descr descr, const double alpha, const double *b, double *x, double *y)#
aoclsparse_status aoclsparse_csymgs_mv(aoclsparse_operation trans, aoclsparse_matrix A, const aoclsparse_mat_descr descr, const aoclsparse_float_complex alpha, const aoclsparse_float_complex *b, aoclsparse_float_complex *x, aoclsparse_float_complex *y)#
aoclsparse_status aoclsparse_zsymgs_mv(aoclsparse_operation trans, aoclsparse_matrix A, const aoclsparse_mat_descr descr, const aoclsparse_double_complex alpha, const aoclsparse_double_complex *b, aoclsparse_double_complex *x, aoclsparse_double_complex *y)#

Symmetric Gauss Seidel Preconditioner followed by SPMV for single and double precision datatypes.

For full details refer to aoclsparse_?symgs().

This variation of SYMGS, namely with a suffix of _mv, performs matrix-vector multiplication between the sparse matrix \(A\) and the Gauss Seidel solution vector \(x\).

  1/* ************************************************************************
  2 * Copyright (c) 2024 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 <cmath>
 27#include <iomanip>
 28#include <iostream>
 29#include <limits>
 30#include <vector>
 31
 32int main(void)
 33{
 34    std::cout
 35        << "------------------------------------------------" << std::endl
 36        << " Symmetric Gauss Seidel Preconditioner followed by sparse matrix-vector multiplication"
 37        << std::endl
 38        << " sample program for real data type" << std::endl
 39        << "------------------------------------------------" << std::endl
 40        << std::endl;
 41
 42    /*
 43     * This example illustrates how to use Symmetric Gauss Seidel API to solve
 44     * linear system of equations. The real matrix used is
 45     *     | 111  2  0  0 |
 46     * A = | 2  111  2  0 |
 47     *     | 0  2  111  2 |
 48     *     | 0  0  2  111 |
 49     */
 50
 51    // Create a tri-diagonal matrix in CSR format
 52    const aoclsparse_int n = 4, m = 4, nnz = 10;
 53    aoclsparse_int       icrow[5] = {0, 2, 5, 8, 10};
 54    aoclsparse_int       icol[10] = {0, 1, 0, 1, 2, 1, 2, 3, 2, 3};
 55    double               aval[10] = {111, 2, 2, 111, 2, 2, 111, 2, 2, 111};
 56    aoclsparse_status    status;
 57    bool                 oki, ok;
 58
 59    const double macheps      = std::numeric_limits<double>::epsilon();
 60    const double safe_macheps = (double)2.0 * macheps;
 61    double       exp_tol      = 20.0 * sqrt(safe_macheps);
 62
 63    // Create aoclsparse matrix and its descriptor
 64    aoclsparse_matrix     A;
 65    aoclsparse_index_base base = aoclsparse_index_base_zero;
 66    aoclsparse_mat_descr  descr_a;
 67    aoclsparse_operation  trans;
 68    status = aoclsparse_create_dcsr(&A, base, m, n, nnz, icrow, icol, aval);
 69    if(status != aoclsparse_status_success)
 70    {
 71        std::cerr << "Error returned from aoclsparse_create_dcsr, status = " << status << "."
 72                  << std::endl;
 73        return 3;
 74    }
 75    aoclsparse_create_mat_descr(&descr_a);
 76
 77    double alpha = 1.0;
 78    double x[m]  = {1, 1, 1, 1};
 79    double y[m]  = {0, 0, 0, 0};
 80
 81    double              b[m] = {115, 230, 345, 450};
 82    std::vector<double> xref;
 83    xref.assign({1.000006, 1.999687, 2.999374, 3.999038});
 84    std::vector<double> yref;
 85    yref.assign({115.0000000, 229.9639753, 344.9279505, 449.8919266});
 86
 87    // Indicate to access the lower part of matrix A
 88    aoclsparse_set_mat_type(descr_a, aoclsparse_matrix_type_symmetric);
 89    aoclsparse_set_mat_fill_mode(descr_a, aoclsparse_fill_mode_lower);
 90    trans = aoclsparse_operation_none;
 91
 92    status = aoclsparse_set_symgs_hint(A, trans, descr_a, 1);
 93    if(status != aoclsparse_status_success)
 94    {
 95        std::cerr << "Error returned from aoclsparse_set_symgs_hint, status = " << status << "."
 96                  << std::endl;
 97        return 3;
 98    }
 99    status = aoclsparse_optimize(A);
100    if(status != aoclsparse_status_success)
101    {
102        std::cerr << "Error returned from aoclsparse_optimize, status = " << status << "."
103                  << std::endl;
104        return 3;
105    }
106    // Solve
107    status = aoclsparse_dsymgs_mv(trans, A, descr_a, alpha, b, x, y);
108    if(status != aoclsparse_status_success)
109    {
110        std::cerr << "Error returned from aoclsparse_dsymgs_mv, status = " << status << "."
111                  << std::endl;
112        return 3;
113    }
114    // Print and check the result
115    std::cout << "Solving A x = alpha b" << std::endl;
116    std::cout << "  where x and b are dense vectors" << std::endl;
117    std::cout << "  Solution vector x = " << std::endl;
118    std::cout << std::fixed;
119    std::cout.precision(1);
120    ok = true;
121    for(int i = 0; i < m; ++i)
122    {
123        oki = std::abs(x[i] - xref[i]) <= exp_tol;
124        std::cout << x[i] << (oki ? "  " : "! ");
125        ok &= oki;
126    }
127    std::cout << std::endl;
128    std::cout << "  and Product vector y = " << std::endl;
129    for(int i = 0; i < m; ++i)
130    {
131        oki = std::abs(y[i] - yref[i]) <= exp_tol;
132        std::cout << y[i] << (oki ? "  " : "! ");
133        ok &= oki;
134    }
135    std::cout << std::endl;
136
137    // Destroy the aoclsparse memory
138    aoclsparse_destroy_mat_descr(descr_a);
139    aoclsparse_destroy(&A);
140
141    return ok ? 0 : 5;
142}


  1/* ************************************************************************
  2 * Copyright (c) 2024 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 <cmath>
 27#include <complex>
 28#include <iomanip>
 29#include <iostream>
 30#include <limits>
 31#include <vector>
 32
 33int main(void)
 34{
 35    std::cout
 36        << "------------------------------------------------" << std::endl
 37        << " Symmetric Gauss Seidel Preconditioner followed by sparse matrix-vector multiplication"
 38        << std::endl
 39        << " sample program for complex data type" << std::endl
 40        << "------------------------------------------------" << std::endl
 41        << std::endl;
 42
 43    /*
 44     * This example illustrates how to use Symmetric Gauss Seidel API to solve
 45     * linear system of equations. The complex matrix used is
 46     *     | 111+11.1i     2+0.2i          0          0 |
 47     * A = |    2+0.2i  111+11.1i     2+0.2i          0 |
 48     *     |         0     2+0.2i  111+11.1i     2+0.2i |
 49     *     |         0          0     2+0.2i  111+11.1i |
 50     */
 51
 52    // Create a tri-diagonal matrix in CSR format
 53    const aoclsparse_int                   n = 4, m = 4, nnz = 10;
 54    aoclsparse_int                         icrow[5] = {0, 2, 5, 8, 10};
 55    aoclsparse_int                         icol[10] = {0, 1, 0, 1, 2, 1, 2, 3, 2, 3};
 56    std::vector<aoclsparse_double_complex> aval;
 57    // clang-format off
 58    aval.assign({{111,11.1}, {2,0.2}, {2,0.2}, {111,11.1},
 59                 {2,0.2}, {2,0.2}, {111,11.1}, {2,0.2},
 60                 {2,0.2}, {111,11.1}});
 61    // clang-format on
 62    aoclsparse_status status;
 63    bool              oki, ok;
 64
 65    const double macheps      = std::numeric_limits<double>::epsilon();
 66    const double safe_macheps = (double)2.0 * macheps;
 67    double       exp_tol      = 20.0 * sqrt(safe_macheps);
 68
 69    // Create aoclsparse matrix and its descriptor
 70    aoclsparse_matrix     A;
 71    aoclsparse_index_base base = aoclsparse_index_base_zero;
 72    aoclsparse_mat_descr  descr_a;
 73    aoclsparse_operation  trans;
 74    status = aoclsparse_create_zcsr(&A, base, m, n, nnz, icrow, icol, &aval[0]);
 75    if(status != aoclsparse_status_success)
 76    {
 77        std::cerr << "Error returned from aoclsparse_create_dcsr, status = " << status << "."
 78                  << std::endl;
 79        return 3;
 80    }
 81    aoclsparse_create_mat_descr(&descr_a);
 82
 83    aoclsparse_double_complex alpha = {1.0, 0.};
 84
 85    std::vector<aoclsparse_double_complex> x, b, y, xref, yref;
 86    x.assign({{1, 0.1}, {1, 0.1}, {1, 0.1}, {1, 0.1}});
 87    b.assign({{113.850000, 23.000000},
 88              {227.700000, 46.000000},
 89              {341.550000, 69.000000},
 90              {445.500000, 90.000000}});
 91    y.assign({{0, 0.0}, {0, 0.0}, {0, 0.0}, {0, 0.0}});
 92    xref.assign({{1.0000056, 0.1000006},
 93                 {1.9996866, 0.1999687},
 94                 {2.9993739, 0.2999374},
 95                 {3.9990376, 0.3999038}});
 96    yref.assign({{113.8500000, 23.0000000},
 97                 {227.6643355, 45.9927951},
 98                 {341.4786710, 68.9855901},
 99                 {445.3930073, 89.9783853}});
100
101    // Indicate to access the lower part of matrix A
102    aoclsparse_set_mat_type(descr_a, aoclsparse_matrix_type_symmetric);
103    aoclsparse_set_mat_fill_mode(descr_a, aoclsparse_fill_mode_lower);
104    trans = aoclsparse_operation_none;
105
106    status = aoclsparse_set_symgs_hint(A, trans, descr_a, 1);
107    if(status != aoclsparse_status_success)
108    {
109        std::cerr << "Error returned from aoclsparse_set_symgs_hint, status = " << status << "."
110                  << std::endl;
111        return 3;
112    }
113    status = aoclsparse_optimize(A);
114    if(status != aoclsparse_status_success)
115    {
116        std::cerr << "Error returned from aoclsparse_optimize, status = " << status << "."
117                  << std::endl;
118        return 3;
119    }
120    // Solve
121    status = aoclsparse_zsymgs_mv(trans, A, descr_a, alpha, &b[0], &x[0], &y[0]);
122    if(status != aoclsparse_status_success)
123    {
124        std::cerr << "Error returned from aoclsparse_zsymgs, status = " << status << "."
125                  << std::endl;
126        return 3;
127    }
128    // Print and check the result
129    std::cout << "Solving A x = alpha b" << std::endl;
130    std::cout << "  where x and b are dense vectors" << std::endl;
131    std::cout << "  Solution vector x = " << std::endl;
132    std::cout << std::fixed;
133    std::cout.precision(1);
134    ok = true;
135    for(int i = 0; i < m; ++i)
136    {
137        oki = std::abs(x[i].real - xref[i].real) <= exp_tol
138              && std::abs(x[i].imag - xref[i].imag) <= exp_tol;
139        std::cout << "(" << x[i].real << ", " << x[i].imag << "i) " << (oki ? "  " : "! ");
140        ok &= oki;
141    }
142    std::cout << std::endl;
143    std::cout << "  Product vector x = " << std::endl;
144    for(int i = 0; i < m; ++i)
145    {
146        oki = std::abs(y[i].real - yref[i].real) <= exp_tol
147              && std::abs(y[i].imag - yref[i].imag) <= exp_tol;
148        std::cout << "(" << y[i].real << ", " << y[i].imag << "i) " << (oki ? "  " : "! ");
149        ok &= oki;
150    }
151    std::cout << std::endl;
152    // Destroy the aoclsparse memory
153    aoclsparse_destroy_mat_descr(descr_a);
154    aoclsparse_destroy(&A);
155
156    return ok ? 0 : 5;
157}
Parameters:

  • trans – [in] matrix operation to perform on \(A\). Possible values are aoclsparse_operation_none, aoclsparse_operation_transpose, and aoclsparse_operation_conjugate_transpose.

  • A – [in] sparse matrix \(A\) of size \(m\).

  • descr – [in] descriptor of the sparse matrix \(A\).

  • alpha – [in] scalar \(\alpha\).

  • b – [in] dense vector, of size \(m\).

  • x – [out] solution vector \(x,\) dense vector of size \(m\).

  • y – [out] sparse-product vector \(y,\) dense vector of size \(m\).

Return values:

  • aoclsparse_status_success – indicates that the operation completed successfully.

  • aoclsparse_status_invalid_size – informs that either m, n or nnz is invalid. The error code also informs if the given sparse matrix \(A\) is not square.

  • aoclsparse_status_invalid_value – informs that either base, trans, matrix type descr->type or fill mode descr->fill_mode is invalid. If the sparse matrix \(A\) is not of full rank, the error code is returned to indicate that the linear system cannot be solved.

  • aoclsparse_status_invalid_pointer – informs that either descr, A, b, x or y pointer is invalid.

  • aoclsparse_status_not_implemented – this error occurs when the provided matrix’s aoclsparse_fill_mode is aoclsparse_diag_type_unit or the input format is not aoclsparse_csr_mat, or when aoclsparse_matrix_type is aoclsparse_matrix_type_general and trans is aoclsparse_operation_conjugate_transpose.