aoclsparse::trsv() - 5.2 English - 68552

  1/* ************************************************************************
  2 * Copyright (c) 2023-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    Description: This is an example for invoking the C++ interfaces of TRSV
 25                 for complex types.
 26*/
 27#include "aoclsparse.hpp"
 28
 29#include <assert.h>
 30#include <complex>
 31#include <iomanip>
 32#include <iostream>
 33#include <limits>
 34#include <vector>
 35
 36template <typename T>
 37auto get_tolerance()
 38{
 39    if constexpr(std::is_same_v<T, std::complex<float>>)
 40        return 20 * std::numeric_limits<float>::epsilon();
 41    else
 42        return 20 * std::numeric_limits<double>::epsilon();
 43}
 44
 45template <typename T>
 46aoclsparse_int trsv_driver()
 47{
 48    /* Solve two complex linear systems of equations
 49     * Lx = alpha b, and U^Tx = alpha b,
 50     * where L is the lower triangular part of A and
 51     * U is the upper triangular part of A. U^H is
 52     * the conjugate transpose of U (a lower triangular matrix).
 53     * The complex matrix A is the tri-diagonal matrix in CSR format
 54     *
 55     * | 1+3i  2+5i     0     0 |
 56     * | 3     1+2i  2-2i     0 |
 57     * |    0     3     1  2+3i |
 58     * |    0     0  3+1i  1+1i |
 59     */
 60    const aoclsparse_int        n = 4, m = 4, nnz = 10;
 61    std::vector<aoclsparse_int> icrow = {0, 2, 5, 8, 10};
 62    std::vector<aoclsparse_int> icol  = {0, 1, 0, 1, 2, 1, 2, 3, 2, 3};
 63    std::vector<T>              aval  = {{1., 3.},
 64                                         {2., 5.},
 65                                         {3., 0.},
 66                                         {1., 2.},
 67                                         {2., -2.},
 68                                         {3., 0.},
 69                                         {1., 0.},
 70                                         {2., 3.},
 71                                         {3., 1.},
 72                                         {1., 1.}};
 73    aoclsparse_status           status;
 74    std::vector<T>              ref;
 75
 76    // create aoclsparse matrix and its descriptor
 77    aoclsparse_matrix     A;
 78    aoclsparse_index_base base = aoclsparse_index_base_zero;
 79    aoclsparse_mat_descr  descr_a;
 80    aoclsparse_operation  trans = aoclsparse_operation_none;
 81
 82    status = aoclsparse::create_csr(&A, base, m, n, nnz, icrow.data(), icol.data(), aval.data());
 83
 84    if(status != aoclsparse_status_success)
 85    {
 86        std::cerr << "Error creating the matrix, status = " << status << std::endl;
 87        return 1;
 88    }
 89    aoclsparse_create_mat_descr(&descr_a);
 90
 91    /* Solve the lower triangular system Lx = b,
 92     * here alpha=1 and b = [1+i, 4+2i, 4+i, 4].
 93     */
 94    std::vector<T> b = {{1., 1.}, {4., 2.}, {4., 1}, {4., 0}};
 95    aoclsparse_set_mat_type(descr_a, aoclsparse_matrix_type_triangular);
 96    aoclsparse_set_mat_fill_mode(descr_a, aoclsparse_fill_mode_lower);
 97    status = aoclsparse_set_sv_hint(A, trans, descr_a, 1);
 98    if(status != aoclsparse_status_success)
 99    {
100        std::cerr << "Error setting SV hint, status = " << status << std::endl;
101        return 1;
102    }
103    status = aoclsparse_optimize(A);
104    if(status != aoclsparse_status_success)
105    {
106        std::cerr << "Error returned from aoclsparse_optimize, status = " << status << "."
107                  << std::endl;
108        return 2;
109    }
110
111    // Call the triangle solver
112    T alpha = {1.0, 0.};
113    T x[n];
114    status = aoclsparse::trsv(trans, alpha, A, descr_a, b.data(), 1, x, 1);
115    if(status != aoclsparse_status_success)
116    {
117        std::cerr << "Error returned from trsv, status = " << status << "." << std::endl;
118        return 3;
119    }
120
121    // Print and check the result
122    auto tol = get_tolerance<T>();
123    ref.assign({{0.4, -0.2}, {1.6, -0.6}, {-0.8, 2.8}, {0.8, -8.4}});
124    std::cout << "Solving Lx = alpha b: " << std::endl << "  x = ";
125    std::cout << std::fixed;
126    std::cout.precision(1);
127    bool oki, ok = true;
128    for(aoclsparse_int i = 0; i < n; i++)
129    {
130        oki = std::abs(x[i].real() - ref[i].real()) <= tol
131              && std::abs(x[i].imag() - ref[i].imag()) <= tol;
132        std::cout << "(" << x[i].real() << ", " << x[i].imag() << "i) " << (oki ? " " : "!  ");
133        ok &= oki;
134    }
135    std::cout << std::endl;
136
137    /* Solve the lower triangular system U^Hx = b,
138     * here alpha=1-i and b is unchanged.
139     */
140    alpha = {1, -1};
141    // Indicate to use only the conjugate transpose of the upper part of A
142    trans = aoclsparse_operation_conjugate_transpose;
143    aoclsparse_set_mat_fill_mode(descr_a, aoclsparse_fill_mode_upper);
144    status = aoclsparse::trsv(trans, alpha, A, descr_a, b.data(), 1, x, 1);
145    if(status != aoclsparse_status_success)
146    {
147        std::cerr << "Error returned from trsv, status = " << status << "." << std::endl;
148        return 3;
149    }
150
151    // Print and check the result
152    ref.assign({{0.2, 0.6}, {1.4, 0.6}, {3.4, -7.0}, {-1.0, 19.2}});
153    std::cout << std::endl;
154    std::cout << "Solving U^Hx = alpha*b: " << std::endl << "  x = ";
155    for(aoclsparse_int i = 0; i < n; i++)
156    {
157        oki = std::abs(x[i].real() - ref[i].real()) <= tol
158              && std::abs(x[i].imag() - ref[i].imag()) <= tol;
159        std::cout << "(" << x[i].real() << ", " << x[i].imag() << "i) " << (oki ? " " : "!  ");
160        ok &= oki;
161    }
162    std::cout << std::endl;
163
164    // Destroy the aoclsparse memory
165    aoclsparse_destroy_mat_descr(descr_a);
166    aoclsparse_destroy(&A);
167
168    if(ok)
169        return 0;
170    else
171        return 1;
172}
173
174int main()
175{
176    std::cout << "-------------------------------" << std::endl
177              << " Triangle solve sample program" << std::endl
178              << "-------------------------------" << std::endl
179              << std::endl;
180
181    aoclsparse_int ctrsv_err = trsv_driver<std::complex<float>>();
182    aoclsparse_int ztrsv_err = trsv_driver<std::complex<double>>();
183
184    if(ctrsv_err != 0)
185    {
186        std::cerr << "CTRSV driver failed with error code: " << ctrsv_err << std::endl;
187        exit(ctrsv_err);
188    }
189
190    if(ztrsv_err != 0)
191    {
192        std::cerr << "ZTRSV driver failed with error code: " << ztrsv_err << std::endl;
193        exit(ztrsv_err);
194    }
195
196    return 0;
197}