Sparse BLAS level 1, 2, and 3 functions#

aoclsparse_functions.h provides AMD CPU hardware optimized level 1, 2, and 3 Sparse Linear Algebra Subprograms (Sparse BLAS).

Level 1#

The sparse level 1 routines describe operations between a vector in sparse format and a vector in dense format.

This section describes all provided level 1 sparse linear algebra functions.

aoclsparse_?axpyi()#

aoclsparse_status aoclsparse_saxpyi(const aoclsparse_int nnz, const float a, const float *x, const aoclsparse_int *indx, float *y)#

aoclsparse_status aoclsparse_daxpyi(const aoclsparse_int nnz, const double a, const double *x, const aoclsparse_int *indx, double *y)#

aoclsparse_status aoclsparse_caxpyi(const aoclsparse_int nnz, const void *a, const void *x, const aoclsparse_int *indx, void *y)#

aoclsparse_status aoclsparse_zaxpyi(const aoclsparse_int nnz, const void *a, const void *x, const aoclsparse_int *indx, void *y)#

A variant of sparse vector-vector addition between compressed sparse vector and dense vector.

aoclsparse_?axpyi adds a scalar multiple of compressed sparse vector to a dense vector.

Let \(y\in R^m\) (or \(C^m\)) be a dense vector, \(x\) be a compressed sparse vector and \(I_x\) be the nonzero indices set for x of length at least nnz described by indx, then

\[ y_{I_{x_{i}}} = a\,x_i + y_{I_{x_{i}}}, \quad i\in\{1,\ldots,{\bf\mathsf{nnz}}\}. \]

Example (tests/examples/sample_axpyi.cpp)

/* ************************************************************************
 * Copyright (c) 2023-2024 Advanced Micro Devices, Inc.
 *
 * Permission is hereby granted, free of charge, to any person obtaining a copy
 * of this software and associated documentation files (the "Software"), to deal
 * in the Software without restriction, including without limitation the rights
 * to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
 * copies of the Software, and to permit persons to whom the Software is
 * furnished to do so, subject to the following conditions:
 *
 * The above copyright notice and this permission notice shall be included in
 * all copies or substantial portions of the Software.
 *
 * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
 * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
 * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
 * AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
 * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
 * OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
 * THE SOFTWARE.
 *
 * ************************************************************************ */

#include "aoclsparse.h"

#include <cfloat>
#include <cmath>
#include <complex>
#include <iomanip>
#include <iostream>
#include <vector>

/* Sample program to illustrate the usage of axpyi API which performs vector-vector
 * addition operation of the form y = ax + y, where
 * a: scalar value
 * x: compressed sparse vector
 * y: dense vector
 * Below example is for complex double precision. Using other precisions are also similar.
 */
int main(void)
{
    std::cout << "--------------------------------" << std::endl
              << "----- axpyi sample program -----" << std::endl
              << "--------------------------------" << std::endl
              << std::endl;

    // Number of non-zeros of the sparse vector
    const aoclsparse_int nnz = 3;
    // Scalar value
    const std::complex<double> a = {10, 20};

    // Sparse index vector (does not need to be ordered)
    std::vector<aoclsparse_int> indx = {0, 6, 3};
    // Sparse value vector in compressed form
    std::vector<std::complex<double>> x{{2, 3}, {4, 1}, {5, 1}};

    // Output vector
    std::vector<std::complex<double>> y{{0, 0}, {0, 0}, {0, 0}, {0, 0}, {0, 0}, {0, 0}, {65, -999}};

    // Expected output vector
    std::vector<std::complex<double>> y_exp{
        {-40, 70}, {0, 0}, {0, 0}, {30, 110}, {0, 0}, {0, 0}, {85, -909}};

    aoclsparse_int ny = y.size();

    aoclsparse_status status;

    std::cout << "Invoking aoclsparse_zaxpyi...\n";
    //Invoke complex axpyi
    status = aoclsparse_zaxpyi(nnz, &a, x.data(), indx.data(), y.data());
    if(status != aoclsparse_status_success)
    {
        std::cerr << "Error returned from aoclsparse_zaxpyi, status = " << status << "."
                  << std::endl;
        return 3;
    }

    // Check and print the result
    std::cout << "vector-vector addition of y: " << std::endl
              << std::setw(10) << "y"
              << "   " << std::setw(20) << "expected y" << std::endl;
    std::cout << std::fixed;
    std::cout.precision(4);
    bool oki, ok = true;
    for(aoclsparse_int i = 0; i < ny; i++)
    {
        oki = y[i] == y_exp[i];
        ok &= oki;
        std::cout << std::setw(15) << y[i] << std::setw(5) << "" << std::setw(15) << y_exp[i]
                  << std::setw(5) << (oki ? "" : " ! Error") << std::endl;
    }
    return (ok ? 0 : 6);
}

Note

The contents of the vectors are not checked for NaNs.

Parameters

nnz – [in] The number of elements in \(x\) and indx.
a – [in] Scalar value.
x – [in] Sparse vector stored in compressed form of at least nnz elements.
indx – [in] Nonzero indices set, \(I_x\), of x described by this array of length at least nnz. The elements in this vector are only checked for non-negativity. The caller should make sure that all indices are less than the size of y. Array is assumed to be in zero base.
y – [inout] Array of at least \(\max(I_{x_i}, i \in \{ 1,\ldots,{\bf\mathsf{nnz}}\})\) elements.

Return values

aoclsparse_status_success – The operation completed successfully.
aoclsparse_status_invalid_pointer – At least one of the pointers x, indx, y is invalid.
aoclsparse_status_invalid_size – Indicates that provided nnz is less than zero.
aoclsparse_status_invalid_index_value – At least one of the indices in indx is negative.

aoclsparse_?dotci()#

aoclsparse_status aoclsparse_cdotci(const aoclsparse_int nnz, const void *x, const aoclsparse_int *indx, const void *y, void *dot)#

aoclsparse_status aoclsparse_zdotci(const aoclsparse_int nnz, const void *x, const aoclsparse_int *indx, const void *y, void *dot)#

Sparse conjugate dot product for single and double data precision complex types.

aoclsparse_cdotci() (complex float) and aoclsparse_zdotci() (complex double) compute the dot product of the conjugate of a complex vector stored in a compressed format and a complex dense vector. Let \(x\) and \(y\) be respectively a sparse and dense vectors in \(C^m\) with indx ( \(I_x\)) the nonzero indices array of x of length at least nnz that is used to index into the entries of dense vector \(y\), then these functions return

\[ {\bf\mathsf{dot}} = \sum_{i=0}^{{\bf\mathsf{nnz}}-1} \overline{\,x_i\,} \cdot y_{I_{x_i}}. \]

Example (tests/examples/sample_dotp.cpp)

/* ************************************************************************
 * Copyright (c) 2023 Advanced Micro Devices, Inc.
 *
 * Permission is hereby granted, free of charge, to any person obtaining a copy
 * of this software and associated documentation files (the "Software"), to deal
 * in the Software without restriction, including without limitation the rights
 * to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
 * copies of the Software, and to permit persons to whom the Software is
 * furnished to do so, subject to the following conditions:
 *
 * The above copyright notice and this permission notice shall be included in
 * all copies or substantial portions of the Software.
 *
 * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
 * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
 * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
 * AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
 * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
 * OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
 * THE SOFTWARE.
 *
 * ************************************************************************ */

#include <complex>
#include <iostream>
#include <vector>

/* One can use std::complex<double> instead of aoclsparse_double_complex
#define STD_CMPLX  (or pass it as a compilation flag -DSTD_CMPLX)
#ifdef STD_CMPLX
#define aoclsparse_double_complex std::complex<double>
#endif
*/

#include "aoclsparse.h"

int main(void)
{
    aoclsparse_status                      status;
    const aoclsparse_int                   nnz  = 3;
    std::vector<aoclsparse_int>            indx = {0, 1, 2};
    std::vector<aoclsparse_double_complex> x{{1, 1}, {2, 2}, {3, 3}};
    std::vector<aoclsparse_double_complex> y{
        {1, 1}, {1, 2}, {2, 3}, {1, 0}, {4, 1}, {1.2, -1}, {7, -1}, {0, 2}, {-2, 3}};
    aoclsparse_double_complex dot_exp{-5, 23};
    aoclsparse_double_complex dotc_exp{23, 5};
    aoclsparse_double_complex dotp;
    bool                      ok;

    std::cout << "Invoking aoclsparse_zdotui...\n";

    //Invoke complex double dot product
    status = aoclsparse_zdotui(nnz, x.data(), indx.data(), y.data(), &dotp);
    if(status != aoclsparse_status_success)
    {
        std::cerr << "Error returned from aoclsparse_cdotui, status = " << status << "."
                  << std::endl;
        return 1;
    }

#ifdef STD_CMPLX
    std::cout << "Dot product: " << dotp << ", Expected dot product: " << dot_exp << std::endl;
    ok = (dotp == dot_exp);
#else
    std::cout << "Dot product: (" << dotp.real << ", " << dotp.imag << "), Expected dot product: ("
              << dot_exp.real << ", " << dot_exp.imag << ")" << std::endl;
    ok = ((dotp.real == dot_exp.real) && (dotp.imag == dot_exp.imag));
#endif

    std::cout << "Invoking aoclsparse_zdotci...\n";

    //Invoke complex double dot product
    status = aoclsparse_zdotci(nnz, x.data(), indx.data(), y.data(), &dotp);
    if(status != aoclsparse_status_success)
    {
        std::cerr << "Error returned from aoclsparse_cdotci, status = " << status << "."
                  << std::endl;
        return 2;
    }
#ifdef STD_CMPLX
    std::cout << "Conjugated dot product: " << dotp
              << ", Expected conjugated dot product: " << dotc_exp << std::endl;
    ok &= (dotp == dotc_exp);
    return ok ? 0 : 3;

#else
    std::cout << "Dot product: (" << dotp.real << ", " << dotp.imag << "), Expected dot product: ("
              << dotc_exp.real << ", " << dotc_exp.imag << ")" << std::endl;
    ok &= ((dotp.real == dotc_exp.real) && (dotp.imag == dotc_exp.imag));
    return ok ? 0 : 4;
#endif
}

Note

The contents of the vectors are not checked for NaNs.

Parameters

nnz – [in] The number of elements (length) of vectors x and indx.
x – [in] Array of at least nnz complex elements.
indx – [in] Nonzero indices set, \(I_x\), of x described by this array of length at least nnz. Each entry must contain a valid index into y and be unique. The entries of indx are not checked for validity.
y – [in] Array of at least \(\max(I_{x_i}, i \in \{ 1,\ldots,{\bf\mathsf{nnz}}\})\) complex elements.
dot – [out] The dot product of the conjugate of \(x\) and \(y\) when nnz \(> 0\). If nnz \(\le 0\), dot is set to 0.

Return values

aoclsparse_status_success – The operation completed successfully.
aoclsparse_status_invalid_pointer – At least one of the pointers x, indx, y, dot is invalid.
aoclsparse_status_invalid_size – Indicates that the provided nnz is not positive.

aoclsparse_?dotui()#

aoclsparse_status aoclsparse_cdotui(const aoclsparse_int nnz, const void *x, const aoclsparse_int *indx, const void *y, void *dot)#

aoclsparse_status aoclsparse_zdotui(const aoclsparse_int nnz, const void *x, const aoclsparse_int *indx, const void *y, void *dot)#

Sparse dot product for single and double data precision complex types.

aoclsparse_cdotui() (complex float) and aoclsparse_zdotui() (complex double) compute the dot product of a complex vector stored in a compressed format and a complex dense vector. Let \(x\) and \(y\) be respectively a sparse and dense vectors in \(C^m\) with indx ( \(I_x\)) the nonzero indices array of x of length at least nnz that is used to index into the entries of dense vector \(y\), then these functions return

\[ {\bf\mathsf{dot}} = \sum_{i=0}^{{\bf\mathsf{nnz}}-1} x_{i} \cdot y_{I_{x_i}}. \]

Example (tests/examples/sample_dotp.cpp)

/* ************************************************************************
 * Copyright (c) 2023 Advanced Micro Devices, Inc.
 *
 * Permission is hereby granted, free of charge, to any person obtaining a copy
 * of this software and associated documentation files (the "Software"), to deal
 * in the Software without restriction, including without limitation the rights
 * to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
 * copies of the Software, and to permit persons to whom the Software is
 * furnished to do so, subject to the following conditions:
 *
 * The above copyright notice and this permission notice shall be included in
 * all copies or substantial portions of the Software.
 *
 * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
 * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
 * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
 * AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
 * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
 * OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
 * THE SOFTWARE.
 *
 * ************************************************************************ */

#include <complex>
#include <iostream>
#include <vector>

/* One can use std::complex<double> instead of aoclsparse_double_complex
#define STD_CMPLX  (or pass it as a compilation flag -DSTD_CMPLX)
#ifdef STD_CMPLX
#define aoclsparse_double_complex std::complex<double>
#endif
*/

#include "aoclsparse.h"

int main(void)
{
    aoclsparse_status                      status;
    const aoclsparse_int                   nnz  = 3;
    std::vector<aoclsparse_int>            indx = {0, 1, 2};
    std::vector<aoclsparse_double_complex> x{{1, 1}, {2, 2}, {3, 3}};
    std::vector<aoclsparse_double_complex> y{
        {1, 1}, {1, 2}, {2, 3}, {1, 0}, {4, 1}, {1.2, -1}, {7, -1}, {0, 2}, {-2, 3}};
    aoclsparse_double_complex dot_exp{-5, 23};
    aoclsparse_double_complex dotc_exp{23, 5};
    aoclsparse_double_complex dotp;
    bool                      ok;

    std::cout << "Invoking aoclsparse_zdotui...\n";

    //Invoke complex double dot product
    status = aoclsparse_zdotui(nnz, x.data(), indx.data(), y.data(), &dotp);
    if(status != aoclsparse_status_success)
    {
        std::cerr << "Error returned from aoclsparse_cdotui, status = " << status << "."
                  << std::endl;
        return 1;
    }

#ifdef STD_CMPLX
    std::cout << "Dot product: " << dotp << ", Expected dot product: " << dot_exp << std::endl;
    ok = (dotp == dot_exp);
#else
    std::cout << "Dot product: (" << dotp.real << ", " << dotp.imag << "), Expected dot product: ("
              << dot_exp.real << ", " << dot_exp.imag << ")" << std::endl;
    ok = ((dotp.real == dot_exp.real) && (dotp.imag == dot_exp.imag));
#endif

    std::cout << "Invoking aoclsparse_zdotci...\n";

    //Invoke complex double dot product
    status = aoclsparse_zdotci(nnz, x.data(), indx.data(), y.data(), &dotp);
    if(status != aoclsparse_status_success)
    {
        std::cerr << "Error returned from aoclsparse_cdotci, status = " << status << "."
                  << std::endl;
        return 2;
    }
#ifdef STD_CMPLX
    std::cout << "Conjugated dot product: " << dotp
              << ", Expected conjugated dot product: " << dotc_exp << std::endl;
    ok &= (dotp == dotc_exp);
    return ok ? 0 : 3;

#else
    std::cout << "Dot product: (" << dotp.real << ", " << dotp.imag << "), Expected dot product: ("
              << dotc_exp.real << ", " << dotc_exp.imag << ")" << std::endl;
    ok &= ((dotp.real == dotc_exp.real) && (dotp.imag == dotc_exp.imag));
    return ok ? 0 : 4;
#endif
}

Note

The contents of the vectors are not checked for NaNs.

Parameters

nnz – [in] The number of elements (length) of vectors \(x\) and \(indx\).
x – [in] Array of at least nnz complex elements.
indx – [in] Nonzero indices set, \(I_x\), of x described by this array of length at least nnz. Each entry must contain a valid index into y and be unique. The entries of indx are not checked for validity.
y – [in] Array of at least \(\max(I_{x_i}, i \in \{ 1,\ldots,{\bf\mathsf{nnz}}\})\) complex elements.
dot – [out] The dot product of x and y when nnz \(> 0\). If nnz \(\le 0\), dot is set to 0.

Return values

aoclsparse_status_success – The operation completed successfully.
aoclsparse_status_invalid_pointer – At least one of the pointers x, indx, y, dot is invalid.
aoclsparse_status_invalid_size – Indicates that the provided nnz is not positive.

aoclsparse_?doti()#

float aoclsparse_sdoti(const aoclsparse_int nnz, const float *x, const aoclsparse_int *indx, const float *y)#

double aoclsparse_ddoti(const aoclsparse_int nnz, const double *x, const aoclsparse_int *indx, const double *y)#

Sparse dot product for single and double data precision real types.

aoclsparse_sdoti() and aoclsparse_ddoti() compute the dot product of a real vector stored in a compressed format and a real dense vector. Let \(x\) and \(y\) be respectively a sparse and dense vectors in \(R^m\) with indx ( \(I_x\)) an indices array of length at least nnz that is used to index into the entries of dense vector \(y\), then these functions return

\[ {\bf\mathsf{dot}} = \sum_{i=0}^{{\bf\mathsf{nnz}}-1} x_{i} \cdot y_{I_{x_i}}. \]

Note

The contents of the vectors are not checked for NaNs.

Parameters

nnz – [in] The number of elements to access in vectors x and indx.
x – [in] Array of at least nnz elements.
indx – [in] Nonzero indices set, \(I_x\), of x described by this array of length at least nnz. Each entry must contain a valid index into y and be unique. The entries of indx are not checked for validity.
y – [in] Array of at least \(\max(I_{x_i}, i \in \{ 1,\ldots,{\bf\mathsf{nnz}}\})\) elements.

Return values

dot – Value of the dot product if nnz is positive, otherwise returns 0.

aoclsparse_?sctr()#

aoclsparse_status aoclsparse_ssctr(const aoclsparse_int nnz, const float *x, const aoclsparse_int *indx, float *y)#

aoclsparse_status aoclsparse_dsctr(const aoclsparse_int nnz, const double *x, const aoclsparse_int *indx, double *y)#

aoclsparse_status aoclsparse_csctr(const aoclsparse_int nnz, const void *x, const aoclsparse_int *indx, void *y)#

aoclsparse_status aoclsparse_zsctr(const aoclsparse_int nnz, const void *x, const aoclsparse_int *indx, void *y)#

Sparse scatter for single and double precision real and complex types.

aoclsparse_?sctr scatter the elements of a compressed sparse vector into a dense vector.

Let \(y\in R^m\) (or \(C^m\)) be a dense vector, and \(x\) be a compressed sparse vector with \(I_x\) be its nonzero indices set of length at least nnz and described by the array indx, then

\[ y_{I_{x_{i}}} = x_i, \quad i\in\{1,\ldots,{\bf\mathsf{nnz}}\}. \]

Example (tests/examples/sample_sctr.cpp)

/* ************************************************************************
 * Copyright (c) 2023 Advanced Micro Devices, Inc.
 *
 * Permission is hereby granted, free of charge, to any person obtaining a copy
 * of this software and associated documentation files (the "Software"), to deal
 * in the Software without restriction, including without limitation the rights
 * to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
 * copies of the Software, and to permit persons to whom the Software is
 * furnished to do so, subject to the following conditions:
 *
 * The above copyright notice and this permission notice shall be included in
 * all copies or substantial portions of the Software.
 *
 * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
 * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
 * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
 * AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
 * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
 * OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
 * THE SOFTWARE.
 *
 * ************************************************************************ */

#include "aoclsparse.h"

#include <cfloat>
#include <cmath>
#include <complex>
#include <iomanip>
#include <iostream>
#include <vector>

// Sample program to illustrate the usage of scatter of a compressed sparse vector to the full storage form (complex double precision)
// Using other precisions are also similar to below example
int main(void)
{
    std::cout << "--------------------------------" << std::endl
              << "---- Scatter sample program ----" << std::endl
              << "--------------------------------" << std::endl
              << std::endl;

    // Number of non-zeros of the sparse vector
    const aoclsparse_int nnz = 3;

    // Sparse index vector (does not need to be ordered)
    std::vector<aoclsparse_int> indx = {0, 3, 6};
    // Sparse value vector in compressed form
    std::vector<std::complex<double>> x{{1.01, -1.13}, {2.4, -2.0}, {-0.3, 1.3}};

    // Output vector
    std::vector<std::complex<double>> y{{0, 0}, {0, 0}, {0, 0}, {0, 0}, {0, 0}, {0, 0}, {0, 0}};

    // Expected output vector
    std::vector<std::complex<double>> y_exp{
        {1.01, -1.13}, {0, 0}, {0, 0}, {2.4, -2.0}, {0, 0}, {0, 0}, {-0.3, 1.3}};

    aoclsparse_int ny = y.size();

    aoclsparse_status status;

    std::cout << "Invoking aoclsparse_zsctr...\n";
    //Invoke complex scatter
    status = aoclsparse_zsctr(nnz, x.data(), indx.data(), y.data());
    if(status != aoclsparse_status_success)
    {
        std::cerr << "Error returned from aoclsparse_zsctr, status = " << status << "."
                  << std::endl;
        return 3;
    }

    // Check and print the result
    std::cout << "Scatter to vector y: " << std::endl
              << std::setw(10) << "y"
              << "   " << std::setw(17) << "expected y" << std::endl;
    std::cout << std::fixed;
    std::cout.precision(4);
    bool oki, ok = true;
    for(aoclsparse_int i = 0; i < ny; i++)
    {
        oki = y[i] == y_exp[i];
        ok &= oki;
        std::cout << std::setw(11) << y[i] << std::setw(3) << "" << std::setw(11) << y_exp[i]
                  << std::setw(2) << (oki ? "" : " ! Error") << std::endl;
    }
    return (ok ? 0 : 6);
}

Note

The contents of the vectors are not checked for NaNs.

Parameters

nnz – [in] The number of elements to use from \(x\) and \({\bf\mathsf{indx}}\).
x – [in] Dense array of at least size \({\bf\mathsf{nnz}}\). The first \({\bf\mathsf{nnz}}\) elements are to be scattered.
indx – [in] Nonzero index set for \(x\) of size at least \({\bf\mathsf{nnz}}\). The first \({\bf\mathsf{nnz}}\) indices are used for the scattering. The elements in this vector are only checked for non-negativity. The user should make sure that index is less than the size of y.
y – [out] Array of at least \(\max(I_{x_i}, i \in \{ 1,\ldots,{\bf\mathsf{nnz}}\})\) elements.

Return values

aoclsparse_status_success – The operation completed successfully.
aoclsparse_status_invalid_pointer – At least one of the pointers x, indx, y is invalid.
aoclsparse_status_invalid_size – Indicates that provided nnz is less than zero.
aoclsparse_status_invalid_index_value – At least one of the indices in indx is negative.

sparse_?sctrs()#

aoclsparse_status aoclsparse_ssctrs(const aoclsparse_int nnz, const float *x, aoclsparse_int stride, float *y)#

aoclsparse_status aoclsparse_dsctrs(const aoclsparse_int nnz, const double *x, aoclsparse_int stride, double *y)#

aoclsparse_status aoclsparse_csctrs(const aoclsparse_int nnz, const void *x, aoclsparse_int stride, void *y)#

aoclsparse_status aoclsparse_zsctrs(const aoclsparse_int nnz, const void *x, aoclsparse_int stride, void *y)#

Sparse scatter with stride for real/complex single and double data precisions.

aoclsparse_?sctrs scatters the elements of a compressed sparse vector into a dense vector using a stride.

Let \(y\) be a dense vector of length \(n>0\), \(x\) be a compressed sparse vector with nnz > 0 nonzeros, and stride be a striding distance, then

\[ y_{{\bf\mathsf{stride}} \times i} = x_i,\quad i\in\{1,\ldots,{\bf\mathsf{nnz}}\}.\]

Note

Contents of the vector x are accessed but not checked.

Parameters

nnz – [in] Number of nonzero elements to access in \(x\).
x – [in] Array of at least nnz elements. The first nnz elements are to be scattered into y.
stride – [in] (Positive) striding distance used to store elements in vector y.
y – [out] Array of size at least stride \(\times\) nnz.

Return values

aoclsparse_status_success – The operation completed successfully.
aoclsparse_status_invalid_pointer – At least one of the pointers x, y is invalid.
aoclsparse_status_invalid_size – Indicates that one or more of the values provided in nnz or stride is not positive.

aoclsparse_?roti()#

aoclsparse_status aoclsparse_sroti(const aoclsparse_int nnz, float *x, const aoclsparse_int *indx, float *y, const float c, const float s)#

aoclsparse_status aoclsparse_droti(const aoclsparse_int nnz, double *x, const aoclsparse_int *indx, double *y, const double c, const double s)#

Apply Givens rotation to single or double precision real vectors.

aoclsparse_sroti() and aoclsparse_droti() apply the Givens rotation on elements of two real vectors.

Let \(y\in R^m\) be a vector in full storage form, \(x\) be a vector in a compressed form and \(I_x\) its nonzero indices set of length at least nnz described by the array indx, then

\[ x_i = {\bf\mathsf{c}} * x_i + {\bf\mathsf{s}} * y_{I_{x_{i}}}, \]

\[ y_{I_{x_{i}}} = {\bf\mathsf{c}} * y_{I_{x_{i}}} - {\bf\mathsf{s}} * x_i, \]

for \(i\in 1, \ldots, {\bf\mathsf{nnz}}\). The elements c, s are scalars.

Example (tests/examples/sample_roti.cpp)

/* ************************************************************************
 * Copyright (c) 2023 Advanced Micro Devices, Inc.
 *
 * Permission is hereby granted, free of charge, to any person obtaining a copy
 * of this software and associated documentation files (the "Software"), to deal
 * in the Software without restriction, including without limitation the rights
 * to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
 * copies of the Software, and to permit persons to whom the Software is
 * furnished to do so, subject to the following conditions:
 *
 * The above copyright notice and this permission notice shall be included in
 * all copies or substantial portions of the Software.
 *
 * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
 * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
 * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
 * AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
 * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
 * OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
 * THE SOFTWARE.
 *
 * ************************************************************************ */

#include "aoclsparse.h"

#include <cfloat>
#include <cmath>
#include <iomanip>
#include <iostream>
#include <limits>
#include <vector>

// Sample program to illustrate the usage of Givens rotation on a compressed sparse vector and a full storage vector(double precision)
int main(void)
{
    std::cout << "-----------------------------" << std::endl
              << "---- roti sample program ----" << std::endl
              << "-----------------------------" << std::endl
              << std::endl;

    // Input data

    // Number of non-zeros of the sparse vector
    const aoclsparse_int nnz = 5;

    // Sparse index vector (does not need to be ordered)
    std::vector<aoclsparse_int> indx{0, 3, 6, 7, 9};

    // Sparse value vector in compressed form
    std::vector<double> x{-0.75, 4, -9.5, 46, 1.25};

    // Output vector
    std::vector<double> y{-0.75, 0, 0, 4, 0, 0, -9.5, 46, 0, 1.25};

    // Expected output vectors
    std::vector<double> y_exp{0, 0, 0, 0, 0, 0, 0, 0, 0, 0};
    std::vector<double> x_exp{-3, 16, -38, 184, 5};

    const double c = 2;
    const double s = 2;

    aoclsparse_int ny = y.size();
    aoclsparse_int nx = x.size();

    aoclsparse_status status;
    std::cout << "Invoking aoclsparse_droti...\n";
    //Invoke aoclsparse_droti to apply Givens rotation
    status = aoclsparse_droti(nnz, x.data(), indx.data(), y.data(), c, s);

    if(status != aoclsparse_status_success)
    {
        std::cerr << "Error returned from aoclsparse_droti, status = " << status << "."
                  << std::endl;
        return 3;
    }

    // Check and print the result
    std::cout << "The vector y after Givens rotation: " << std::endl
              << std::setw(11) << "y"
              << "   " << std::setw(11) << "expected" << std::endl;
    std::cout << std::fixed;

    std::cout.precision(2);

    bool oki, ok_y = true, ok_x = true;

    //Initializing precision tolerance range for double
    const double tol = std::sqrt(std::numeric_limits<double>::epsilon());

    for(aoclsparse_int i = 0; i < ny; i++)
    {
        oki = std::abs(y[i] - y_exp[i]) <= tol;
        ok_y &= oki;
        std::cout << std::setw(11) << y[i] << std::setw(3) << "" << std::setw(11) << y_exp[i]
                  << std::setw(2) << (oki ? "" : " !") << std::endl;
    }
    std::cout << std::endl;
    std::cout << "The vector x after Givens rotation: " << std::endl
              << std::setw(11) << "x"
              << "   " << std::setw(11) << "expected" << std::endl;
    std::cout << std::fixed;
    for(aoclsparse_int i = 0; i < nx; i++)
    {
        oki = std::abs(x[i] - x_exp[i]) <= tol;
        ok_x &= oki;
        std::cout << std::setw(11) << x[i] << std::setw(3) << "" << std::setw(11) << x_exp[i]
                  << std::setw(2) << (oki ? "" : " !") << std::endl;
    }
    return ((ok_y && ok_x) ? 0 : 6);
}

Note

The contents of the vectors are not checked for NaNs.

Parameters

nnz – [in] The number of elements to use from \(x\) and \({\bf\mathsf{indx}}\).
x – [inout] Array \(x\) of at least \({\bf\mathsf{nnz}}\) elements in compressed form. The elements of the array are updated after applying the Givens rotation.
indx – [in] Nonzero index set of \(x\), \(I_x\), with at least \({\bf\mathsf{nnz}}\) elements. The first \({\bf\mathsf{nnz}}\) elements are used to apply the Givens rotation. The elements in this vector are only checked for non-negativity. The caller should make sure that each entry is less than the size of y and are all distinct.
y – [inout] Dense array of at least \(\max(I_{x_i}, \text{ for } i \in \{ 1,\ldots, {\bf\mathsf{nnz}}\})\) elements in full storage form. The elements of the array are updated after applying the Givens rotation.
c – [in] A scalar.
s – [in] A scalar.

Return values

aoclsparse_status_success – The operation completed successfully.
aoclsparse_status_invalid_pointer – At least one of the pointers x, indx, y is invalid.
aoclsparse_status_invalid_size – Indicates that provided nnz is less than zero.
aoclsparse_status_invalid_index_value – At least one of the indices in indx is negative. With this error, the values of vectors x and y are undefined.

aoclsparse_?gthr()#

aoclsparse_status aoclsparse_sgthr(aoclsparse_int nnz, const float *y, float *x, const aoclsparse_int *indx)#

aoclsparse_status aoclsparse_dgthr(aoclsparse_int nnz, const double *y, double *x, const aoclsparse_int *indx)#

aoclsparse_status aoclsparse_cgthr(aoclsparse_int nnz, const void *y, void *x, const aoclsparse_int *indx)#

aoclsparse_status aoclsparse_zgthr(aoclsparse_int nnz, const void *y, void *x, const aoclsparse_int *indx)#

Gather elements from a dense vector and store them into a sparse vector.

The aoclsparse_?gthr is a group of functions that gather the elements indexed in indx from the dense vector y into the sparse vector x.

Let \(y\in R^m\) (or \(C^m\)) be a dense vector, \(x\) be a sparse vector from the same space and \(I_x\) be a set of indices of size \(0<\) nnz \(\le m\) described by indx, then

\[ x_i = y_{I_{x_i}}, i\in\{1,\ldots,{\bf\mathsf{nnz}}\}. \]

For double precision complex vectors use aoclsparse_zgthr() and for single precision complex vectors use aoclsparse_cgthr().

Example - Complex space (tests/examples/sample_zgthr.cpp)

/* ************************************************************************
 * Copyright (c) 2023 Advanced Micro Devices, Inc.
 *
 * Permission is hereby granted, free of charge, to any person obtaining a copy
 * of this software and associated documentation files (the "Software"), to deal
 * in the Software without restriction, including without limitation the rights
 * to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
 * copies of the Software, and to permit persons to whom the Software is
 * furnished to do so, subject to the following conditions:
 *
 * The above copyright notice and this permission notice shall be included in
 * all copies or substantial portions of the Software.
 *
 * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
 * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
 * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
 * AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
 * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
 * OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
 * THE SOFTWARE.
 *
 * ************************************************************************ */

#include "aoclsparse.h"

#include <cfloat>
#include <cmath>
#include <complex>
#include <iomanip>
#include <iostream>
#include <vector>

int main()
{
    std::cout << "-------------------------------" << std::endl
              << "---- Gather sample program ----" << std::endl
              << "-------------------------------" << std::endl
              << std::endl;

    // Input data

    // Number of non-zeros of the sparse vector
    const aoclsparse_int nnz = 3;

    // Sparse index vector (does not need to be ordered)
    std::vector<aoclsparse_int> indx = {0, 5, 3};

    // Sparse value vector
    std::vector<std::complex<double>> x(3);

    // Dense vector
    std::vector<std::complex<double>> y{
        {1, 1}, {1, 2}, {2, 3}, {1, 0}, {4, 1}, {1.2, -1}, {7, -1}, {0, 2}, {-2, 3}};

    // Expected sparse value vector
    std::vector<std::complex<double>> x_exp{{1, 1}, {1.2, -1}, {1, 0}};

    // Expected dense vector
    std::vector<std::complex<double>> y_exp{
        {0, 0}, {1, 2}, {2, 3}, {0, 0}, {4, 1}, {0, 0}, {7, -1}, {0, 2}, {-2, 3}};

    aoclsparse_int ny = y.size();

    aoclsparse_status status;

    // Call the gather function
    status = aoclsparse_zgthr(nnz, y.data(), x.data(), indx.data());
    if(status != aoclsparse_status_success)
    {
        std::cerr << "Error returned from aoclsparse_zgthr, status = " << status << "."
                  << std::endl;
        return 1;
    }

    // Check and print the result
    std::cout << "Gather from vector y: " << std::endl
              << std::setw(11) << "x"
              << "   " << std::setw(11) << "expected" << std::endl;
    std::cout << std::fixed;
    std::cout.precision(1);
    bool oki, ok = true;
    for(aoclsparse_int i = 0; i < nnz; i++)
    {
        oki = x[i] == x_exp[i];
        ok &= oki;
        std::cout << std::setw(11) << x[i] << std::setw(3) << "" << std::setw(11) << x_exp[i]
                  << std::setw(2) << (oki ? "" : " !") << std::endl;
    }

    // Call the gather with zero function
    status = aoclsparse_zgthrz(nnz, y.data(), x.data(), indx.data());
    if(status != aoclsparse_status_success)
    {
        std::cerr << "Error returned from aoclsparse_zgthr, status = " << status << "."
                  << std::endl;
        return 3;
    }

    // Check and print the result
    std::cout << std::endl
              << "Gather from vector y: " << std::endl
              << std::setw(11) << "y" << std::setw(3) << "" << std::setw(11) << "expected"
              << std::setw(3) << "" << std::setw(11) << "x" << std::setw(3) << "" << std::setw(11)
              << "expected" << std::endl;
    std::cout << std::fixed;
    std::cout.precision(1);
    bool oky;
    for(aoclsparse_int i = 0; i < ny; i++)
    {
        oky = y[i] == y_exp[i];
        ok &= oky;
        std::cout << std::setw(11) << y[i] << std::setw(3) << "" << std::setw(11) << y_exp[i]
                  << std::setw(2) << (oky ? "" : " !");
        if(i < nnz)
        {
            oki = x[i] == x_exp[i];
            ok &= oki;
            std::cout << " " << std::setw(11) << x[i] << std::setw(3) << "" << std::setw(11)
                      << x_exp[i] << std::setw(2) << (oki ? "" : " !");
        }
        std::cout << std::endl;
    }

    return (ok ? 0 : 6);
}

Note

These functions assume that the indices stored in indx are less than \(m\) without duplicate elements, and that x and indx are pointers to vectors of size at least nnz.

Parameters

nnz – [in] number of non-zero entries of \(x\). If nnz is zero, then none of the entries of vectors x, y, and indx are touched.
y – [in] pointer to dense vector \(y\) of size at least \(m\).
x – [out] pointer to sparse vector \(x\) with at least nnz non-zero elements.
indx – [in] index vector of size nnz, containing the indices of the non-zero values of \(x\). Indices should range from 0 to \(m-1\), need not be ordered. The elements in this vector are only checked for non-negativity. The user should make sure that no index is out-of-bound and that it does not contains any duplicates.

Return values

aoclsparse_status_success – the operation completed successfully
aoclsparse_status_invalid_size – nnz parameter value is negative
aoclsparse_status_invalid_pointer – at least one of the pointers y, x or indx is invalid
aoclsparse_status_invalid_index_value – at least one of the indices in indx is negative

aoclsparse_?gthrz()#

aoclsparse_status aoclsparse_sgthrz(aoclsparse_int nnz, float *y, float *x, const aoclsparse_int *indx)#

aoclsparse_status aoclsparse_dgthrz(aoclsparse_int nnz, double *y, double *x, const aoclsparse_int *indx)#

aoclsparse_status aoclsparse_cgthrz(aoclsparse_int nnz, void *y, void *x, const aoclsparse_int *indx)#

aoclsparse_status aoclsparse_zgthrz(aoclsparse_int nnz, void *y, void *x, const aoclsparse_int *indx)#

Gather and zero out elements from a dense vector and store them into a sparse vector.

The aoclsparse_?gthrz is a group of functions that gather the elements

indexed in indx from the dense vector y into the sparse vector x. The gathered elements in \(y\) are replaced by zero.

Let \(y\in R^m\) (or \(C^m\)) be a dense vector, \(x\) be a sparse vector from the same space and \(I_x\) be a set of indices of size \(0<\) nnz \(\le m\) described by indx, then

\[ x_i = y_{I_{x_i}}, i\in\{1,\ldots,{\bf\mathsf{nnz}}\}, \text{ and after the assignment, } y_{I_{x_i}}=0, i\in\{1,\ldots,{\bf\mathsf{nnz}}\}. \]

For double precision complex vectors use aoclsparse_zgthrz() and for single precision complex vectors use aoclsparse_cgthrz().

Note

These functions assume that the indices stored in indx are less than \(m\) without duplicate elements, and that x and indx are pointers to vectors of size at least nnz.

Parameters

nnz – [in] number of non-zero entries of \(x\). If nnz is zero, then none of the entries of vectors x, y, and indx are touched.
y – [in] pointer to dense vector \(y\) of size at least \(m\).
x – [out] pointer to sparse vector \(x\) with at least nnz non-zero elements.
indx – [in] index vector of size nnz, containing the indices of the non-zero values of \(x\). Indices should range from 0 to \(m-1\), need not be ordered. The elements in this vector are only checked for non-negativity. The user should make sure that no index is out-of-bound and that it does not contains any duplicates.

Return values

aoclsparse_status_success – the operation completed successfully
aoclsparse_status_invalid_size – nnz parameter value is negative
aoclsparse_status_invalid_pointer – at least one of the pointers y, x or indx is invalid
aoclsparse_status_invalid_index_value – at least one of the indices in indx is negative

aoclsparse_?gthrs()#

aoclsparse_status aoclsparse_sgthrs(aoclsparse_int nnz, const float *y, float *x, aoclsparse_int stride)#

aoclsparse_status aoclsparse_dgthrs(aoclsparse_int nnz, const double *y, double *x, aoclsparse_int stride)#

aoclsparse_status aoclsparse_cgthrs(aoclsparse_int nnz, const void *y, void *x, aoclsparse_int stride)#

aoclsparse_status aoclsparse_zgthrs(aoclsparse_int nnz, const void *y, void *x, aoclsparse_int stride)#

Gather elements from a dense vector using a stride and store them into a sparse vector.

The aoclsparse_?gthrs is a group of functions that gather the elements from the dense vector y using a fixed stride distance and copies them into the sparse vector x.

Let \(y\in R^m\) (or \(C^m\)) be a dense vector, \(x\) be a sparse vector from the same space and stride be a (positive) striding distance, then \( x_i = y_{\text{stride} \times i}, \quad i\in\{1,\ldots,{\bf\mathsf{nnz}}\}. \)

Parameters

nnz – [in] Number of non-zero entries of \(x\). If nnz is zero, then none of the entries of vectors x and y are accessed. Note that nnz must be such that stride \(\times\) nnz must be less or equal to \(m\).
y – [in] Pointer to dense vector \(y\) of size at least \(m\).
x – [out] Pointer to sparse vector \(x\) with at least nnz non-zero elements.
stride – [in] Striding distance used to access elements in the dense vector y. It must be such that stride \(\times\) nnz is less or equal to \(m\).

Return values

aoclsparse_status_success – the operation completed successfully.
aoclsparse_status_invalid_size – at least one of the parameters nnz or stride has a negative value.
aoclsparse_status_invalid_pointer – at least one of the pointers y, or x is invalid.

Level 2#

This module holds all sparse level 2 routines.

The sparse level 2 routines describe operations between a matrix in sparse format and a vector in dense or sparse format.

aoclsparse_?mv()#

aoclsparse_status aoclsparse_smv(aoclsparse_operation op, const float *alpha, aoclsparse_matrix A, const aoclsparse_mat_descr descr, const float *x, const float *beta, float *y)#

aoclsparse_status aoclsparse_dmv(aoclsparse_operation op, const double *alpha, aoclsparse_matrix A, const aoclsparse_mat_descr descr, const double *x, const double *beta, double *y)#

aoclsparse_status aoclsparse_cmv(aoclsparse_operation op, const aoclsparse_float_complex *alpha, aoclsparse_matrix A, const aoclsparse_mat_descr descr, const aoclsparse_float_complex *x, const aoclsparse_float_complex *beta, aoclsparse_float_complex *y)#

aoclsparse_status aoclsparse_zmv(aoclsparse_operation op, const aoclsparse_double_complex *alpha, aoclsparse_matrix A, const aoclsparse_mat_descr descr, const aoclsparse_double_complex *x, const aoclsparse_double_complex *beta, aoclsparse_double_complex *y)#

Compute sparse matrix-vector multiplication for real/complex single and double data precisions.

The aoclsparse_?mv perform sparse matrix-vector products of the form

\[ y = \alpha \, op(A) \, x + \beta \, y, \]

where, \(x\) and \(y\) are dense vectors, \(\alpha\) and \(\beta\) are scalars, and \(A\) is a sparse matrix structure. The matrix operation \(op()\) is defined as:

\[\begin{split} op(A) = \left\{ \begin{array}{ll} A, & \text{if } {\bf\mathsf{op}} = \text{aoclsparse}\_\text{operation}\_\text{none} \\ A^T, & \text{if } {\bf\mathsf{op}} = \text{aoclsparse}\_\text{operation}\_\text{transpose} \\ A^H, & \text{if } {\bf\mathsf{op}} = \text{aoclsparse}\_\text{operation}\_\text{conjugate}\_\text{transpose} \end{array} \right. \end{split}\]

Example - C++ (tests/examples/sample_spmv.cpp)

/* ************************************************************************
 * Copyright (c) 2020-2023 Advanced Micro Devices, Inc.
 *
 * Permission is hereby granted, free of charge, to any person obtaining a copy
 * of this software and associated documentation files (the "Software"), to deal
 * in the Software without restriction, including without limitation the rights
 * to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
 * copies of the Software, and to permit persons to whom the Software is
 * furnished to do so, subject to the following conditions:
 *
 * The above copyright notice and this permission notice shall be included in
 * all copies or substantial portions of the Software.
 *
 * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
 * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
 * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
 * AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
 * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
 * OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
 * THE SOFTWARE.
 *
 * ************************************************************************ */

#include "aoclsparse.h"

#include <iostream>

#define M 5
#define N 5
#define NNZ 8

int main(void)
{
    aoclsparse_operation trans = aoclsparse_operation_none;

    double alpha = 1.0;
    double beta  = 0.0;

    // Print aoclsparse version
    std::cout << aoclsparse_get_version() << std::endl;

    // Create matrix descriptor
    aoclsparse_mat_descr descr;
    // aoclsparse_create_mat_descr set aoclsparse_matrix_type to aoclsparse_matrix_type_general
    // and aoclsparse_index_base to aoclsparse_index_base_zero.
    aoclsparse_create_mat_descr(&descr);

    aoclsparse_index_base base = aoclsparse_index_base_zero;

    // Initialise matrix
    //  1  0  0  2  0
    //  0  3  0  0  0
    //  0  0  4  0  0
    //  0  5  0  6  7
    //  0  0  0  0  8
    aoclsparse_int    csr_row_ptr[M + 1] = {0, 2, 3, 4, 7, 8};
    aoclsparse_int    csr_col_ind[NNZ]   = {0, 3, 1, 2, 1, 3, 4, 4};
    double            csr_val[NNZ]       = {1, 2, 3, 4, 5, 6, 7, 8};
    aoclsparse_matrix A;
    aoclsparse_create_dcsr(&A, base, M, N, NNZ, csr_row_ptr, csr_col_ind, csr_val);

    // Initialise vectors
    double x[N] = {1.0, 2.0, 3.0, 4.0, 5.0};
    double y[M];

    //to identify hint id(which routine is to be executed, destroyed later)
    aoclsparse_set_mv_hint(A, trans, descr, 1);

    // Optimize the matrix, "A"
    aoclsparse_optimize(A);

    std::cout << "Invoking aoclsparse_dmv..";
    //Invoke SPMV API (double precision)
    aoclsparse_dmv(trans, &alpha, A, descr, x, &beta, y);
    std::cout << "Done." << std::endl;
    std::cout << "Output Vector:" << std::endl;
    for(aoclsparse_int i = 0; i < M; i++)
        std::cout << y[i] << std::endl;

    aoclsparse_destroy_mat_descr(descr);
    aoclsparse_destroy(&A);
    return 0;
}

Example - C (tests/examples/sample_spmv_c.c)

/* ************************************************************************
 * Copyright (c) 2023 Advanced Micro Devices, Inc.
 *
 * Permission is hereby granted, free of charge, to any person obtaining a copy
 * of this software and associated documentation files (the "Software"), to deal
 * in the Software without restriction, including without limitation the rights
 * to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
 * copies of the Software, and to permit persons to whom the Software is
 * furnished to do so, subject to the following conditions:
 *
 * The above copyright notice and this permission notice shall be included in
 * all copies or substantial portions of the Software.
 *
 * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
 * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
 * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
 * AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
 * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
 * OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
 * THE SOFTWARE.
 *
 * ************************************************************************ */

#include "aoclsparse.h"

#include <math.h>
#include <stdio.h>

#define M 5
#define N 5
#define NNZ 8

int main(void)
{
    aoclsparse_status     status;
    aoclsparse_matrix     A     = NULL;
    aoclsparse_mat_descr  descr = NULL;
    aoclsparse_operation  trans = aoclsparse_operation_none;
    aoclsparse_index_base base  = aoclsparse_index_base_zero;

    double alpha = 1.0;
    double beta  = 0.0;

    // Input matrix
    //  1  0  0  2  0
    //  0  3  0  0  0
    //  0  0  4  0  0
    //  0  5  0  6  7
    //  0  0  0  0  8
    aoclsparse_int csr_row_ptr[M + 1] = {0, 2, 3, 4, 7, 8};
    aoclsparse_int csr_col_ind[NNZ]   = {0, 3, 1, 2, 1, 3, 4, 4};
    double         csr_val[NNZ]       = {1, 2, 3, 4, 5, 6, 7, 8};

    // Input vectors
    double x[N] = {1.0, 2.0, 3.0, 4.0, 5.0};
    double y[M];

    // Expected reference result & tolerance
    double         y_exp[M] = {9.0, 6.0, 12.0, 69.0, 40.0};
    double         tol      = 1e-12;
    aoclsparse_int oki, ok = 1;

    // Print aoclsparse version
    printf("%s\n", aoclsparse_get_version());

    // Create matrix descriptor
    // aoclsparse_create_mat_descr sets aoclsparse_matrix_type to aoclsparse_matrix_type_general
    // and aoclsparse_index_base to aoclsparse_index_base_zero.
    aoclsparse_create_mat_descr(&descr);

    // Initialise sparse matrix
    status = aoclsparse_create_dcsr(&A, base, M, N, NNZ, csr_row_ptr, csr_col_ind, csr_val);
    if(status != aoclsparse_status_success)
    {
        printf("Error while creating a sparse matrix, status = %i\n", status);
        return 1;
    }

    // Hint the system what operation to expect // to identify hint id(which routine is to be executed, destroyed later)
    status = aoclsparse_set_mv_hint(A, trans, descr, 1);
    if(status != aoclsparse_status_success)
    {
        printf("Error while hinting operation, status = %i\n", status);
        return 1;
    }

    // Optimize the matrix
    status = aoclsparse_optimize(A);
    if(status != aoclsparse_status_success)
    {
        printf("Error while optimizing the matrix, status = %i\n", status);
        return 1;
    }

    // Invoke SPMV API (double precision)
    printf("Invoking aoclsparse_dmv...\n");
    status = aoclsparse_dmv(trans, &alpha, A, descr, x, &beta, y);
    if(status != aoclsparse_status_success)
    {
        printf("Error while computing SPMV, status = %i\n", status);
        return 1;
    }

    // Print and check the results
    printf("Output vector:\n");
    for(aoclsparse_int i = 0; i < M; i++)
    {
        oki = fabs(y[i] - y_exp[i]) <= tol;
        printf("  %lf %c\n", y[i], oki ? ' ' : '!');
        ok = ok && oki;
    }

    aoclsparse_destroy_mat_descr(descr);
    aoclsparse_destroy(&A);
    return ok ? 0 : 2;
}

Parameters

op – [in] Matrix operation, op can be one of aoclsparse_operation_none, aoclsparse_operation_conjugate_transpose, or aoclsparse_operation_conjugate_transpose.
alpha – [in] Scalar \(\alpha\).
A – [in] The sparse matrix created using e.g. aoclsparse_create_scsr() or other variant. Matrix is considered of size \(m\) by \(n\).
descr – [in] Descriptor of the matrix. These functions support the following aoclsparse_matrix_type types: aoclsparse_matrix_type_general, aoclsparse_matrix_type_triangular, aoclsparse_matrix_type_symmetric, and aoclsparse_matrix_type_hermitian. Both base-zero and base-one are supported, however, the index base needs to match with the one defined in matrix A.
x – [in] An array of n elements if \(op(A) = A\); or of m elements if \(op(A) = A^T\) or \(op(A) = A^H\).
beta – [in] Scalar \(\beta\).
y – [inout] An array of m elements if \(op(A) = A\); or of n elements if \(op(A) = A^T\) or \(op(A) = A^H\).

Return values

aoclsparse_status_success – The operation completed successfully.
aoclsparse_status_invalid_size – The value of m, n or nnz is invalid.
aoclsparse_status_invalid_pointer – descr, alpha, internal structures related to the sparse matrix A, x, beta or y has an invalid pointer.
aoclsparse_status_not_implemented – The requested functionality is not implemented.

aoclsparse_?trsv()#

aoclsparse_status aoclsparse_strsv(aoclsparse_operation trans, const float alpha, aoclsparse_matrix A, const aoclsparse_mat_descr descr, const float *b, float *x)#

aoclsparse_status aoclsparse_dtrsv(aoclsparse_operation trans, const double alpha, aoclsparse_matrix A, const aoclsparse_mat_descr descr, const double *b, double *x)#

aoclsparse_status aoclsparse_ctrsv(aoclsparse_operation trans, const aoclsparse_float_complex alpha, aoclsparse_matrix A, const aoclsparse_mat_descr descr, const aoclsparse_float_complex *b, aoclsparse_float_complex *x)#

aoclsparse_status aoclsparse_ztrsv(aoclsparse_operation trans, const aoclsparse_double_complex alpha, aoclsparse_matrix A, const aoclsparse_mat_descr descr, const aoclsparse_double_complex *b, aoclsparse_double_complex *x)#

Sparse triangular solver for real/complex single and double data precisions.

The function aoclsparse_strsv() and variants solve sparse lower (or upper) triangular linear system of equations. The system is defined by the sparse \(m \times m\) matrix \(A\), the dense solution \(m\)-vector \(x\), and the right-hand side dense \(m\)-vector \(b\). Vector \(b\) is multiplied by \(\alpha\). The solution \(x\) is estimated by solving

\[ op(L) \cdot x = \alpha \cdot b, \quad \text{ or } \quad op(U) \cdot x = \alpha \cdot b, \]

where \(L = \text{tril}(A)\) is the lower triangle of matrix \(A\), similarly, \(U = \text{triu}(A)\) is the upper triangle of matrix \(A\). The operator \(op()\) is regarded as the matrix linear operation,

\[\begin{split} op(A) = \left\{ \begin{array}{ll} A, & \text{if } {\bf\mathsf{trans}} = \text{aoclsparse}\_\text{operation}\_\text{none} \\ A^T, & \text{if } {\bf\mathsf{trans}} = \text{aoclsparse}\_\text{operation}\_\text{transpose} \\ A^H, & \text{if } {\bf\mathsf{trans}} = \text{aoclsparse}\_\text{operation}\_\text{conjugate}\_\text{transpose} \end{array} \right. \end{split}\]

Notes

1. This routine supports only sparse matrices in CSR format.

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

3. The input matrix need not be (upper or lower) triangular matrix, in the descr, the fill_mode entity specifies which triangle to consider, namely, if fill_mode = aoclsparse_fill_mode_lower, then

\[ op(L) \cdot x = \alpha \cdot b, \]

otherwise, if fill_mode = aoclsparse_fill_mode_upper, then

\[ op(U) \cdot x = \alpha \cdot b, \]

is solved.

4. To increase performance and if the matrix \(A\) is to be used more than once to solve for different right-hand sides b, then it is encouraged to provide hints using aoclsparse_set_sv_hint() and aoclsparse_optimize(), otherwise, the optimization for the matrix will be done by the solver on entry.

5. There is a kid (Kernel ID) variation of TRSV, namely with a suffix of _kid, aoclsparse_strsv_kid() (and variations) where it is possible to specify the TRSV kernel to use (if possible).

Example - Real space (tests/examples/sample_dtrsv.cpp)

/* ************************************************************************
 * Copyright (c) 2022-2023 Advanced Micro Devices, Inc.
 *
 * Permission is hereby granted, free of charge, to any person obtaining a copy
 * of this software and associated documentation files (the "Software"), to deal
 * in the Software without restriction, including without limitation the rights
 * to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
 * copies of the Software, and to permit persons to whom the Software is
 * furnished to do so, subject to the following conditions:
 *
 * The above copyright notice and this permission notice shall be included in
 * all copies or substantial portions of the Software.
 *
 * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
 * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
 * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
 * AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
 * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
 * OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
 * THE SOFTWARE.
 *
 * ************************************************************************ */

#include "aoclsparse.h"

#include <iomanip>
#include <iostream>

int main()
{
    std::cout << "-------------------------------" << std::endl
              << " Triangle solve sample program" << std::endl
              << "-------------------------------" << std::endl
              << std::endl;

    // Create a tri-diagonal matrix in CSR format
    // | 1  2  0  0 |
    // | 3  1  2  0 |
    // | 0  3  1  2 |
    // | 0  0  3  1 |
    aoclsparse_int    n = 4, m = 4, nnz = 10;
    aoclsparse_int    icrow[5] = {0, 2, 5, 8, 10};
    aoclsparse_int    icol[18] = {0, 1, 0, 1, 2, 1, 2, 3, 2, 3};
    double            aval[18] = {1, 2, 3, 1, 2, 3, 1, 2, 3, 1};
    aoclsparse_status status;

    // create aoclsparse matrix and its descriptor
    aoclsparse_matrix     A;
    aoclsparse_index_base base = aoclsparse_index_base_zero;
    aoclsparse_mat_descr  descr_a;
    aoclsparse_operation  trans = aoclsparse_operation_none;
    aoclsparse_create_dcsr(&A, base, m, n, nnz, icrow, icol, aval);
    aoclsparse_create_mat_descr(&descr_a);

    // A = L + D + U where:
    // - L and U are the strict lower and upper triangle of A
    // - D is the diagonal
    // Use the matrix descriptor to solve (L+D)x = b
    // b = [1, 4, 4, 4]^t
    double b[4] = {1, 4, 4, 4};
    aoclsparse_set_mat_type(descr_a, aoclsparse_matrix_type_triangular);
    aoclsparse_set_mat_fill_mode(descr_a, aoclsparse_fill_mode_lower);
    aoclsparse_set_mat_index_base(descr_a, base);
    status = aoclsparse_set_sv_hint(A, trans, descr_a, 1);
    if(status != aoclsparse_status_success)
    {
        std::cerr << "Error setting SV hint, status = " << status << std::endl;
        return 1;
    }
    status = aoclsparse_optimize(A);
    if(status != aoclsparse_status_success)
    {
        std::cerr << "Error returned from aoclsparse_optimize, status = " << status << "."
                  << std::endl;
        return 2;
    }
    // Call the triangle solver
    double alpha = 1.0;
    double x[n];
    status = aoclsparse_dtrsv(trans, alpha, A, descr_a, b, x);
    if(status != aoclsparse_status_success)
    {
        std::cerr << "Error returned from aoclsparse_dtrsv, status = " << status << "."
                  << std::endl;
        return 3;
    }

    // Print the result
    std::cout << "Solving (L+D)x = b: " << std::endl << "  x = ";
    std::cout << std::fixed;
    std::cout.precision(1);
    for(aoclsparse_int i = 0; i < n; i++)
        std::cout << x[i] << "  ";
    std::cout << std::endl;

    // The same method can be used to solve (U+D)x = b
    aoclsparse_set_mat_fill_mode(descr_a, aoclsparse_fill_mode_upper);
    status = aoclsparse_dtrsv(trans, alpha, A, descr_a, b, x);
    if(status != aoclsparse_status_success)
    {
        std::cerr << "Error returned from aoclsparse_dtrsv, status = " << status << "."
                  << std::endl;
        return 3;
    }

    // Print the result
    std::cout << std::endl;
    std::cout << "Solving (U+D)x = b: " << std::endl << "  x = ";
    for(aoclsparse_int i = 0; i < n; i++)
        std::cout << x[i] << "  ";
    std::cout << std::endl;

    // destroy the aoclsparse memory
    aoclsparse_destroy_mat_descr(descr_a);
    aoclsparse_destroy(&A);

    return 0;
}

Example - Complex space (tests/examples/sample_ztrsv.cpp)

/* ************************************************************************
 * Copyright (c) 2023-2024 Advanced Micro Devices, Inc.
 *
 * Permission is hereby granted, free of charge, to any person obtaining a copy
 * of this software and associated documentation files (the "Software"), to deal
 * in the Software without restriction, including without limitation the rights
 * to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
 * copies of the Software, and to permit persons to whom the Software is
 * furnished to do so, subject to the following conditions:
 *
 * The above copyright notice and this permission notice shall be included in
 * all copies or substantial portions of the Software.
 *
 * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
 * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
 * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
 * AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
 * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
 * OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
 * THE SOFTWARE.
 *
 * ************************************************************************ */

#include "aoclsparse.h"

#include <assert.h>
#include <complex>
#include <iomanip>
#include <iostream>
#include <limits>
#include <vector>

int main()
{
    std::cout << "-------------------------------" << std::endl
              << " Triangle solve sample program" << std::endl
              << "-------------------------------" << std::endl
              << std::endl;

    /* Solve two complex linear systems of equations
     * Lx = alpha b, and U^Tx = alpha b,
     * where L is the lower triangular part of A and
     * U is the upper triangular part of A. U^H is
     * the conjugate transpose of U (a lower triangular matrix).
     * The complex matrix A is the tri-diagonal matrix in CSR format
     *
     * | 1+3i  2+5i     0     0 |
     * | 3     1+2i  2-2i     0 |
     * |    0     3     1  2+3i |
     * |    0     0  3+1i  1+1i |
     */
    const aoclsparse_int                   n = 4, m = 4, nnz = 10;
    aoclsparse_int                         icrow[5] = {0, 2, 5, 8, 10};
    aoclsparse_int                         icol[18] = {0, 1, 0, 1, 2, 1, 2, 3, 2, 3};
    std::vector<aoclsparse_double_complex> aval     = {{1., 3.},
                                                       {2., 5.},
                                                       {3., 0.},
                                                       {1., 2.},
                                                       {2., -2.},
                                                       {3., 0.},
                                                       {1., 0.},
                                                       {2., 3.},
                                                       {3., 1.},
                                                       {1., 1.}};
    aoclsparse_status                      status;
    std::vector<aoclsparse_double_complex> ref;

    // create aoclsparse matrix and its descriptor
    aoclsparse_matrix     A;
    aoclsparse_index_base base = aoclsparse_index_base_zero;
    aoclsparse_mat_descr  descr_a;
    aoclsparse_operation  trans = aoclsparse_operation_none;
    status = aoclsparse_create_zcsr(&A, base, m, n, nnz, icrow, icol, aval.data());
    if(status != aoclsparse_status_success)
    {
        std::cerr << "Error creating the matrix, status = " << status << std::endl;
        return 1;
    }
    aoclsparse_create_mat_descr(&descr_a);

    /* Solve the lower triangular system Lx = b,
     * here alpha=1 and b = [1+i, 4+2i, 4+i, 4].
     */
    std::vector<aoclsparse_double_complex> b = {{1., 1.}, {4., 2.}, {4., 1}, {4., 0}};
    aoclsparse_set_mat_type(descr_a, aoclsparse_matrix_type_triangular);
    aoclsparse_set_mat_fill_mode(descr_a, aoclsparse_fill_mode_lower);
    status = aoclsparse_set_sv_hint(A, trans, descr_a, 1);
    if(status != aoclsparse_status_success)
    {
        std::cerr << "Error setting SV hint, status = " << status << std::endl;
        return 1;
    }
    status = aoclsparse_optimize(A);
    if(status != aoclsparse_status_success)
    {
        std::cerr << "Error returned from aoclsparse_optimize, status = " << status << "."
                  << std::endl;
        return 2;
    }

    // Call the triangle solver
    aoclsparse_double_complex alpha = {1.0, 0.};
    aoclsparse_double_complex x[n];
    status = aoclsparse_ztrsv(trans, alpha, A, descr_a, b.data(), x);
    if(status != aoclsparse_status_success)
    {
        std::cerr << "Error returned from aoclsparse_ztrsv, status = " << status << "."
                  << std::endl;
        return 3;
    }

    // Print and check the result
    double tol = 20 * std::numeric_limits<double>::epsilon();
    ref.assign({{0.4, -0.2}, {1.6, -0.6}, {-0.8, 2.8}, {0.8, -8.4}});
    std::cout << "Solving Lx = alpha b: " << std::endl << "  x = ";
    std::cout << std::fixed;
    std::cout.precision(1);
    bool oki, ok = true;
    for(aoclsparse_int i = 0; i < n; i++)
    {
        oki = std::abs(x[i].real - ref[i].real) <= tol && std::abs(x[i].imag - ref[i].imag) <= tol;
        std::cout << "(" << x[i].real << ", " << x[i].imag << "i) " << (oki ? " " : "!  ");
        ok &= oki;
    }
    std::cout << std::endl;

    /* Solve the lower triangular system U^Hx = b,
     * here alpha=1-i and b is unchanged.
     */
    alpha = {1, -1};
    // Indicate to use only the conjugate transpose of the upper part of A
    trans = aoclsparse_operation_conjugate_transpose;
    aoclsparse_set_mat_fill_mode(descr_a, aoclsparse_fill_mode_upper);
    status = aoclsparse_ztrsv(trans, alpha, A, descr_a, b.data(), x);
    if(status != aoclsparse_status_success)
    {
        std::cerr << "Error returned from aoclsparse_ztrsv, status = " << status << "."
                  << std::endl;
        return 3;
    }

    // Print and check the result
    ref.assign({{0.2, 0.6}, {1.4, 0.6}, {3.4, -7.0}, {-1.0, 19.2}});
    std::cout << std::endl;
    std::cout << "Solving U^Hx = alpha*b: " << std::endl << "  x = ";
    for(aoclsparse_int i = 0; i < n; i++)
    {
        oki = std::abs(x[i].real - ref[i].real) <= tol && std::abs(x[i].imag - ref[i].imag) <= tol;
        std::cout << "(" << x[i].real << ", " << x[i].imag << "i) " << (oki ? " " : "!  ");
        ok &= oki;
    }
    std::cout << std::endl;

    // Destroy the aoclsparse memory
    aoclsparse_destroy_mat_descr(descr_a);
    aoclsparse_destroy(&A);

    return ok ? 0 : 4;
}

Parameters

trans – [in] matrix operation type, either aoclsparse_operation_none, aoclsparse_operation_transpose, or aoclsparse_operation_conjugate_transpose.
alpha – [in] scalar \(\alpha\), used to multiply right-hand side vector \(b\).
A – [inout] matrix containing data used to represent the \(m \times m\) triangular linear system to solve.
descr – [in] matrix descriptor. Supported matrix types are aoclsparse_matrix_type_symmetric and aoclsparse_matrix_type_triangular.
b – [in] array of m elements, storing the right-hand side.
x – [out] array of m elements, storing the solution if solver returns aoclsparse_status_success.

Return values

aoclsparse_status_success – the operation completed successfully and \(x\) contains the solution to the linear system of equations.
aoclsparse_status_invalid_size – matrix \(A\) or \(op(A)\) is invalid.
aoclsparse_status_invalid_pointer – One or more of A, descr, x, b are invalid pointers.
aoclsparse_status_internal_error – an internal error occurred.
aoclsparse_status_not_implemented – the requested operation is not yet implemented.
other – possible failure values from a call to aoclsparse_optimize.

aoclsparse_status aoclsparse_strsv_strided(aoclsparse_operation trans, const float alpha, aoclsparse_matrix A, const aoclsparse_mat_descr descr, const float *b, const aoclsparse_int incb, float *x, const aoclsparse_int incx)#

aoclsparse_status aoclsparse_dtrsv_strided(aoclsparse_operation trans, const double alpha, aoclsparse_matrix A, const aoclsparse_mat_descr descr, const double *b, const aoclsparse_int incb, double *x, const aoclsparse_int incx)#

aoclsparse_status aoclsparse_ctrsv_strided(aoclsparse_operation trans, const aoclsparse_float_complex alpha, aoclsparse_matrix A, const aoclsparse_mat_descr descr, const aoclsparse_float_complex *b, const aoclsparse_int incb, aoclsparse_float_complex *x, const aoclsparse_int incx)#

aoclsparse_status aoclsparse_ztrsv_strided(aoclsparse_operation trans, const aoclsparse_double_complex alpha, aoclsparse_matrix A, const aoclsparse_mat_descr descr, const aoclsparse_double_complex *b, const aoclsparse_int incb, aoclsparse_double_complex *x, const aoclsparse_int incx)#

This is a variation of TRSV, namely with a suffix of _strided, allows to set the stride for the dense vectors b and x.

For full details refer to aoclsparse_?trsv().

Parameters

trans – [in] matrix operation type, either aoclsparse_operation_none, aoclsparse_operation_transpose, or aoclsparse_operation_conjugate_transpose.
alpha – [in] scalar \(\alpha\), used to multiply right-hand side vector \(b\).
A – [inout] matrix containing data used to represent the \(m \times m\) triangular linear system to solve.
descr – [in] matrix descriptor. Supported matrix types are aoclsparse_matrix_type_symmetric and aoclsparse_matrix_type_triangular.
b – [in] array of m elements, storing the right-hand side.
incb – [in] a positive integer holding the stride value for \(b\) vector.
x – [out] array of m elements, storing the solution if solver returns aoclsparse_status_success.
incx – [in] a positive integer holding the stride value for \(x\) vector.

aoclsparse_status aoclsparse_strsv_kid(aoclsparse_operation trans, const float alpha, aoclsparse_matrix A, const aoclsparse_mat_descr descr, const float *b, float *x, aoclsparse_int kid)#

aoclsparse_status aoclsparse_dtrsv_kid(aoclsparse_operation trans, const double alpha, aoclsparse_matrix A, const aoclsparse_mat_descr descr, const double *b, double *x, aoclsparse_int kid)#

aoclsparse_status aoclsparse_ctrsv_kid(aoclsparse_operation trans, const aoclsparse_float_complex alpha, aoclsparse_matrix A, const aoclsparse_mat_descr descr, const aoclsparse_float_complex *b, aoclsparse_float_complex *x, aoclsparse_int kid)#

aoclsparse_status aoclsparse_ztrsv_kid(aoclsparse_operation trans, const aoclsparse_double_complex alpha, aoclsparse_matrix A, const aoclsparse_mat_descr descr, const aoclsparse_double_complex *b, aoclsparse_double_complex *x, aoclsparse_int kid)#

Sparse triangular solver for real/complex single and double data precisions (kernel flag variation).

For full details refer to aoclsparse_?trsv().

This variation of TRSV, namely with a suffix of _kid, allows to choose which TRSV kernel to use (if possible). Currently the possible choices are:

kid=0: Reference implementation (No explicit AVX instructions).
kid=1: Alias to kid=2 (Kernel Template AVX 256-bit implementation)
kid=2: Kernel Template version using AVX2 extensions.
kid=3: Kernel Template version using AVX512F+ CPU extensions.

Any other Kernel ID value will default to kid = 0.

Parameters

trans – [in] matrix operation type, either aoclsparse_operation_none, aoclsparse_operation_transpose, or aoclsparse_operation_conjugate_transpose.
alpha – [in] scalar \(\alpha\), used to multiply right-hand side vector \(b\).
A – [inout] matrix containing data used to represent the \(m \times m\) triangular linear system to solve.
descr – [in] matrix descriptor. Supported matrix types are aoclsparse_matrix_type_symmetric and aoclsparse_matrix_type_triangular.
b – [in] array of m elements, storing the right-hand side.
x – [out] array of m elements, storing the solution if solver returns aoclsparse_status_success.
kid – [in] Kernel ID, hints a request on which TRSV kernel to use.

aoclsparse_?dotmv()#

aoclsparse_status aoclsparse_sdotmv(const aoclsparse_operation op, const float alpha, aoclsparse_matrix A, const aoclsparse_mat_descr descr, const float *x, const float beta, float *y, float *d)#

aoclsparse_status aoclsparse_ddotmv(const aoclsparse_operation op, const double alpha, aoclsparse_matrix A, const aoclsparse_mat_descr descr, const double *x, const double beta, double *y, double *d)#

aoclsparse_status aoclsparse_cdotmv(const aoclsparse_operation op, const aoclsparse_float_complex alpha, aoclsparse_matrix A, const aoclsparse_mat_descr descr, const aoclsparse_float_complex *x, const aoclsparse_float_complex beta, aoclsparse_float_complex *y, aoclsparse_float_complex *d)#

aoclsparse_status aoclsparse_zdotmv(const aoclsparse_operation op, const aoclsparse_double_complex alpha, aoclsparse_matrix A, const aoclsparse_mat_descr descr, const aoclsparse_double_complex *x, const aoclsparse_double_complex beta, aoclsparse_double_complex *y, aoclsparse_double_complex *d)#

Performs sparse matrix-vector multiplication followed by vector-vector multiplication.

aoclsparse_?dotmv multiplies the scalar \(\alpha\) with a sparse \(m \times n\) matrix, defined in a sparse storage format, and the dense vector \(x\) and adds the result to the dense vector \(y\) that is multiplied by the scalar \(\beta\), such that

\[\begin{split} y := \alpha \, op(A) \, x + \beta \, y, \quad \text{ with } \quad op(A) = \left\{ \begin{array}{ll} A, & \text{if } {\bf\mathsf{op}} = \text{aoclsparse}\_\text{operation}\_\text{none} \\ A^T, & \text{if } {\bf\mathsf{op}} = \text{aoclsparse}\_\text{operation}\_\text{transpose} \\ A^H, & \text{if } {\bf\mathsf{op}} = \text{aoclsparse}\_\text{operation}\_\text{conjugate}\_\text{transpose} \end{array} \right. \end{split}\]

followed by dot product of dense vectors \(x\) and \(y\) such that

\[\begin{split} {\bf\mathsf{d}} = \left\{ \begin{array}{ll} \sum_{i=0}^{\min(m,n)-1} x_{i} \; y_{i}, & \text{real case} \\ \sum_{i=0}^{\min(m,n)-1} \overline{x_i} \; y_{i}, & \text{complex case} \end{array} \right. \end{split}\]

Example (tests/examples/sample_dotmv.cpp)

/* ************************************************************************
 * Copyright (c) 2023-2024 Advanced Micro Devices, Inc.
 *
 * Permission is hereby granted, free of charge, to any person obtaining a copy
 * of this software and associated documentation files (the "Software"), to deal
 * in the Software without restriction, including without limitation the rights
 * to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
 * copies of the Software, and to permit persons to whom the Software is
 * furnished to do so, subject to the following conditions:
 *
 * The above copyright notice and this permission notice shall be included in
 * all copies or substantial portions of the Software.
 *
 * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
 * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
 * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
 * AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
 * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
 * OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
 * THE SOFTWARE.
 *
 * ************************************************************************ */

#include "aoclsparse.h"

#include <iostream>
#include <limits>

#define M 6
#define N 5
#define NNZ 8

int main(void)
{
    aoclsparse_operation op = aoclsparse_operation_none;

    double alpha = 1.0;
    double beta  = 0.0;

    // Print aoclsparse version
    std::cout << aoclsparse_get_version() << std::endl;

    // Create matrix descriptor
    aoclsparse_mat_descr descr;
    // aoclsparse_create_mat_descr set aoclsparse_matrix_type to aoclsparse_matrix_type_general
    // and aoclsparse_index_base to aoclsparse_index_base_zero.
    aoclsparse_create_mat_descr(&descr);

    aoclsparse_index_base base = aoclsparse_index_base_zero;

    // Initialise matrix
    //  1  0  0  2  0
    //  0  3  0  0  0
    //  0  0  4  0  0
    //  0  5  0  6  7
    //  0  0  0  0  8
    //  0  0  0  0  0
    aoclsparse_int    csr_row_ptr[M + 1] = {0, 2, 3, 4, 7, 8, 8};
    aoclsparse_int    csr_col_ind[NNZ]   = {0, 3, 1, 2, 1, 3, 4, 4};
    double            csr_val[NNZ]       = {1, 2, 3, 4, 5, 6, 7, 8};
    aoclsparse_matrix A;
    aoclsparse_create_dcsr(&A, base, M, N, NNZ, csr_row_ptr, csr_col_ind, csr_val);

    // Initialise vectors
    double x[N]     = {1.0, 1.0, 1.0, 1.0, 1.0};
    double y[M]     = {1.0, 1.0, 1.0, 1.0, 1.0, 1.0};
    double d        = 0;
    double y_exp[M] = {3.0, 3.0, 4.0, 18.0, 8.0};
    double d_exp    = 36;

    //to identify hint id(which routine is to be executed, destroyed later)
    aoclsparse_set_dotmv_hint(A, op, descr, 1);

    // Optimize the matrix, "A"
    aoclsparse_optimize(A);

    std::cout << "Invoking aoclsparse_ddotmv..";
    //Invoke DOTMV API (double precision)
    aoclsparse_ddotmv(op, alpha, A, descr, x, beta, y, &d);
    std::cout << "Done." << std::endl;

    // Print and check the result
    std::cout << std::endl << "y = ";
    bool   oki, ok = true;
    double tol = 20 * std::numeric_limits<double>::epsilon();
    for(aoclsparse_int i = 0; i < M; i++)
    {
        oki = std::abs(y[i] - y_exp[i]) <= tol;
        std::cout << y[i] << (oki ? " " : "!  ");
        ok &= oki;
    }

    ok &= (std::abs(d - d_exp) <= tol);
    std::cout << std::endl << "Output dot product: " << d << (ok ? " " : "!  ") << std::endl;

    aoclsparse_destroy_mat_descr(descr);
    aoclsparse_destroy(&A);
    return ok ? 0 : 4;
}

Parameters

op – [in] matrix operation type.
alpha – [in] scalar \(\alpha\).
A – [in] the sparse \(m \times n\) matrix structure that is created using aoclsparse_create_scsr() or other variation.
descr – [in] descriptor of the sparse CSR matrix. Both base-zero and base-one are supported, however, the index base needs to match the one used when aoclsparse_matrix was created.
x – [in] array of atleast n elements if \(op(A) = A\) or at least m elements if \(op(A) = A^T\) or \(A^H\).
beta – [in] scalar \(\beta\).
y – [inout] array of atleast m elements if \(op(A) = A\) or at least n elements if \(op(A) = A^T\) or \(A^H\).
d – [out] dot product of y and x.

Return values

aoclsparse_status_success – the operation completed successfully.
aoclsparse_status_invalid_size – m, n or nnz is invalid.
aoclsparse_status_invalid_value – (base index is neither aoclsparse_index_base_zero nor aoclsparse_index_base_one, or matrix base index and descr base index values do not match.
aoclsparse_status_invalid_pointer – descr, internal structures related to the sparse matrix A, x, y or d are invalid pointer.
aoclsparse_status_wrong_type – matrix data type is not supported.
aoclsparse_status_not_implemented – aoclsparse_matrix_type is aoclsparse_matrix_type_hermitian or, aoclsparse_matrix_format_type is not aoclsparse_csr_mat

aoclsparse_?ellmv()#

aoclsparse_status aoclsparse_sellmv(aoclsparse_operation trans, const float *alpha, aoclsparse_int m, aoclsparse_int n, aoclsparse_int nnz, const float *ell_val, const aoclsparse_int *ell_col_ind, aoclsparse_int ell_width, const aoclsparse_mat_descr descr, const float *x, const float *beta, float *y)#

aoclsparse_status aoclsparse_dellmv(aoclsparse_operation trans, const double *alpha, aoclsparse_int m, aoclsparse_int n, aoclsparse_int nnz, const double *ell_val, const aoclsparse_int *ell_col_ind, aoclsparse_int ell_width, const aoclsparse_mat_descr descr, const double *x, const double *beta, double *y)#

Real single and double precision sparse matrix vector product using ELL storage format.

aoclsparse_?ellmv multiplies the scalar \(\alpha\) with a sparse \(m \times n\) matrix, defined in ELL storage format, and the dense vector \(x\) and adds the result to the dense vector \(y\) that is multiplied by the scalar \(\beta\), such that

\[ y = \alpha \, op(A) \, x + \beta \, y, \]

with

Note

Currently, only trans = aoclsparse_operation_none is supported.

Parameters

trans – [in] matrix operation type.
alpha – [in] scalar \(\alpha\).
m – [in] number of rows of the sparse ELL matrix.
n – [in] number of columns of the sparse ELL matrix.
nnz – [in] number of non-zero entries of the sparse ELL matrix.
descr – [in] descriptor of the sparse ELL matrix. Both, base-zero and base-one input arrays of ELL matrix are supported
ell_val – [in] array that contains the elements of the sparse ELL matrix. Padded elements should be zero.
ell_col_ind – [in] array that contains the column indices of the sparse ELL matrix. Padded column indices should be -1.
ell_width – [in] number of non-zero elements per row of the sparse ELL matrix.
x – [in] array of n elements ( \(op(A) = A\)) or m elements ( \(op(A) = A^T\) or \(op(A) = A^H\)).
beta – [in] scalar \(\beta\).
y – [inout] array of m elements ( \(op(A) = A\)) or n elements ( \(op(A) = A^T\) or \(op(A) = A^H\)).

Return values

aoclsparse_status_success – the operation completed successfully.
aoclsparse_status_invalid_size – m, n or ell_width is invalid.
aoclsparse_status_invalid_pointer – descr, alpha, ell_val, ell_col_ind, x, beta or y pointer is invalid.
aoclsparse_status_not_implemented – trans is not aoclsparse_operation_none, or aoclsparse_matrix_type is not aoclsparse_matrix_type_general.

aoclsparse_?diamv()#

aoclsparse_status aoclsparse_sdiamv(aoclsparse_operation trans, const float *alpha, aoclsparse_int m, aoclsparse_int n, aoclsparse_int nnz, const float *dia_val, const aoclsparse_int *dia_offset, aoclsparse_int dia_num_diag, const aoclsparse_mat_descr descr, const float *x, const float *beta, float *y)#

aoclsparse_status aoclsparse_ddiamv(aoclsparse_operation trans, const double *alpha, aoclsparse_int m, aoclsparse_int n, aoclsparse_int nnz, const double *dia_val, const aoclsparse_int *dia_offset, aoclsparse_int dia_num_diag, const aoclsparse_mat_descr descr, const double *x, const double *beta, double *y)#

Real single and double precision sparse matrix vector product using DIA storage format.

aoclsparse_?diamv multiplies the scalar \(\alpha\) with a sparse \(m \times n\) matrix, defined in DIA storage format, and the dense vector \(x\) and adds the result to the dense vector \(y\) that is multiplied by the scalar \(\beta\), such that

\[ y = \alpha \, op(A) \, x + \beta \, y, \]

with

Note

Currently, only trans = aoclsparse_operation_none is supported.

Parameters

trans – [in] matrix operation type.
alpha – [in] scalar \(\alpha\).
m – [in] number of rows of the matrix.
n – [in] number of columns of the matrix.
nnz – [in] number of non-zero entries of the matrix.
descr – [in] descriptor of the sparse DIA matrix.
dia_val – [in] array that contains the elements of the matrix. Padded elements should be zero.
dia_offset – [in] array that contains the offsets of each diagonal of the matrix.
dia_num_diag – [in] number of diagonals in the matrix.
x – [in] array of n elements ( \(op(A) = A\)) or m elements ( \(op(A) = A^T\) or \(op(A) = A^H\)).
beta – [in] scalar \(\beta\).
y – [inout] array of m elements ( \(op(A) = A\)) or n elements ( \(op(A) = A^T\) or \(op(A) = A^H\)).

Return values

aoclsparse_status_success – the operation completed successfully.
aoclsparse_status_invalid_size – m, n or ell_width is invalid.
aoclsparse_status_invalid_pointer – descr, alpha, ell_val, ell_col_ind, x, beta or y pointer is invalid.
aoclsparse_status_not_implemented – trans is not aoclsparse_operation_none, or aoclsparse_matrix_type is not aoclsparse_matrix_type_general.

aoclsparse_?bsrmv()#

aoclsparse_status aoclsparse_sbsrmv(aoclsparse_operation trans, const float *alpha, aoclsparse_int mb, aoclsparse_int nb, aoclsparse_int bsr_dim, const float *bsr_val, const aoclsparse_int *bsr_col_ind, const aoclsparse_int *bsr_row_ptr, const aoclsparse_mat_descr descr, const float *x, const float *beta, float *y)#

aoclsparse_status aoclsparse_dbsrmv(aoclsparse_operation trans, const double *alpha, aoclsparse_int mb, aoclsparse_int nb, aoclsparse_int bsr_dim, const double *bsr_val, const aoclsparse_int *bsr_col_ind, const aoclsparse_int *bsr_row_ptr, const aoclsparse_mat_descr descr, const double *x, const double *beta, double *y)#

Real single and double precision matrix vector product using BSR storage format.

aoclsparse_?bsrmv multiplies the scalar \(\alpha\) with a sparse mb times bsr_dim by nb times bsr_dim matrix, defined in BSR storage format, and the dense vector \(x\) and adds the result to the dense vector \(y\) that is multiplied by the scalar \(\beta\), such that

\[ y = \alpha \, op(A) \, x + \beta \, y, \]

with

Note

Only trans = aoclsparse_operation_none is supported.

Parameters

trans – [in] matrix operation type.
mb – [in] number of block rows of the sparse BSR matrix.
nb – [in] number of block columns of the sparse BSR matrix.
alpha – [in] scalar \(\alpha\).
descr – [in] descriptor of the sparse BSR matrix. Both, base-zero and base-one input arrays of BSR matrix are supported.
bsr_val – [in] array of nnzb blocks of the sparse BSR matrix.
bsr_row_ptr – [in] array of mb+1 elements that point to the start of every block row of the sparse BSR matrix.
bsr_col_ind – [in] array of nnz containing the block column indices of the sparse BSR matrix.
bsr_dim – [in] block dimension of the sparse BSR matrix.
x – [in] array of nb times bsr_dim elements ( \(op(A) = A\)) or mb times bsr_dim elements ( \(op(A) = A^T\) or \(op(A) = A^H\)).
beta – [in] scalar \(\beta\).
y – [inout] array of mb times bsr_dim elements ( \(op(A) = A\)) or nb times bsr_dim elements ( \(op(A) = A^T\) or \(op(A) = A^H\)).

Return values

aoclsparse_status_success – the operation completed successfully.
aoclsparse_status_invalid_handle – the library context was not initialized.
aoclsparse_status_invalid_size – mb, nb, nnzb or bsr_dim is invalid.
aoclsparse_status_invalid_pointer – descr, alpha, bsr_val, bsr_row_ind, bsr_col_ind, x, beta or y pointer is invalid.
aoclsparse_status_arch_mismatch – the device is not supported.
aoclsparse_status_not_implemented – trans is not aoclsparse_operation_none, or aoclsparse_matrix_type is not aoclsparse_matrix_type_general.

aoclsparse_?csrmv()#

aoclsparse_status aoclsparse_scsrmv(aoclsparse_operation trans, const float *alpha, aoclsparse_int m, aoclsparse_int n, aoclsparse_int nnz, const float *csr_val, const aoclsparse_int *csr_col_ind, const aoclsparse_int *csr_row_ptr, const aoclsparse_mat_descr descr, const float *x, const float *beta, float *y)#

aoclsparse_status aoclsparse_dcsrmv(aoclsparse_operation trans, const double *alpha, aoclsparse_int m, aoclsparse_int n, aoclsparse_int nnz, const double *csr_val, const aoclsparse_int *csr_col_ind, const aoclsparse_int *csr_row_ptr, const aoclsparse_mat_descr descr, const double *x, const double *beta, double *y)#

Real single and double precision sparse matrix-vector multiplication using CSR storage format.

aoclsparse_?csrmv multiplies the scalar \(\alpha\) with a sparse \(m \times n\) matrix, defined in CSR storage format, and the dense vector \(x\) and adds the result to the dense vector \(y\) that is multiplied by the scalar \(\beta\), such that

\[ y = \alpha \, op(A) \, x + \beta \, y, \]

with

Parameters

trans – [in] matrix operation type.
alpha – [in] scalar \(\alpha\).
m – [in] number of rows of the sparse CSR matrix.
n – [in] number of columns of the sparse CSR matrix.
nnz – [in] number of non-zero entries of the sparse CSR matrix.
csr_val – [in] array of nnz elements of the sparse CSR matrix.
csr_col_ind – [in] array of nnz elements containing the column indices of the sparse CSR matrix.
csr_row_ptr – [in] array of m +1 elements that point to the start of every row of the sparse CSR matrix.
descr – [in] descriptor of the sparse CSR matrix. Currently, only aoclsparse_matrix_type_general and aoclsparse_matrix_type_symmetric is supported.
x – [in] array of n elements ( \(op(A) = A\)) or m elements ( \(op(A) = A^T\) or \(op(A) = A^H\)).
beta – [in] scalar \(\beta\).
y – [inout] array of m elements ( \(op(A) = A\)) or n elements ( \(op(A) = A^T\) or \(op(A) = A^H\)).

Return values

aoclsparse_status_success – the operation completed successfully.
aoclsparse_status_invalid_size – m, n or nnz is invalid.
aoclsparse_status_invalid_pointer – descr, alpha, csr_val, csr_row_ptr, csr_col_ind, x, beta or y pointer is invalid.
aoclsparse_status_not_implemented – trans is not aoclsparse_operation_none and trans is not aoclsparse_operation_transpose. aoclsparse_matrix_type is not aoclsparse_matrix_type_general, or aoclsparse_matrix_type is not aoclsparse_matrix_type_symmetric.

aoclsparse_?csrsv()#

aoclsparse_status aoclsparse_scsrsv(aoclsparse_operation trans, const float *alpha, aoclsparse_int m, const float *csr_val, const aoclsparse_int *csr_col_ind, const aoclsparse_int *csr_row_ptr, const aoclsparse_mat_descr descr, const float *x, float *y)#

aoclsparse_status aoclsparse_dcsrsv(aoclsparse_operation trans, const double *alpha, aoclsparse_int m, const double *csr_val, const aoclsparse_int *csr_col_ind, const aoclsparse_int *csr_row_ptr, const aoclsparse_mat_descr descr, const double *x, double *y)#

Sparse triangular solve using CSR storage format for single and double data precisions.

Deprecated:: This API is superseded by aoclsparse_strsv() and aoclsparse_dtrsv().

aoclsparse_?csrsv solves a sparse triangular linear system of a sparse \(m \times m\) matrix, defined in CSR storage format, a dense solution vector \(y\) and the right-hand side \(x\) that is multiplied by \(\alpha\), such that

\[ op(A) \, y = \alpha \, x, \]

with

Note

Only trans = aoclsparse_operation_none is supported.

Note

The input matrix has to be sparse upper or lower triangular matrix with unit or non-unit main diagonal. Matrix has to be sorted. No diagonal element can be omitted from a sparse storage if the solver is called with the non-unit indicator.

Parameters

trans – [in] matrix operation type.
alpha – [in] scalar \(\alpha\).
m – [in] number of rows of the sparse CSR matrix.
csr_val – [in] array of nnz elements of the sparse CSR matrix.
csr_row_ptr – [in] array of m+1 elements that point to the start of every row of the sparse CSR matrix.
csr_col_ind – [in] array of nnz elements containing the column indices of the sparse CSR matrix.
descr – [in] descriptor of the sparse CSR matrix.
x – [in] array of m elements, holding the right-hand side.
y – [out] array of m elements, holding the solution.

Return values

aoclsparse_status_success – the operation completed successfully.
aoclsparse_status_invalid_size – m is invalid.
aoclsparse_status_invalid_pointer – descr, alpha, csr_val, csr_row_ptr, csr_col_ind, x or y pointer is invalid.
aoclsparse_status_internal_error – an internal error occurred.
aoclsparse_status_not_implemented – trans = aoclsparse_operation_conjugate_transpose or trans = aoclsparse_operation_transpose or aoclsparse_matrix_type is not aoclsparse_matrix_type_general.

Level 3#

This module holds all sparse level 3 routines.

The sparse level 3 routines describe operations between matrices.

aoclsparse_?trsm()#

aoclsparse_status aoclsparse_strsm(const aoclsparse_operation trans, const float alpha, aoclsparse_matrix A, const aoclsparse_mat_descr descr, aoclsparse_order order, const float *B, aoclsparse_int n, aoclsparse_int ldb, float *X, aoclsparse_int ldx)#

aoclsparse_status aoclsparse_dtrsm(const aoclsparse_operation trans, const double alpha, aoclsparse_matrix A, const aoclsparse_mat_descr descr, aoclsparse_order order, const double *B, aoclsparse_int n, aoclsparse_int ldb, double *X, aoclsparse_int ldx)#

aoclsparse_status aoclsparse_ctrsm(aoclsparse_operation trans, const aoclsparse_float_complex alpha, aoclsparse_matrix A, const aoclsparse_mat_descr descr, aoclsparse_order order, const aoclsparse_float_complex *B, aoclsparse_int n, aoclsparse_int ldb, aoclsparse_float_complex *X, aoclsparse_int ldx)#

aoclsparse_status aoclsparse_ztrsm(aoclsparse_operation trans, const aoclsparse_double_complex alpha, aoclsparse_matrix A, const aoclsparse_mat_descr descr, aoclsparse_order order, const aoclsparse_double_complex *B, aoclsparse_int n, aoclsparse_int ldb, aoclsparse_double_complex *X, aoclsparse_int ldx)#

Solve sparse triangular linear system of equations with multiple right hand sides for real/complex single and double data precisions.

aoclsparse_?trsm solves a sparse triangular linear system of equations with multiple right hand sides, of the form

\[ op(A)\; X = \alpha B, \]

where \(A\) is a sparse matrix of size \(m\), \(op()\) is a linear operator, \(X\) and \(B\) are rectangular dense matrices of appropiate size, while \(\alpha\) is a scalar. The sparse matrix \(A\) can be interpreted either as a lower triangular or upper triangular. This is indicated by fill_mode from the matrix descriptor descr where either upper or lower triangular portion of the matrix is only referenced. The matrix can also be of class symmetric in which case only the selected triangular part 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 Hermitian transposition operations. By default, no transposition is performed. The right-hand-side matrix \(B\) and the solution matrix \(X\) are dense and must be of the correct size, that is \(m\) by \(n\), see ldb and ldx input parameters for further details.

Explicitly, this kernel solves

\[ op(A)\; X = \alpha \; B\text{, with solution } X = \alpha \; (op(A)^{-1}) \; B, \]

where

If a linear operator is applied then, the possible problems solved are

\[ A^T \; X = \alpha \; B\text{, with solution } X = \alpha \; A^{-T} \; B\text{, and } \; A^H \; X = \alpha \; B\text{, with solution } X = \alpha \; A^{-H} \; B. \]

Notes

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.
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.
This set of APIs allocates work array of size \(m\) for each case where the matrices \(B\) or \(X\) are stored in row-major format (ref aoclsparse_order_row).
A subset of kernels are parallel (on parallel builds) and can be expected potential acceleration in the solve. These kernels are available when both dense matrices \(X\) and \(B\) are stored in column-major format (ref aoclsparse_order_column) and thread count is greater than 1 on a parallel build.
There is aoclsparse_trsm_kid (Kernel ID) variation of TRSM, namely with a suffix of _kid, this solver allows to choose which underlying TRSV kernels to use (if possible). Currently, all the existing kernels avilable in aoclsparse_trsv_kid are supported.
This routine supports only sparse matrices in CSR format.

Example - Real space (tests/examples/sample_dtrsm.cpp)

/* ************************************************************************
 * Copyright (c) 2023-2024 Advanced Micro Devices, Inc.
 *
 * Permission is hereby granted, free of charge, to any person obtaining a copy
 * of this software and associated documentation files (the "Software"), to deal
 * in the Software without restriction, including without limitation the rights
 * to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
 * copies of the Software, and to permit persons to whom the Software is
 * furnished to do so, subject to the following conditions:
 *
 * The above copyright notice and this permission notice shall be included in
 * all copies or substantial portions of the Software.
 *
 * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
 * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
 * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
 * AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
 * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
 * OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
 * THE SOFTWARE.
 *
 * ************************************************************************ */

#include "aoclsparse.h"

#include <iomanip>
#include <iostream>
#include <limits>
#include <vector>

int main(void)
{
    std::cout << "------------------------------------------------" << std::endl
              << " Triangle Solve with Multiple Right Hand Sides" << std::endl
              << " sample program for real data type" << std::endl
              << "------------------------------------------------" << std::endl
              << std::endl;

    /*
     * This example illustrates how to solve two triangular
     * linear system of equations. The real matrix used is
     *     | 1  2  0  0 |
     * A = | 3  1  2  0 |
     *     | 0  3  1  2 |
     *     | 0  0  3  1 |
     *
     * The first linear system is L X = alpha * B with L = tril(A), X and B
     * are two dense matrices of size 4x2 stored in column-major format.
     *
     * The second linear system is U X = alpha * B with U = triu(A), X and B
     * are two dense matrices of size 4x2 stored in row-major format.
     *
     */

    // Create a tri-diagonal matrix in CSR format
    const aoclsparse_int n = 4, m = 4, nnz = 10;
    aoclsparse_int       icrow[5] = {0, 2, 5, 8, 10};
    aoclsparse_int       icol[10] = {0, 1, 0, 1, 2, 1, 2, 3, 2, 3};
    double               aval[10] = {1, 2, 3, 1, 2, 3, 1, 2, 3, 1};
    aoclsparse_status    status;
    const aoclsparse_int k = 2;
    aoclsparse_int       ldb, ldx;
    bool                 oki, ok;
    double               tol = 20.0 * std::numeric_limits<double>::epsilon();

    // Create aoclsparse matrix and its descriptor
    aoclsparse_matrix     A;
    aoclsparse_index_base base = aoclsparse_index_base_zero;
    aoclsparse_order      order;
    aoclsparse_mat_descr  descr_a;
    aoclsparse_operation  trans;
    status = aoclsparse_create_dcsr(&A, base, m, n, nnz, icrow, icol, aval);
    if(status != aoclsparse_status_success)
    {
        std::cerr << "Error returned from aoclsparse_create_dcsr, status = " << status << "."
                  << std::endl;
        return 3;
    }
    aoclsparse_create_mat_descr(&descr_a);

    /* Case 1: Solving the lower triangular system L^T X = alpha*B.
     *     |1, -1|
     * B = |2, -2| stored in column-major layout
     *     |3, -6| k = 2 and ldb = 4
     *     |4, -8|
     *
     * Linear system L X = B, with solution matrix
     *     | 86,-167|
     * X = |-29,  56| stored in column-major layout
     *     |  9, -18| k = 2 and ldx = 4
     *     | -4,   8|
     */

    double alpha = -1.0;
    double X[m * k];

    double              B[m * k] = {1, 2, 3, 4, -1, -2, -6, -8};
    std::vector<double> XRef;
    XRef.assign({86, -29, 9, -4, -167, 56, -18, 8});

    // Indicate to only access the lower triangular part of matrix A
    aoclsparse_set_mat_type(descr_a, aoclsparse_matrix_type_triangular);
    aoclsparse_set_mat_fill_mode(descr_a, aoclsparse_fill_mode_lower);
    trans = aoclsparse_operation_transpose;

    order  = aoclsparse_order_column;
    ldb    = m;
    ldx    = m;
    status = aoclsparse_set_sm_hint(A, trans, descr_a, order, 1);
    if(status != aoclsparse_status_success)
    {
        std::cerr << "Error returned from aoclsparse_set_sm_hint, status = " << status << "."
                  << std::endl;
        return 3;
    }
    status = aoclsparse_optimize(A);
    if(status != aoclsparse_status_success)
    {
        std::cerr << "Error returned from aoclsparse_optimize, status = " << status << "."
                  << std::endl;
        return 3;
    }
    // Solve
    status = aoclsparse_dtrsm(trans, alpha, A, descr_a, order, &B[0], k, ldb, &X[0], ldx);
    if(status != aoclsparse_status_success)
    {
        std::cerr << "Error returned from aoclsparse_dtrsm, status = " << status << "."
                  << std::endl;
        return 3;
    }

    // Print and check the result
    std::cout << "Solving L^T X = alpha B, where  L=tril(A), and" << std::endl;
    std::cout << "  X and B are dense rectangular martices (column-major layout)" << std::endl;
    std::cout << "  Solution matrix X = " << std::endl;
    std::cout << std::fixed;
    std::cout.precision(1);
    ok = true;
    size_t idx;
    for(int row = 0; row < m; ++row)
    {
        for(int col = 0; col < k; ++col)
        {
            idx = row + col * ldx;
            oki = std::abs(X[idx] - XRef[idx]) <= tol;
            std::cout << X[idx] << (oki ? "  " : "! ");
            ok &= oki;
        }
        std::cout << std::endl;
    }
    std::cout << std::endl;

    /* Case 2: Solving the lower triangular system U X = alpha*B.
     * In this case we "reinterpret" B as a col-major layout, having the
     * form
     *     | 1,  2|
     * B = | 3,  4| stored in row-major layout
     *     |-1, -2| k = 2 and ldb = 2
     *     |-6, -8|
     *
     * Linear system L X = B, with solution matrix
     *     | 39,  50|
     * X = |-19, -24| stored in row-major layout
     *     | 11,  14| k = 2 and ldx = 2
     *     | -6,  -8|
     */

    // Store same XRef in row-major format
    XRef.assign({39, 50, -19, -24, 11, 14, -6, -8});
    alpha = 1.0;
    ldb   = k;
    ldx   = k;
    // Indicate to use only the conjugate transpose of the upper part of A
    aoclsparse_set_mat_fill_mode(descr_a, aoclsparse_fill_mode_upper);
    trans  = aoclsparse_operation_none;
    order  = aoclsparse_order_row;
    status = aoclsparse_set_sm_hint(A, trans, descr_a, order, 1);
    if(status != aoclsparse_status_success)
    {
        std::cerr << "Error returned from aoclsparse_set_sm_hint, status = " << status << "."
                  << std::endl;
        return 3;
    }
    status = aoclsparse_optimize(A);
    if(status != aoclsparse_status_success)
    {
        std::cerr << "Error returned from aoclsparse_optimize, status = " << status << "."
                  << std::endl;
        return 3;
    }
    // Solve
    status = aoclsparse_dtrsm(trans, alpha, A, descr_a, order, &B[0], k, ldb, &X[0], ldx);
    if(status != aoclsparse_status_success)
    {
        std::cerr << "Error returned from aoclsparse_dtrsm, status = " << status << "."
                  << std::endl;
        return 3;
    }
    // Print and check the result
    std::cout << "Solving U X = alpha B, where  U=triu(A), and" << std::endl;
    std::cout << "  X and B are dense rectangular martices (row-major layout)" << std::endl;
    std::cout << "  Solution matrix X = " << std::endl;
    for(int row = 0; row < m; ++row)
    {
        for(int col = 0; col < k; col++)
        {
            idx = row * ldx + col;
            oki = std::abs(X[idx] - XRef[idx]) <= tol;
            std::cout << X[idx] << (oki ? " " : "!  ");
            ok &= oki;
        }
        std::cout << std::endl;
    }
    std::cout << std::endl;

    // Destroy the aoclsparse memory
    aoclsparse_destroy_mat_descr(descr_a);
    aoclsparse_destroy(&A);

    return ok ? 0 : 5;
}

Example - Complex space (tests/examples/sample_ztrsm.cpp)

/* ************************************************************************
 * Copyright (c) 2023-2024 Advanced Micro Devices, Inc.
 *
 * Permission is hereby granted, free of charge, to any person obtaining a copy
 * of this software and associated documentation files (the "Software"), to deal
 * in the Software without restriction, including without limitation the rights
 * to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
 * copies of the Software, and to permit persons to whom the Software is
 * furnished to do so, subject to the following conditions:
 *
 * The above copyright notice and this permission notice shall be included in
 * all copies or substantial portions of the Software.
 *
 * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
 * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
 * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
 * AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
 * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
 * OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
 * THE SOFTWARE.
 *
 * ************************************************************************ */

#include "aoclsparse.h"

#include <complex>
#include <iomanip>
#include <iostream>
#include <limits>
#include <vector>

int main(void)
{
    std::cout << "------------------------------------------------" << std::endl
              << " Triangle Solve with Multiple Right Hand Sides" << std::endl
              << " sample program for complex data type" << std::endl
              << "------------------------------------------------" << std::endl
              << std::endl;

    /*
     * This example illustrates how to solve two triangular
     * linear system of equations. The complex matrix used is
     *     | 1+3i  2+5i     0     0 |
     * A = | 3     1+2i  2-2i     0 |
     *     |    0     3     1  2+3i |
     *     |    0     0  3+1i  1+1i |
     *
     * The first linear system is L X = alpha * B with L = tril(A), X and B
     * are two dense matrices of size 4x2 stored in column-major format.
     *
     * The second linear system is U^H X = alpha * B with U = triu(A), X and B
     * are two dense matrices of size 4x2 stored in row-major format.
     *
     */

    // Create a tri-diagonal matrix in CSR format
    const aoclsparse_int                   n = 4, m = 4, nnz = 10;
    aoclsparse_int                         icrow[5] = {0, 2, 5, 8, 10};
    aoclsparse_int                         icol[10] = {0, 1, 0, 1, 2, 1, 2, 3, 2, 3};
    std::vector<aoclsparse_double_complex> aval;
    // clang-format off
    aval.assign({{1, 3}, {2, 5}, {3, 0}, {1, 2}, {2, -2}, {3, 0},
                 {1, 0}, {2, 3}, {3, 1}, {1, 1}});
    // clang-format on
    aoclsparse_status status;
    aoclsparse_int    ldb, ldx;
    bool              oki, ok;
    double            tol = 20 * std::numeric_limits<double>::epsilon();

    aoclsparse_matrix     A;
    const aoclsparse_int  k     = 2; // number of columns of X and B
    aoclsparse_index_base base  = aoclsparse_index_base_zero;
    aoclsparse_order      order = aoclsparse_order_column;
    aoclsparse_mat_descr  descr_a;
    aoclsparse_operation  trans = aoclsparse_operation_none;
    status = aoclsparse_create_zcsr(&A, base, m, n, nnz, icrow, icol, aval.data());
    if(status != aoclsparse_status_success)
    {
        std::cerr << "Error returned from aoclsparse_create_zcsr, status = " << status << "."
                  << std::endl;
        return 3;
    }
    aoclsparse_create_mat_descr(&descr_a);

    /* Case 1: Solving the lower triangular system L X = B.
     *     |-2+4i,   -3+i|
     * B = |3+13i,  4+11i| stored in column-major layout
     *     |16+7i,  16+2i| k = 2 and ldb = 4
     *     |15+11i,10+16i|
     *
     * Linear system L X = B, with solution matrix
     *     |1+i,     i|
     * X = |4+2i,    4| stored in column-major layout
     *     |4+i,  4+2i| k = 2 and ldx = 4
     *     |4,    3+3i|
     */

    aoclsparse_double_complex alpha = {1.0, 0.};
    aoclsparse_double_complex X[m * k];

    std::vector<aoclsparse_double_complex> B;
    B.assign({{-2, 4}, {3, 13}, {16, 7}, {15, 11}, {-3, 1}, {4, 11}, {16, 2}, {10, 16}});
    std::vector<aoclsparse_double_complex> XRef;
    XRef.assign({{1, 1}, {4, 2}, {4, 1}, {4, 0}, {0, 1}, {4, 0}, {4, 2}, {3, 3}});

    // indicate to only access the lower triangular part of matrix A
    aoclsparse_set_mat_type(descr_a, aoclsparse_matrix_type_triangular);
    aoclsparse_set_mat_fill_mode(descr_a, aoclsparse_fill_mode_lower);

    // Prepare and call solver
    order  = aoclsparse_order_column;
    ldb    = m;
    ldx    = m;
    status = aoclsparse_set_sm_hint(A, trans, descr_a, order, 1);
    if(status != aoclsparse_status_success)
    {
        std::cerr << "Error returned from aoclsparse_set_sm_hint, status = " << status << "."
                  << std::endl;
        return 3;
    }
    status = aoclsparse_optimize(A);
    if(status != aoclsparse_status_success)
    {
        std::cerr << "Error returned from aoclsparse_optimize, status = " << status << "."
                  << std::endl;
        return 3;
    }

    // Solve
    status = aoclsparse_ztrsm(trans, alpha, A, descr_a, order, B.data(), k, ldb, &X[0], ldx);
    if(status != aoclsparse_status_success)
    {
        std::cerr << "Error from aoclsparse_ztrsm, status = " << status << std::endl;
        return 3;
    }

    // Print and check the result
    std::cout << "Solving L X = alpha B, where  L=tril(A), and" << std::endl;
    std::cout << "  X and B are dense rectangular martices (column-major layout)" << std::endl;
    std::cout << "  Solution matrix X = " << std::endl;
    std::cout << std::fixed;
    std::cout.precision(1);
    ok = true;
    size_t idx;
    for(int row = 0; row < m; ++row)
    {
        for(int col = 0; col < k; ++col)
        {
            idx = row + col * ldx;
            oki = std::abs(X[idx].real - XRef[idx].real) <= tol
                  && std::abs(X[idx].imag - XRef[idx].imag) <= tol;
            std::cout << "(" << X[idx].real << ", " << X[idx].imag << "i) " << (oki ? "  " : "! ");
            ok &= oki;
        }
        std::cout << std::endl;
    }
    std::cout << std::endl;

    /* Case 2: Solving the lower triangular system U^H X = alpha*B.
     *     |2+4i,    -1+3i|
     * B = |9+15i,    6+9i| stored in row-major layout
     *     |-13+8i,-10+12i| k = 2 and ldx = 2
     *     |14+15i,  8+20i|
     *
     * Linear system L X = B, with solution matrix
     *     |1+i,     i|
     * X = |4+2i,    4| stored in row-major layout
     *     |4+i,  4+2i| k = 2 and ldx = 2
     *     |4,    3+3i|
     */

    B.assign({{2, 4}, {-1, 3}, {9, 15}, {6, 9}, {-13, 8}, {-10, 12}, {14, 15}, {8, 20}});
    // Store same XRef in row-major format
    XRef.assign({{1, 1}, {0, 1}, {4, 2}, {4, 0}, {4, 1}, {4, 2}, {4, 0}, {3, 3}});
    alpha = {0, -1};
    ldb   = k;
    ldx   = k;
    // Indicate to use only the conjugate transpose of the upper part of A
    trans = aoclsparse_operation_conjugate_transpose;
    order = aoclsparse_order_row;
    aoclsparse_set_mat_fill_mode(descr_a, aoclsparse_fill_mode_upper);
    status = aoclsparse_set_sm_hint(A, trans, descr_a, order, 1);
    if(status != aoclsparse_status_success)
    {
        std::cerr << "Error returned from aoclsparse_set_sm_hint, status = " << status << "."
                  << std::endl;
        return 3;
    }
    status = aoclsparse_optimize(A);
    if(status != aoclsparse_status_success)
    {
        std::cerr << "Error returned from aoclsparse_optimize, status = " << status << "."
                  << std::endl;
        return 3;
    }
    // Solve
    status = aoclsparse_ztrsm(trans, alpha, A, descr_a, order, B.data(), k, ldb, &X[0], ldx);
    if(status != aoclsparse_status_success)
    {
        std::cerr << "Error from aoclsparse_ztrsm, status = " << status << std::endl;
        return 3;
    }
    // Print and check the result
    std::cout << "Solving U^H X = alpha B, where  U=triu(A), and" << std::endl;
    std::cout << "  X and B are dense rectangular martices (row-major layout)" << std::endl;
    std::cout << "  Solution matrix X = " << std::endl;
    for(int row = 0; row < m; ++row)
    {
        for(int col = 0; col < k; col++)
        {
            idx = row * ldx + col;
            oki = std::abs(X[idx].real - XRef[idx].real) <= tol
                  && std::abs(X[idx].imag - XRef[idx].imag) <= tol;
            std::cout << "(" << X[idx].real << ", " << X[idx].imag << "i) " << (oki ? "  " : "! ");
            ok &= oki;
        }
        std::cout << std::endl;
    }
    std::cout << std::endl;
    // Destroy the aoclsparse memory
    aoclsparse_destroy_mat_descr(descr_a);
    aoclsparse_destroy(&A);

    return ok ? 0 : 5;
}

Parameters

trans – [in] matrix operation to perform on \(A\). Possible values are aoclsparse_operation_none, aoclsparse_operation_transpose, and aoclsparse_operation_conjugate_transpose.
alpha – [in] scalar \(\alpha\).
A – [in] sparse matrix \(A\) of size \(m\).
descr – [in] descriptor of the sparse matrix \(A\).
order – [in] storage order of dense matrices \(B\) and \(X\). Possible options are aoclsparse_order_row and aoclsparse_order_column.
B – [in] dense matrix, potentially rectangular, of size \(m \times n\).
n – [in] \(n,\) number of columns of the dense matrix \(B\).
ldb – [in] leading dimension of \(B\). Eventhough the matrix \(B\) is considered of size \(m \times n\), its memory layout may correspond to a larger matrix (ldb by \(N>n\)) in which only the submatrix \(B\) is of interest. In this case, this parameter provides means to access the correct elements of \(B\) within the larger layout.

matrix layout

row count

column count

aoclsparse_order_row

\(m\)

ldb with ldb \(\ge n\)

aoclsparse_order_column

ldb with ldb \(\ge m\)

\(n\)
X – [out] solution matrix \(X,\) dense and potentially rectangular matrix of size \(m \times n\).
ldx – [in] leading dimension of \(X\). Eventhough the matrix \(X\) is considered of size \(m \times n\), its memory layout may correspond to a larger matrix (ldx by \(N>n\)) in which only the submatrix \(X\) is of interest. In this case, this parameter provides means to access the correct elements of \(X\) within the larger layout.

matrix layout

row count

column count

aoclsparse_order_row

\(m\)

ldx with ldx \(\ge n\)

aoclsparse_order_column

ldx with ldx \(\ge m\)

\(n\)

Return values

matrix layout	row count	column count
aoclsparse_order_row	\(m\)	`ldb` with `ldb` \(\ge n\)
aoclsparse_order_column	`ldb` with `ldb` \(\ge m\)	\(n\)

matrix layout	row count	column count
aoclsparse_order_row	\(m\)	`ldx` with `ldx` \(\ge n\)
aoclsparse_order_column	`ldx` with `ldx` \(\ge m\)	\(n\)

aoclsparse_status_success – indicates that the operation completed successfully.
aoclsparse_status_invalid_size – informs that either m, n, nnz, ldb or ldx is invalid.
aoclsparse_status_invalid_pointer – informs that either descr, alpha, A, B, or X pointer is invalid.
aoclsparse_status_not_implemented – this error occurs when the provided matrix aoclsparse_matrix_type is aoclsparse_matrix_type_general or aoclsparse_matrix_type_hermitian or when matrix A is not in CSR format.

aoclsparse_status aoclsparse_strsm_kid(const aoclsparse_operation trans, const float alpha, aoclsparse_matrix A, const aoclsparse_mat_descr descr, aoclsparse_order order, const float *B, aoclsparse_int n, aoclsparse_int ldb, float *X, aoclsparse_int ldx, const aoclsparse_int kid)#

aoclsparse_status aoclsparse_dtrsm_kid(const aoclsparse_operation trans, const double alpha, aoclsparse_matrix A, const aoclsparse_mat_descr descr, aoclsparse_order order, const double *B, aoclsparse_int n, aoclsparse_int ldb, double *X, aoclsparse_int ldx, const aoclsparse_int kid)#

aoclsparse_status aoclsparse_ctrsm_kid(aoclsparse_operation trans, const aoclsparse_float_complex alpha, aoclsparse_matrix A, const aoclsparse_mat_descr descr, aoclsparse_order order, const aoclsparse_float_complex *B, aoclsparse_int n, aoclsparse_int ldb, aoclsparse_float_complex *X, aoclsparse_int ldx, const aoclsparse_int kid)#

aoclsparse_status aoclsparse_ztrsm_kid(aoclsparse_operation trans, const aoclsparse_double_complex alpha, aoclsparse_matrix A, const aoclsparse_mat_descr descr, aoclsparse_order order, const aoclsparse_double_complex *B, aoclsparse_int n, aoclsparse_int ldb, aoclsparse_double_complex *X, aoclsparse_int ldx, const aoclsparse_int kid)#

Solve sparse triangular linear system of equations with multiple right hand sides for real/complex single and double data precisions (kernel flag variation).

For full details refer to aoclsparse_?trsm().

This variation of TRSM, namely with a suffix of _kid, allows to choose which underlying TRSV kernels to use (if possible). Currently, all the existing kernels supported by aoclsparse_?trsv_kid() are available here as well.

Parameters

trans – [in] matrix operation to perform on \(A\). Possible values are aoclsparse_operation_none, aoclsparse_operation_transpose, and aoclsparse_operation_conjugate_transpose.
alpha – [in] scalar \(\alpha\).
A – [in] sparse matrix \(A\) of size \(m\).
descr – [in] descriptor of the sparse matrix \(A\).
order – [in] storage order of dense matrices \(B\) and \(X\). Possible options are aoclsparse_order_row and aoclsparse_order_column.
B – [in] dense matrix, potentially rectangular, of size \(m \times n\).
n – [in] \(n,\) number of columns of the dense matrix \(B\).
ldb – [in] leading dimension of \(B\). Eventhough the matrix \(B\) is considered of size \(m \times n\), its memory layout may correspond to a larger matrix (ldb by \(N>n\)) in which only the submatrix \(B\) is of interest. In this case, this parameter provides means to access the correct elements of \(B\) within the larger layout.

matrix layout

row count

column count

aoclsparse_order_row

\(m\)

ldb with ldb \(\ge n\)

aoclsparse_order_column

ldb with ldb \(\ge m\)

\(n\)
X – [out] solution matrix \(X,\) dense and potentially rectangular matrix of size \(m \times n\).
ldx – [in] leading dimension of \(X\). Eventhough the matrix \(X\) is considered of size \(m \times n\), its memory layout may correspond to a larger matrix (ldx by \(N>n\)) in which only the submatrix \(X\) is of interest. In this case, this parameter provides means to access the correct elements of \(X\) within the larger layout.

matrix layout

row count

column count

aoclsparse_order_row

\(m\)

ldx with ldx \(\ge n\)

aoclsparse_order_column

ldx with ldx \(\ge m\)

\(n\)
kid – [in] kernel ID, hints which kernel to use.

matrix layout	row count	column count
aoclsparse_order_row	\(m\)	`ldb` with `ldb` \(\ge n\)
aoclsparse_order_column	`ldb` with `ldb` \(\ge m\)	\(n\)

matrix layout	row count	column count
aoclsparse_order_row	\(m\)	`ldx` with `ldx` \(\ge n\)
aoclsparse_order_column	`ldx` with `ldx` \(\ge m\)	\(n\)

aoclsparse_sp2m()#

aoclsparse_status aoclsparse_sp2m(aoclsparse_operation opA, const aoclsparse_mat_descr descrA, const aoclsparse_matrix A, aoclsparse_operation opB, const aoclsparse_mat_descr descrB, const aoclsparse_matrix B, const aoclsparse_request request, aoclsparse_matrix *C)#

Sparse matrix Sparse matrix multiplication for real and complex datatypes.

aoclsparse_?sp2m multiplies two sparse matrices in CSR storage format. The result is stored in a newly allocated sparse matrix in CSR format, such that

\[ C = op(A) \, op(B), \]

with

\[\begin{split} op(A) = \left\{ \begin{array}{ll} A, & \text{if } {\bf\mathsf{opA}} = \text{aoclsparse}\_\text{operation}\_\text{none} \\ A^T, & \text{if } {\bf\mathsf{opA}} = \text{aoclsparse}\_\text{operation}\_\text{transpose} \\ A^H, & \text{if } {\bf\mathsf{opA}} = \text{aoclsparse}\_\text{operation}\_\text{conjugate}\_\text{transpose} \end{array} \right. \end{split}\]

and

\[\begin{split} op(B) = \left\{ \begin{array}{ll} B, & \text{if } {\bf\mathsf{opB}} = \text{aoclsparse}\_\text{operation}\_\text{none} \\ B^T, & \text{if } {\bf\mathsf{opB}} = \text{aoclsparse}\_\text{operation}\_\text{transpose} \\ B^H, & \text{if } {\bf\mathsf{opB}} = \text{aoclsparse}\_\text{operation}\_\text{conjugate}\_\text{transpose} \end{array} \right. \end{split}\]

where \(A\) is a \(m \times k\) matrix , \(B\) is a \(k \times n\) matrix, resulting in \(m \times n\) matrix \(C\), for opA and opB = aoclsparse_operation_none. \(A\) is a \(k \times m\) matrix when opA = aoclsparse_operation_transpose or aoclsparse_operation_conjugate_transpose and \(B\) is a \(n \times k\) matrix when opB = aoclsparse_operation_transpose or aoclsparse_operation_conjugate_transpose

aoclsparse_sp2m can be run in single-stage or two-stage. The single-stage algorithm allocates and computes the entire output matrix in a single stage aoclsparse_stage_full_computation. Whereas, in two-stage algorithm, the first stage aoclsparse_stage_nnz_count allocates memory for the output matrix and computes the number of entries of the matrix. The second stage aoclsparse_stage_finalize computes column indices of non-zero elements and values of the output matrix. The second stage has to be invoked only after the first stage. But, it can be also be invoked multiple times consecutively when the sparsity structure of input matrices remains unchanged, with only the values getting updated.

Example - Complex space (tests/examples/sample_zsp2m.cpp)

/* ************************************************************************
 * Copyright (c) 2023 Advanced Micro Devices, Inc. All rights reserved.
 *
 * Permission is hereby granted, free of charge, to any person obtaining a copy
 * of this software and associated documentation files (the "Software"), to deal
 * in the Software without restriction, including without limitation the rights
 * to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
 * copies of the Software, and to permit persons to whom the Software is
 * furnished to do so, subject to the following conditions:
 *
 * The above copyright notice and this permission notice shall be included in
 * all copies or substantial portions of the Software.
 *
 * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
 * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
 * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
 * AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
 * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
 * OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
 * THE SOFTWARE.
 *
 * ************************************************************************ */

#include "aoclsparse.h"

#include <complex>
#include <iomanip>
#include <iostream>

/* Computes multiplication of 2 sparse matrices A and B in CSR format. */

// Comment this out for single stage computation
//#define TWO_STAGE_COMPUTATION

int main(void)
{
    std::cout << "-------------------------------" << std::endl
              << "----- SP2M sample program -----" << std::endl
              << "-------------------------------" << std::endl
              << std::endl;

    aoclsparse_status     status;
    aoclsparse_int        nnz_C;
    aoclsparse_request    request;
    aoclsparse_index_base base = aoclsparse_index_base_zero;
    aoclsparse_operation  opA  = aoclsparse_operation_transpose;
    aoclsparse_operation  opB  = aoclsparse_operation_none;

    // Print aoclsparse version
    std::cout << aoclsparse_get_version() << std::endl;

    // Initialise matrix descriptor and csr matrix structure of inputs A and B
    aoclsparse_mat_descr descrA;
    aoclsparse_mat_descr descrB;
    aoclsparse_matrix    csrA;
    aoclsparse_matrix    csrB;

    // Create matrix descriptor of input matrices
    // aoclsparse_create_mat_descr set aoclsparse_matrix_type to aoclsparse_matrix_type_general
    // and aoclsparse_index_base to aoclsparse_index_base_zero.
    aoclsparse_create_mat_descr(&descrA);
    aoclsparse_create_mat_descr(&descrB);

    // Matrix sizes
    aoclsparse_int m = 5, n = 5, k = 5;
    aoclsparse_int nnz_A = 10, nnz_B = 10;
    // Matrix A
    aoclsparse_int            row_ptr_A[] = {0, 1, 2, 5, 9, 10};
    aoclsparse_int            col_ind_A[] = {0, 0, 1, 2, 4, 0, 1, 2, 3, 4};
    aoclsparse_double_complex val_A[]     = {{-0.86238, 0.454626},
                                             {-2.62138, -0.442597},
                                             {-0.875679, 0.137933},
                                             {-0.661939, -1.09106},
                                             {0.0501717, -2.37527},
                                             {-1.48812, -0.420546},
                                             {-0.588085, -0.708977},
                                             {0.310933, -0.96569},
                                             {-0.88964, -2.37881},
                                             {-1.23201, 0.213152}};
    aoclsparse_create_zcsr(
        &csrA, base, k, m, nnz_A, row_ptr_A, col_ind_A, (aoclsparse_double_complex *)val_A);

    // Matrix B
    aoclsparse_int            row_ptr_B[] = {0, 4, 4, 7, 8, 10};
    aoclsparse_int            col_ind_B[] = {0, 1, 2, 4, 0, 1, 2, 2, 2, 3};
    aoclsparse_double_complex val_B[]     = {{-1.59204, -0.259325},
                                             {0.467532, -0.980612},
                                             {0.078412, -0.513591},
                                             {-1.52364, 0.403911},
                                             {0.211966, -1.33485},
                                             {-1.37901, -1.44562},
                                             {1.42472, -2.08662},
                                             {-2.26549, -1.0073},
                                             {-1.75098, 0.207783},
                                             {-1.8152, 0.482205}};

    aoclsparse_create_zcsr(
        &csrB, base, k, n, nnz_B, row_ptr_B, col_ind_B, (aoclsparse_double_complex *)val_B);

    aoclsparse_matrix          csrC          = NULL;
    aoclsparse_int            *csr_row_ptr_C = NULL;
    aoclsparse_int            *csr_col_ind_C = NULL;
    aoclsparse_double_complex *csr_val_C     = NULL;
    aoclsparse_int             C_M, C_N;

    // expected output matrices
    aoclsparse_int            C_M_exp = 5, C_N_exp = 5, nnz_C_exp = 15;
    aoclsparse_int            csr_row_ptr_C_exp[] = {0, 4, 7, 10, 11, 15};
    aoclsparse_int            csr_col_ind_C_exp[] = {0, 1, 2, 4, 0, 1, 2, 0, 1, 2, 2, 0, 1, 2, 3};
    aoclsparse_double_complex csr_val_C_exp[]     = {{1.49084, -0.500145},
                                                     {0.0426217, 1.05821},
                                                     {3.11358, 2.93029},
                                                     {1.13033, -1.04101},
                                                     {-0.0014949, 1.19813},
                                                     {1.40697, 1.07569},
                                                     {-0.34163, 4.22229},
                                                     {-1.59671, 0.652319},
                                                     {-0.664445, 2.4615},
                                                     {-4.89686, 1.70132},
                                                     {-0.380707, 6.28532},
                                                     {-3.15998, -0.570448},
                                                     {-3.50293, 3.20298},
                                                     {-2.77187, -4.11799},
                                                     {2.13356, -0.980994}};

#ifdef TWO_STAGE_COMPUTATION
    std::cout << "Invoking aoclsparse_sp2m with aoclsparse_stage_nnz_count..\n";
    // aoclsparse_stage_nnz_count : Only rowIndex array of the CSR matrix
    // is computed internally.
    request = aoclsparse_stage_nnz_count;
    status  = aoclsparse_sp2m(opA, descrA, csrA, opB, descrB, csrB, request, &csrC);
    if(status != aoclsparse_status_success)
        return 1;

    std::cout << "Invoking aoclsparse_sp2m with aoclsparse_stage_finalize..\n";
    // aoclsparse_stage_finalize : Finalize computation of remaining
    // output arrays ( column indices and values of output matrix entries) .
    // Has to be called only after aoclsparse_sp2m call with
    // aoclsparse_stage_nnz_count parameter.
    request = aoclsparse_stage_finalize;
    status  = aoclsparse_sp2m(opA, descrA, csrA, opB, descrB, csrB, request, &csrC);
    if(status != aoclsparse_status_success)
        return 2;

#else // SINGLE STAGE
    std::cout << "Invoking aoclsparse_sp2m with aoclsparse_stage_full_computation..\n";
    // aoclsparse_stage_full_computation :  Whole computation is performed in
    // single step.
    request = aoclsparse_stage_full_computation;
    status  = aoclsparse_sp2m(opA, descrA, csrA, opB, descrB, csrB, request, &csrC);
    if(status != aoclsparse_status_success)
        return 3;

#endif

    aoclsparse_export_zcsr(
        csrC, &base, &C_M, &C_N, &nnz_C, &csr_row_ptr_C, &csr_col_ind_C, &csr_val_C);
    // Check and print the result
    std::cout << std::fixed;
    std::cout.precision(1);
    bool oka, okb, okc, oki, okj, okk, ok = true;
    std::cout << std::endl
              << "Output Matrix C: " << std::endl
              << std::setw(11) << "C_M" << std::setw(3) << "" << std::setw(11) << "expected"
              << std::setw(3) << "" << std::setw(11) << "C_N" << std::setw(3) << "" << std::setw(11)
              << "expected" << std::setw(3) << "" << std::setw(11) << "nnz_C" << std::setw(3) << ""
              << std::setw(11) << "expected" << std::endl;
    oka = C_M == C_M_exp;
    ok &= oka;
    std::cout << std::setw(11) << C_M << std::setw(3) << "" << std::setw(11) << C_M_exp
              << std::setw(2) << (oka ? "" : " !");
    okb = C_N == C_N_exp;
    ok &= okb;
    std::cout << std::setw(11) << C_N << std::setw(3) << "" << std::setw(11) << C_N_exp
              << std::setw(2) << (okb ? "" : " !");
    okc = nnz_C == nnz_C_exp;
    ok &= okc;
    std::cout << std::setw(11) << nnz_C << std::setw(3) << "" << std::setw(11) << nnz_C_exp
              << std::setw(2) << (okc ? "" : " !");
    std::cout << std::endl;
    std::cout << std::endl;
    std::cout << std::setw(11) << "csr_val_C" << std::setw(3) << "" << std::setw(11) << "expected"
              << std::setw(3) << "" << std::setw(11) << "csr_col_ind_C" << std::setw(3) << ""
              << std::setw(11) << "expected" << std::setw(3) << "" << std::setw(11)
              << "csr_row_ptr_C" << std::setw(3) << "" << std::setw(11) << "expected" << std::endl;
    //Initializing precision tolerance range for double
    const double tol = 1e-03;
    for(aoclsparse_int i = 0; i < nnz_C; i++)
    {
        oki = ((std::abs(csr_val_C[i].real - csr_val_C_exp[i].real) <= tol)
               && (std::abs(csr_val_C[i].imag - csr_val_C_exp[i].imag) <= tol));
        ok &= oki;
        std::cout << std::setw(11) << "(" << csr_val_C[i].real << ", " << csr_val_C[i].imag << "i) "
                  << std::setw(3) << "" << std::setw(11) << "(" << csr_val_C_exp[i].real << ", "
                  << csr_val_C_exp[i].imag << "i) " << std::setw(2) << (oki ? "" : " !");
        okj = csr_col_ind_C[i] == csr_col_ind_C_exp[i];
        ok &= okj;
        std::cout << std::setw(11) << csr_col_ind_C[i] << std::setw(3) << "" << std::setw(11)
                  << csr_col_ind_C_exp[i] << std::setw(2) << (okj ? "" : " !");
        if(i < C_M)
        {
            okk = csr_row_ptr_C[i] == csr_row_ptr_C_exp[i];
            ok &= okk;
            std::cout << " " << std::setw(11) << csr_row_ptr_C[i] << std::setw(3) << ""
                      << std::setw(11) << csr_row_ptr_C_exp[i] << std::setw(2) << (okk ? "" : " !");
        }
        std::cout << std::endl;
    }

    aoclsparse_destroy_mat_descr(descrA);
    aoclsparse_destroy_mat_descr(descrB);
    aoclsparse_destroy(&csrA);
    aoclsparse_destroy(&csrB);
    aoclsparse_destroy(&csrC);
    return (ok ? 0 : 6);
}

Parameters

opA – [in] matrix \(A\) operation type.
descrA – [in] descriptor of the sparse CSR matrix \(A\). Currently, only aoclsparse_matrix_type_general is supported.
A – [in] sparse CSR matrix \(A\) .
opB – [in] matrix \(B\) operation type.
descrB – [in] descriptor of the sparse CSR matrix \(B\). Currently, only aoclsparse_matrix_type_general is supported.
B – [in] sparse CSR matrix \(B\) .
request – [in] Specifies full computation or two-stage algorithm aoclsparse_stage_nnz_count , Only rowIndex array of the CSR matrix is computed internally. The output sparse CSR matrix can be extracted to measure the memory required for full operation. aoclsparse_stage_finalize . Finalize computation of remaining output arrays ( column indices and values of output matrix entries) . Has to be called only after aoclsparse_sp2m call with aoclsparse_stage_nnz_count parameter. aoclsparse_stage_full_computation . Perform the entire computation in a single step.
*C – [out] Pointer to sparse CSR matrix \(C\) . Matrix \(C\) arrays will always have zero-based indexing, irrespective of matrix \(A\) or matrix \(B\) being one-based or zero-based indexing. The column indices of the output matrix in CSR format can appear unsorted.

Return values

aoclsparse_status_success – the operation completed successfully.
aoclsparse_status_invalid_pointer – descrA, descrB, A, B, C is invalid.
aoclsparse_status_invalid_size – input size parameters contain an invalid value.
aoclsparse_status_invalid_value – input parameters contain an invalid value.
aoclsparse_status_wrong_type – A and B matrix datatypes dont match.
aoclsparse_status_memory_error – Memory allocation failure.
aoclsparse_status_not_implemented – aoclsparse_matrix_type is not aoclsparse_matrix_type_general or input matrices A or B is not in CSR format

aoclsparse_spmm()#

aoclsparse_status aoclsparse_spmm(aoclsparse_operation opA, const aoclsparse_matrix A, const aoclsparse_matrix B, aoclsparse_matrix *C)#

Sparse matrix Sparse matrix multiplication for real and complex datatypes.

aoclsparse_?spmm multiplies two sparse matrices in CSR storage format. The result is stored in a newly allocated sparse matrix in CSR format, such that

\[\begin{split}C = op(A) \cdot B, \text{ with } op(A) = \left\{ \begin{array}{ll} A, & \text{ if } {\bf\mathsf{opA}} = \text{aoclsparse}\_\text{operation}\_\text{none} \\ A^T, & \text{ if } {\bf\mathsf{opA}} = \text{aoclsparse}\_\text{operation}\_\text{transpose} \\ A^H, & \text{ if } {\bf\mathsf{opA}} = \text{aoclsparse}\_\text{operation}\_\text{conjugate}\_\text{transpose} \end{array} \right.\end{split}\]

where \(A\) is a \(m \times k\) matrix , \(B\) is a \(k \times n\) matrix, resulting in \(m \times n\) matrix \(C\), for opA = aoclsparse_operation_none. \(A\) is a \(k \times m\) matrix when opA = aoclsparse_operation_transpose or aoclsparse_operation_conjugate_transpose

Parameters

opA – [in] matrix \(A\) operation type.
A – [in] sparse CSR matrix \(A\).
B – [in] sparse CSR matrix \(B\).
*C – [out] Pointer to sparse CSR matrix \(C\) . Matrix \(C\) arrays will always have zero-based indexing, irrespective of matrix \(A\) or matrix \(B\) being one-based or zero-based indexing. The column indices of the output matrix in CSR format can appear unsorted.

Return values

aoclsparse_status_success – the operation completed successfully.
aoclsparse_status_invalid_pointer – A, B, C is invalid.
aoclsparse_status_invalid_size – input size parameters contain an invalid value.
aoclsparse_status_invalid_value – input parameters contain an invalid value.
aoclsparse_status_wrong_type – A and B matrix data types do not match.
aoclsparse_status_memory_error – Memory allocation failure.
aoclsparse_status_not_implemented – Input matrices A or B is not in CSR format

aoclsparse_?csrmm()#

aoclsparse_status aoclsparse_scsrmm(aoclsparse_operation op, const float alpha, const aoclsparse_matrix A, const aoclsparse_mat_descr descr, aoclsparse_order order, const float *B, aoclsparse_int n, aoclsparse_int ldb, const float beta, float *C, aoclsparse_int ldc)#

aoclsparse_status aoclsparse_dcsrmm(aoclsparse_operation op, const double alpha, const aoclsparse_matrix A, const aoclsparse_mat_descr descr, aoclsparse_order order, const double *B, aoclsparse_int n, aoclsparse_int ldb, const double beta, double *C, aoclsparse_int ldc)#

aoclsparse_status aoclsparse_ccsrmm(aoclsparse_operation op, const aoclsparse_float_complex alpha, const aoclsparse_matrix A, const aoclsparse_mat_descr descr, aoclsparse_order order, const aoclsparse_float_complex *B, aoclsparse_int n, aoclsparse_int ldb, const aoclsparse_float_complex beta, aoclsparse_float_complex *C, aoclsparse_int ldc)#

aoclsparse_status aoclsparse_zcsrmm(aoclsparse_operation op, const aoclsparse_double_complex alpha, const aoclsparse_matrix A, const aoclsparse_mat_descr descr, aoclsparse_order order, const aoclsparse_double_complex *B, aoclsparse_int n, aoclsparse_int ldb, const aoclsparse_double_complex beta, aoclsparse_double_complex *C, aoclsparse_int ldc)#

Sparse matrix dense matrix multiplication using CSR storage format.

aoclsparse_?csrmm multiplies a scalar \(\alpha\) with a sparse \(m \times k\) matrix \(A\), defined in CSR storage format, and a dense \(k \times n\) matrix \(B\) and adds the result to the dense \(m \times n\) matrix \(C\) that is multiplied by a scalar \(\beta\), such that

\[\begin{split} C = \alpha \, op(A) \, B + \beta \, C, \quad \text{ with } \quad op(A) = \left\{ \begin{array}{ll} A, & \text{ if } {\bf\mathsf{trans}\_\mathsf{A}} = \text{aoclsparse}\_\text{operation}\_\text{none} \\ A^T, & \text{ if } {\bf\mathsf{trans}\_\mathsf{A}} = \text{aoclsparse}\_\text{operation}\_\text{transpose} \\ A^H, & \text{ if } {\bf\mathsf{trans}\_\mathsf{A}} = \text{aoclsparse}\_\text{operation}\_\text{conjugate}\_\text{transpose} \end{array} \right. \end{split}\]

Example (tests/examples/sample_csrmm.cpp)

/* ************************************************************************
 * Copyright (c) 2023 Advanced Micro Devices, Inc.
 *
 * Permission is hereby granted, free of charge, to any person obtaining a copy
 * of this software and associated documentation files (the "Software"), to deal
 * in the Software without restriction, including without limitation the rights
 * to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
 * copies of the Software, and to permit persons to whom the Software is
 * furnished to do so, subject to the following conditions:
 *
 * The above copyright notice and this permission notice shall be included in
 * all copies or substantial portions of the Software.
 *
 * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
 * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
 * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
 * AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
 * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
 * OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
 * THE SOFTWARE.
 *
 * ************************************************************************ */

#include "aoclsparse.h"

#include <complex>
#include <cstring>
#include <iomanip>
#include <iostream>
#include <limits>
#include <vector>

// Example to illustrate the usage of aoclsparse_dcsrmm API
int main(void)
{

    aoclsparse_index_base base = aoclsparse_index_base_zero;
    aoclsparse_matrix     A;
    aoclsparse_status     status;
    aoclsparse_order      order = aoclsparse_order_row;
    aoclsparse_mat_descr  descr;

    // By default aoclsparse_create_mat_descr sets aoclsparse_matrix_type to aoclsparse_matrix_type_general
    status = aoclsparse_create_mat_descr(&descr);
    if(status != aoclsparse_status_success)
        return status;

    // Set aoclsparse_index_base to aoclsparse_index_base_zero.
    status = aoclsparse_set_mat_index_base(descr, base);
    if(status != aoclsparse_status_success)
        return status;

    // Matrix sizes
    aoclsparse_int m = 3, n = 3, k = 3;
    aoclsparse_int nnz = 4;
    // Matrix A
    // [ 0.  42.   0.2]
    // [ 4.6  0.   0. ]
    // [ 0.   0.  -8. ]
    aoclsparse_int row_ptr[] = {0, 2, 3, 4};
    aoclsparse_int col_ind[] = {1, 2, 0, 2};
    double         csr_val[] = {42., 0.2, 4.6, -8};

    //Dense Matrix B
    double B[]   = {-1.0, -2.7, 3.0, 4.5, 5.8, -6.0, 1.0, -2.0, 3.0};
    double alpha = 1, beta = 0;
    status = aoclsparse_create_dcsr(&A, base, m, k, nnz, row_ptr, col_ind, csr_val);
    if(status != aoclsparse_status_success)
    {
        return status;
    }

    // expected output matrices
    double C_exp[] = {189.2, 243.2, -251.4, -4.6, -12.42, 13.8, -8, 16, -24};

    // output matrices;
    double               C[9] = {0};
    aoclsparse_operation op   = aoclsparse_operation_none;

    status = aoclsparse_dcsrmm(op, alpha, A, descr, order, B, n, k, beta, C, n);
    if(status != aoclsparse_status_success)
    {
        return status;
    }

    std::cout << std::fixed;
    std::cout.precision(2);
    bool okij, ok_C = true;
    //Initializing precision tolerance range for double
    const double tol = std::sqrt((std::numeric_limits<double>::epsilon()));

    aoclsparse_int i, j;
    std::cout << "The output C matrix\n";
    for(i = 0; i < m; i++)
    {
        for(j = 0; j < n; j++)
        {
            okij = std::abs(C[i * n + j] - C_exp[i * n + j]) <= tol;
            ok_C &= okij;
            std::cout << std::setw(10) << C[i * n + j] << std::setw(3) << (okij ? "" : " !");
        }
        std::cout << "\n";
    }
    aoclsparse_destroy_mat_descr(descr);
    aoclsparse_destroy(&A);

    return ((ok_C ? 0 : 4));
}

Parameters

op – [in] Matrix \(A\) operation type.
alpha – [in] Scalar \(\alpha\).
A – [in] Sparse CSR matrix \(A\) structure.
descr – [in] descriptor of the sparse CSR matrix \(A\). Currently, supports aoclsparse_matrix_type_general, aoclsparse_matrix_type_symmetric, and aoclsparse_matrix_type_hermitian matrices. Both, base-zero and base-one input arrays of CSR matrix are supported.
order – [in] aoclsparse_order_row / aoclsparse_order_column for dense matrix
B – [in] Array of dimension \(ldb \times n\) or \(ldb \times k\) .
n – [in] Number of columns of the dense matrix \(B\) and \(C\).
ldb – [in] Leading dimension of \(B\), must be at least \(\max{(1, k)}\) for \(op(A) = A\), or \(\max{(1, m)}\) when \(op(A) = A^T\) or \(op(A) = A^H\).
beta – [in] Scalar \(\beta\).
C – [inout] Array of dimension \(ldc \times n\).
ldc – [in] Leading dimension of \(C\), must be at least \(\max{(1, m)}\) for \(op(A) = A\), or \(\max{(1, k)}\) when \(op(A) = A^T\) or \(op(A) = A^H\).

Return values

aoclsparse_status_success – The operation completed successfully.
aoclsparse_status_invalid_size – The value of m, n, k, nnz, ldb or ldc is invalid.
aoclsparse_status_invalid_pointer – The pointer descr, A, B, or C is invalid.
aoclsparse_status_invalid_value – The values of descr->base and A->base do not coincide.
aoclsparse_status_not_implemented – aoclsparse_matrix_type is not one of these: aoclsparse_matrix_type_general, aoclsparse_matrix_type_symmetric, aoclsparse_matrix_type_hermitian or input matrix A is not in CSR format

aoclsparse_?csr2m()#

aoclsparse_status aoclsparse_dcsr2m(aoclsparse_operation trans_A, const aoclsparse_mat_descr descrA, const aoclsparse_matrix csrA, aoclsparse_operation trans_B, const aoclsparse_mat_descr descrB, const aoclsparse_matrix csrB, const aoclsparse_request request, aoclsparse_matrix *csrC)#

aoclsparse_status aoclsparse_scsr2m(aoclsparse_operation trans_A, const aoclsparse_mat_descr descrA, const aoclsparse_matrix csrA, aoclsparse_operation trans_B, const aoclsparse_mat_descr descrB, const aoclsparse_matrix csrB, const aoclsparse_request request, aoclsparse_matrix *csrC)#

Sparse matrix Sparse matrix multiplication using CSR storage format for single and double precision datatypes.

aoclsparse_?csr2m multiplies a sparse \(m \times k\) matrix \(A\), defined in CSR storage format, and the sparse \(k \times n\) matrix \(B\), defined in CSR storage format and stores the result to the sparse \(m \times n\) matrix \(C\), such that

\[ C = op(A) \cdot op(B), \]

with

\[\begin{split} op(A) = \left\{ \begin{array}{ll} A, & \text{if } {\bf\mathsf{trans}\_\mathsf{A}} = \text{aoclsparse}\_\text{operation}\_\text{none} \\ A^T, & \text{if } {\bf\mathsf{trans}\_\mathsf{A}} = \text{aoclsparse}\_\text{operation}\_\text{transpose} \\ A^H, & \text{if } {\bf\mathsf{trans}\_\mathsf{A}} = \text{aoclsparse}\_\text{operation}\_\text{conjugate}\_\text{transpose} \end{array} \right. \end{split}\]

and

\[\begin{split} op(B) = \left\{ \begin{array}{ll} B, & \text{if } {\bf\mathsf{trans}\_\mathsf{B}} = \text{aoclsparse}\_\text{operation}\_\text{none} \\ B^T, & \text{if } {\bf\mathsf{trans}\_\mathsf{B}} = \text{aoclsparse}\_\text{operation}\_\text{transpose} \\ B^H, & \text{if } {\bf\mathsf{trans}\_\mathsf{B}} = \text{aoclsparse}\_\text{operation}\_\text{conjugate}\_\text{transpose} \end{array} \right. \end{split}\]

Example (tests/examples/sample_csr2m.cpp)

/* ************************************************************************
 * Copyright (c) 2021-2023 Advanced Micro Devices, Inc. All rights reserved.
 *
 * Permission is hereby granted, free of charge, to any person obtaining a copy
 * of this software and associated documentation files (the "Software"), to deal
 * in the Software without restriction, including without limitation the rights
 * to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
 * copies of the Software, and to permit persons to whom the Software is
 * furnished to do so, subject to the following conditions:
 *
 * The above copyright notice and this permission notice shall be included in
 * all copies or substantial portions of the Software.
 *
 * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
 * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
 * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
 * AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
 * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
 * OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
 * THE SOFTWARE.
 *
 * ************************************************************************ */

#include "aoclsparse.h"

#include <cstdlib>
#include <cstring>
#include <iostream>

/* Computes multiplication of 2 sparse matrices A and B in CSR format. */

// Comment this out for single stage computation
#define TWO_STAGE_COMPUTATION

#define PRINT_OUTPUT

int main(void)
{
    aoclsparse_status     status;
    aoclsparse_int        nnz_C;
    aoclsparse_request    request;
    aoclsparse_index_base base   = aoclsparse_index_base_zero;
    aoclsparse_operation  transA = aoclsparse_operation_none;
    aoclsparse_operation  transB = aoclsparse_operation_none;

    // Print aoclsparse version
    std::cout << aoclsparse_get_version() << std::endl;

    // Initialise matrix descriptor and csr matrix structure of inputs A and B
    aoclsparse_mat_descr descrA;
    aoclsparse_mat_descr descrB;
    aoclsparse_matrix    csrA;
    aoclsparse_matrix    csrB;

    // Create matrix descriptor of input matrices
    // aoclsparse_create_mat_descr set aoclsparse_matrix_type to aoclsparse_matrix_type_general
    // and aoclsparse_index_base to aoclsparse_index_base_zero.
    aoclsparse_create_mat_descr(&descrA);
    aoclsparse_create_mat_descr(&descrB);

    // Matrix sizes
    aoclsparse_int m = 3, n = 3, k = 3;
    aoclsparse_int nnz_A = 6, nnz_B = 4;
    // Matrix A
    // 	1  0  2
    // 	0  0  3
    // 	4  5  6
    aoclsparse_int row_ptr_A[] = {0, 2, 3, 6};
    aoclsparse_int col_ind_A[] = {0, 2, 2, 0, 1, 2};
    float          val_A[]     = {1, 2, 3, 4, 5, 6};
    aoclsparse_create_scsr(&csrA, base, m, k, nnz_A, row_ptr_A, col_ind_A, val_A);

    // Matrix B
    // 	1  2  0
    // 	0  0  3
    // 	0  4  0
    aoclsparse_int row_ptr_B[] = {0, 2, 3, 4};
    aoclsparse_int col_ind_B[] = {0, 1, 2, 1};
    float          val_B[]     = {1, 2, 3, 4};
    aoclsparse_create_scsr(&csrB, base, k, n, nnz_B, row_ptr_B, col_ind_B, val_B);

    aoclsparse_matrix csrC          = NULL;
    aoclsparse_int   *csr_row_ptr_C = NULL;
    aoclsparse_int   *csr_col_ind_C = NULL;
    float            *csr_val_C     = NULL;
    aoclsparse_int    C_M, C_N;

#ifdef TWO_STAGE_COMPUTATION
    std::cout << "Invoking aoclsparse_scsr2m with aoclsparse_stage_nnz_count..";
    // aoclsparse_stage_nnz_count : Only rowIndex array of the CSR matrix
    // is computed internally. The output sparse CSR matrix can be
    // extracted to measure the memory required for full operation.
    request = aoclsparse_stage_nnz_count;
    status  = aoclsparse_scsr2m(transA, descrA, csrA, transB, descrB, csrB, request, &csrC);
    if(status == aoclsparse_status_success)
        std::cout << "DONE\n";
    else
    {
        std::cout << "ERROR in aoclsparse_scsr2m\n";
        exit(EXIT_FAILURE);
    }

    aoclsparse_export_scsr(
        csrC, &base, &C_M, &C_N, &nnz_C, &csr_row_ptr_C, &csr_col_ind_C, &csr_val_C);

#if 0
    std::cout << "C_M"
              << "\t"
              << "C_N"
              << "\t"
              << "nnz_C" << std::endl;
    std::cout << C_M << "\t" << C_N << "\t" << nnz_C << std::endl;

    std::cout << "csr_row_ptr_C: " << std::endl;
    for(aoclsparse_int i = 0; i < C_M + 1; i++)
        std::cout << csr_row_ptr_C[i] << "\t";
    std::cout << std::endl;
#endif

    std::cout << "Invoking aoclsparse_scsr2m with aoclsparse_stage_finalize..";
    // aoclsparse_stage_finalize : Finalize computation of remaining
    // output arrays ( column indices and values of output matrix entries) .
    // Has to be called only after aoclsparse_scsr2m call with
    // aoclsparse_stage_nnz_count parameter.
    request = aoclsparse_stage_finalize;
    status  = aoclsparse_scsr2m(transA, descrA, csrA, transB, descrB, csrB, request, &csrC);
    if(status == aoclsparse_status_success)
        std::cout << "DONE\n";
    else
    {
        std::cout << "ERROR in aoclsparse_scsr2m\n";
        exit(EXIT_FAILURE);
    }
    aoclsparse_export_scsr(
        csrC, &base, &C_M, &C_N, &nnz_C, &csr_row_ptr_C, &csr_col_ind_C, &csr_val_C);

#ifdef PRINT_OUTPUT
    std::cout << "csr_col_ind_C: " << std::endl;
    for(aoclsparse_int i = 0; i < nnz_C; i++)
        std::cout << csr_col_ind_C[i] << "\t";
    std::cout << std::endl;

    std::cout << "csr_val_C: " << std::endl;
    for(aoclsparse_int i = 0; i < nnz_C; i++)
        std::cout << csr_val_C[i] << "\t";
    std::cout << std::endl;
#endif

#else // SINGLE STAGE
    std::cout << "Invoking aoclsparse_scsr2m with aoclsparse_stage_full_computation..";
    // aoclsparse_stage_full_computation :  Whole computation is performed in
    // single step.
    request = aoclsparse_stage_full_computation;
    status  = aoclsparse_scsr2m(transA, descrA, csrA, transB, descrB, csrB, request, &csrC);
    if(status == aoclsparse_status_success)
        std::cout << "DONE\n";
    else
    {
        std::cout << "ERROR in aoclsparse_scsr2m\n";
        exit(EXIT_FAILURE);
    }

    aoclsparse_export_scsr(
        csrC, &base, &C_M, &C_N, &nnz_C, &csr_row_ptr_C, &csr_col_ind_C, &csr_val_C);

#ifdef PRINT_OUTPUT
    std::cout << "C_M"
              << "\t"
              << "C_N"
              << "\t"
              << "nnz_C" << std::endl;
    std::cout << C_M << "\t" << C_N << "\t" << nnz_C << std::endl;

    std::cout << "csr_row_ptr_C: " << std::endl;
    for(aoclsparse_int i = 0; i < C_M + 1; i++)
        std::cout << csr_row_ptr_C[i] << "\t";
    std::cout << std::endl;

    std::cout << "csr_col_ind_C: " << std::endl;
    for(aoclsparse_int i = 0; i < nnz_C; i++)
        std::cout << csr_col_ind_C[i] << "\t";
    std::cout << std::endl;

    std::cout << "csr_val_C: " << std::endl;
    for(aoclsparse_int i = 0; i < nnz_C; i++)
        std::cout << csr_val_C[i] << "\t";
    std::cout << std::endl;
#endif

#endif
    aoclsparse_destroy_mat_descr(descrA);
    aoclsparse_destroy_mat_descr(descrB);
    aoclsparse_destroy(&csrA);
    aoclsparse_destroy(&csrB);
    aoclsparse_destroy(&csrC);
    return 0;
}

Parameters

trans_A – [in] matrix \(A\) operation type.
descrA – [in] descriptor of the sparse CSR matrix \(A\). Currently, only aoclsparse_matrix_type_general is supported.
csrA – [in] sparse CSR matrix \(A\) structure.
trans_B – [in] matrix \(B\) operation type.
descrB – [in] descriptor of the sparse CSR matrix \(B\). Currently, only aoclsparse_matrix_type_general is supported.
csrB – [in] sparse CSR matrix \(B\) structure.
request – [in] Specifies full computation or two-stage algorithm aoclsparse_stage_nnz_count , Only rowIndex array of the CSR matrix is computed internally. The output sparse CSR matrix can be extracted to measure the memory required for full operation. aoclsparse_stage_finalize . Finalize computation of remaining output arrays ( column indices and values of output matrix entries) . Has to be called only after aoclsparse_dcsr2m() call with aoclsparse_stage_nnz_count parameter. aoclsparse_stage_full_computation. Perform the entire computation in a single step.
*csrC – [out] Pointer to sparse CSR matrix \(C\) structure.

Return values

aoclsparse_status_success – the operation completed successfully.
aoclsparse_status_invalid_size – input parameters contain an invalid value.
aoclsparse_status_invalid_pointer – descrA, csr, descrB, csrB, csrC is invalid.
aoclsparse_status_not_implemented – aoclsparse_matrix_type is not aoclsparse_matrix_type_general or input matrices A or B is not in CSR format

aoclsparse_?add()#

aoclsparse_status aoclsparse_sadd(const aoclsparse_operation op, const aoclsparse_matrix A, const float alpha, const aoclsparse_matrix B, aoclsparse_matrix *C)#

aoclsparse_status aoclsparse_dadd(const aoclsparse_operation op, const aoclsparse_matrix A, const double alpha, const aoclsparse_matrix B, aoclsparse_matrix *C)#

aoclsparse_status aoclsparse_cadd(const aoclsparse_operation op, const aoclsparse_matrix A, const aoclsparse_float_complex alpha, const aoclsparse_matrix B, aoclsparse_matrix *C)#

aoclsparse_status aoclsparse_zadd(const aoclsparse_operation op, const aoclsparse_matrix A, const aoclsparse_double_complex alpha, const aoclsparse_matrix B, aoclsparse_matrix *C)#

Addition of two sparse matrices.

aoclsparse_?add adds two sparse matrices and returns a sparse matrix. Matrices can be either real or complex types but cannot be intermixed. It performs

\[\begin{split} C = \alpha \, op(A) + B \qquad\text{ with } \qquad op(A) = \left\{ \begin{array}{ll} A, & \text{ if } {\bf\mathsf{op}} = \text{aoclsparse}\_\text{operation}\_\text{none} \\ A^T, & \text{ if } {\bf\mathsf{op}} = \text{aoclsparse}\_\text{operation}\_\text{transpose} \\ A^H, & \text{ if } {\bf\mathsf{op}} = \text{aoclsparse}\_\text{operation}\_\text{conjugate}\_\text{transpose} \end{array} \right. \end{split}\]

where \(A\) is a \(m \times n\) matrix and \(B\) is a \(m \times n\) matrix, if op = aoclsparse_operation_none. Otherwise \(A\) is \(n \times m\) and the result matrix \(C\) has the same dimension as \(B\).

Note

Only matrices in CSR format are supported in this release.

Parameters

op – [in] matrix \(A\) operation type.
alpha – [in] scalar with same precision as \(A\) and \(B\) matrix
A – [in] source sparse matrix \(A\)
B – [in] source sparse matrix \(B\)
*C – [out] pointer to the sparse output matrix \(C\)

Return values

aoclsparse_status_success – The operation completed successfully.
aoclsparse_status_invalid_pointer – A or B or C are invalid
aoclsparse_status_invalid_size – The dimensions of A and B are not compatible.
aoclsparse_status_internal_error – Internal Error Occured
aoclsparse_status_memory_error – Memory allocation failure.
aoclsparse_status_not_implemented – Matrices are not in CSR format.

aoclsparse_?spmmd()#

aoclsparse_status aoclsparse_sspmmd(const aoclsparse_operation op, const aoclsparse_matrix A, const aoclsparse_matrix B, const aoclsparse_order layout, float *C, const aoclsparse_int ldc)#

aoclsparse_status aoclsparse_dspmmd(const aoclsparse_operation op, const aoclsparse_matrix A, const aoclsparse_matrix B, const aoclsparse_order layout, double *C, const aoclsparse_int ldc)#

aoclsparse_status aoclsparse_cspmmd(const aoclsparse_operation op, const aoclsparse_matrix A, const aoclsparse_matrix B, const aoclsparse_order layout, aoclsparse_float_complex *C, const aoclsparse_int ldc)#

aoclsparse_status aoclsparse_zspmmd(const aoclsparse_operation op, const aoclsparse_matrix A, const aoclsparse_matrix B, const aoclsparse_order layout, aoclsparse_double_complex *C, const aoclsparse_int ldc)#

Matrix multiplication of two sparse matrices stored in the CSR storage format. The output matrix is stored in a dense format.

aoclsparse_?spmmd multiplies a sparse matrix \(A\) and a sparse matrix \(B\), both stored in the CSR storage format, and saves the result in a dense matrix \(C\), such that

\[ C := op(A) \cdot B, \]

with

\[\begin{split} op(A) = \left\{ \begin{array}{ll} A, & \text{if op} = \text{aoclsparse\_operation\_none} \\ A^T, & \text{if op} = \text{aoclsparse\_operation\_transpose} \\ A^H, & \text{if op} = \text{aoclsparse\_operation\_conjugate\_transpose} \end{array} \right. \end{split}\]

Parameters

op – [in] Operation to perform on matrix \(A\).
A – [in] Matrix structure containing sparse matrix \(A\) of size \(m \times k\).
B – [in] Matrix structure containing sparse matrix \(B\) of size \(k \times n\) if op is aoclsparse_operation_none otherwise of size \(m \times n\).
layout – [in] Ordering of the dense output matrix: valid values are aoclsparse_order_row and aoclsparse_order_column.
C – [inout] Dense output matrix \(C\) of size \(m \times n\) if op is aoclsparse_operation_none, otherwise of size \(k \times n\) containing the matrix-matrix product of \(A\) and \(B\).
ldc – [in] Leading dimension of \(C\), e.g., for C stored in aoclsparse_order_row, ldc must be at least \(\max{(1, m)}\) when \(op(A) = A\), or \(\max{(1, k)}\) if \(op(A) = A^T\) or \(op(A) = A^H\).

Return values

aoclsparse_status_success – The operation completed successfully.
aoclsparse_status_invalid_size – m, n, k, nnz or ldc is not valid.
aoclsparse_status_invalid_pointer – A, B or C pointer is not valid.
aoclsparse_status_wrong_type – aoclsparse_matrix_data_type does not match the precision type.
aoclsparse_status_not_implemented – aoclsparse_matrix_format_type is not aoclsparse_csr_mat.

aoclsparse_?sp2md()#

aoclsparse_status aoclsparse_ssp2md(const aoclsparse_operation opA, const aoclsparse_mat_descr descrA, const aoclsparse_matrix A, const aoclsparse_operation opB, const aoclsparse_mat_descr descrB, const aoclsparse_matrix B, const float alpha, const float beta, float *C, const aoclsparse_order layout, const aoclsparse_int ldc)#

aoclsparse_status aoclsparse_dsp2md(const aoclsparse_operation opA, const aoclsparse_mat_descr descrA, const aoclsparse_matrix A, const aoclsparse_operation opB, const aoclsparse_mat_descr descrB, const aoclsparse_matrix B, const double alpha, const double beta, double *C, const aoclsparse_order layout, const aoclsparse_int ldc)#

aoclsparse_status aoclsparse_csp2md(const aoclsparse_operation opA, const aoclsparse_mat_descr descrA, const aoclsparse_matrix A, const aoclsparse_operation opB, const aoclsparse_mat_descr descrB, const aoclsparse_matrix B, aoclsparse_float_complex alpha, aoclsparse_float_complex beta, aoclsparse_float_complex *C, const aoclsparse_order layout, const aoclsparse_int ldc)#

aoclsparse_status aoclsparse_zsp2md(const aoclsparse_operation opA, const aoclsparse_mat_descr descrA, const aoclsparse_matrix A, const aoclsparse_operation opB, const aoclsparse_mat_descr descrB, const aoclsparse_matrix B, aoclsparse_double_complex alpha, aoclsparse_double_complex beta, aoclsparse_double_complex *C, const aoclsparse_order layout, const aoclsparse_int ldc)#

A variant of matrix multiplication of two sparse matrices stored in the CSR storage format. The output matrix is stored in a dense format. Supports operations on both sparse matrices.

aoclsparse_?sp2md multiplies a sparse matrix \(A\) and a sparse matrix \(B\), both stored in the CSR storage format, and saves the result in a dense matrix \(C\), such that

\[ C := \alpha \cdot op(A) \cdot op(B) + \beta \cdot C, \]

with

\[\begin{split} op(A) = \left\{ \begin{array}{ll} A, & \text{if opA} = \text{aoclsparse\_operation\_none} \\ A^T, & \text{if opA} = \text{aoclsparse\_operation\_transpose} \\ A^H, & \text{if opA} = \text{aoclsparse\_operation\_conjugate\_transpose} \end{array} \right. \end{split}\]

and

\[\begin{split} op(B) = \left\{ \begin{array}{ll} B, & \text{if opB} = \text{aoclsparse\_operation\_none} \\ B^T, & \text{if opB} = \text{aoclsparse\_operation\_transpose} \\ B^H, & \text{if opB} = \text{aoclsparse\_operation\_conjugate\_transpose} \end{array} \right. \end{split}\]

Parameters

opA – [in] Operation to perform on matrix \(A\).
descrA – [in] Descriptor of A. Only aoclsparse_matrix_type_general is supported at present. As a consequence, all other parameters within the descriptor are ignored.
A – [in] Matrix structure containing sparse matrix \(A\) of size \(m \times k\).
opB – [in] Operation to perform on matrix \(B\).
descrB – [in] Descriptor of B. Only aoclsparse_matrix_type_general is supported at present. As a consequence, all other parameters within the descriptor are ignored.
B – [in] Matrix structure containing sparse matrix \(B\) of size \(k \times n\) if op is aoclsparse_operation_none otherwise of size \(m \times n\).
alpha – [in] Value of \( \alpha\).
beta – [in] Value of \( \beta\).
C – [inout] Dense output matrix \(C\).
layout – [in] Ordering of the dense output matrix: valid values are aoclsparse_order_row and aoclsparse_order_column.
ldc – [in] Leading dimension of \(C\), e.g., for C stored in aoclsparse_order_row, ldc must be at least \(\max{(1, m)}\) ( \(op(A) = A\)) or \(\max{(1, k)}\) ( \(op(A) = A^T\) or \(op(A) = A^H\)).

Return values

aoclsparse_status_success – The operation completed successfully.
aoclsparse_status_invalid_size – m, n, k, nnz or ldc is not valid.
aoclsparse_status_invalid_pointer – A, B or C pointer is not valid.
aoclsparse_status_wrong_type – aoclsparse_matrix_data_type does not match the precision type.
aoclsparse_status_not_implemented – aoclsparse_matrix_format_type is not aoclsparse_csr_mat.
aoclsparse_status_internal_error – An internal error occurred.

aoclsparse_syrk()#

aoclsparse_status aoclsparse_syrk(const aoclsparse_operation opA, const aoclsparse_matrix A, aoclsparse_matrix *C)#

Multiplication of a sparse matrix and its transpose (or conjugate transpose) stored as a sparse matrix.

aoclsparse_syrk multiplies a sparse matrix with its transpose (or conjugate transpose) in CSR storage format. The result is stored in a newly allocated sparse matrix in CSR format, such that

\[ C := A \cdot op(A) \]

if \(opA\) is aoclsparse_operation_none.

Otherwise,

\[ C := op(A) \cdot A, \]

where

\[\begin{split} op(A) = \left\{ \begin{array}{ll} A^T, & \text{transpose of } {\bf\mathsf{A}} \text{ for real matrices}\\ A^H, & \text{conjugate transpose of } {\bf\mathsf{A}} \text{ for complex matrices}\\ \end{array} \right. \end{split}\]

where \(A\) is a \(m \times n\) matrix, opA is one of aoclsparse_operation_none, aoclsparse_operation_transpose (for real matrices) or aoclsparse_operation_conjugate_transpose (for complex matrices). The output matrix \(C\) is a sparse symmetric (or Hermitian) matrix stored as an upper triangular matrix in CSR format.

Example (tests/examples/sample_dsyrk.cpp)

/* ************************************************************************
 * Copyright (c) 2024 Advanced Micro Devices, Inc.
 *
 * Permission is hereby granted, free of charge, to any person obtaining a copy
 * of this software and associated documentation files (the "Software"), to deal
 * in the Software without restriction, including without limitation the rights
 * to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
 * copies of the Software, and to permit persons to whom the Software is
 * furnished to do so, subject to the following conditions:
 *
 * The above copyright notice and this permission notice shall be included in
 * all copies or substantial portions of the Software.
 *
 * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
 * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
 * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
 * AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
 * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
 * OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
 * THE SOFTWARE.
 *
 * ************************************************************************ */

#include "aoclsparse.h"

#include <complex>
#include <cstring>
#include <iomanip>
#include <iostream>
#include <limits>
#include <vector>

// An example to illustrate the usage of aoclsparse_syrk
int main(void)
{
    std::cout << "-------------------------------" << std::endl
              << "----- syrk sample program -----" << std::endl
              << "-------------------------------" << std::endl
              << std::endl;

    aoclsparse_matrix     A;
    aoclsparse_matrix     C;
    aoclsparse_status     status;
    aoclsparse_index_base base          = aoclsparse_index_base_zero;
    aoclsparse_int       *csr_row_ptr_C = NULL;
    aoclsparse_int       *csr_col_ind_C = NULL;
    double               *csr_val_C     = NULL;
    aoclsparse_int        C_m, C_n, C_nnz;

    // Matrix sizes
    aoclsparse_int m = 4, k = 3;
    aoclsparse_int nnz = 7;
    // Matrix A
    // [ 0.  -1.2   2.3]
    // [ 4.6  0.    0. ]
    // [ 0.   3.0  -8.1 ]
    // [ 0.3   0.  -5.1 ]
    aoclsparse_int       row_ptr[] = {0, 2, 3, 5, 7};
    aoclsparse_int       col_ind[] = {1, 2, 0, 1, 2, 0, 2};
    double               csr_val[] = {-1.2, 2.3, 4.6, 3.0, -8.1, 0.3, -5.1};
    aoclsparse_operation op        = aoclsparse_operation_none;

    status = aoclsparse_create_dcsr(&A, base, m, k, nnz, row_ptr, col_ind, csr_val);
    if(status != aoclsparse_status_success)
    {
        return status;
    }

    // expected output matrices
    aoclsparse_int C_m_exp = 4, C_n_exp = 4, C_nnz_exp = 8;
    aoclsparse_int csr_row_ptr_C_exp[] = {0, 3, 5, 7, 8};
    aoclsparse_int csr_col_ind_C_exp[] = {0, 2, 3, 1, 3, 2, 3, 3};
    double         csr_val_C_exp[]
        = {6.73000, -22.23000, -11.73000, 21.16000, 1.38000, 74.61000, 41.31000, 26.10000};

    // output matrices;

    status = aoclsparse_syrk(op, A, &C);
    if(status != aoclsparse_status_success)
    {
        return status;
    }
    aoclsparse_export_dcsr(
        C, &base, &C_m, &C_n, &C_nnz, &csr_row_ptr_C, &csr_col_ind_C, &csr_val_C);

    bool oka, okb, okc, oki, okj, okk, ok = true;
    std::cout << std::endl
              << "Output Matrix C: " << std::endl
              << std::setw(11) << "C_m" << std::setw(3) << "" << std::setw(11) << "expected"
              << std::setw(3) << "" << std::setw(11) << "C_n" << std::setw(3) << "" << std::setw(11)
              << "expected" << std::setw(3) << "" << std::setw(11) << "C_nnz" << std::setw(3) << ""
              << std::setw(11) << "expected" << std::endl;
    oka = C_m == C_m_exp;
    ok &= oka;
    std::cout << std::setw(11) << C_m << std::setw(3) << "" << std::setw(11) << C_m_exp
              << std::setw(2) << (oka ? "" : " !");
    okb = C_n == C_n_exp;
    ok &= okb;
    std::cout << std::setw(11) << C_n << std::setw(3) << "" << std::setw(11) << C_n_exp
              << std::setw(2) << (okb ? "" : " !");
    okc = C_nnz == C_nnz_exp;
    ok &= okc;
    std::cout << std::setw(11) << C_nnz << std::setw(3) << "" << std::setw(11) << C_nnz_exp
              << std::setw(2) << (okc ? "" : " !");
    std::cout << std::endl;
    std::cout << std::endl;
    std::cout << std::setw(11) << "csr_val_C" << std::setw(3) << "" << std::setw(11) << "expected"
              << std::setw(3) << "" << std::setw(11) << "csr_col_ind_C" << std::setw(3) << ""
              << std::setw(11) << "expected" << std::setw(3) << "" << std::setw(11)
              << "csr_row_ptr_C" << std::setw(3) << "" << std::setw(11) << "expected" << std::endl;
    //Initializing precision tolerance range for double
    const double tol = 1e-03;
    for(aoclsparse_int i = 0; i < C_nnz; i++)
    {
        oki = ((std::abs(csr_val_C[i] - csr_val_C_exp[i]) <= tol));
        ok &= oki;
        std::cout << std::setw(11) << csr_val_C[i] << "" << std::setw(11) << csr_val_C_exp[i]
                  << (oki ? "" : " !");
        okj = csr_col_ind_C[i] == csr_col_ind_C_exp[i];
        ok &= okj;
        std::cout << std::setw(11) << csr_col_ind_C[i] << std::setw(3) << "" << std::setw(11)
                  << csr_col_ind_C_exp[i] << std::setw(2) << (okj ? "" : " !");
        if(i <= C_m)
        {
            okk = csr_row_ptr_C[i] == csr_row_ptr_C_exp[i];
            ok &= okk;
            std::cout << " " << std::setw(11) << csr_row_ptr_C[i] << std::setw(3) << ""
                      << std::setw(11) << csr_row_ptr_C_exp[i] << std::setw(2) << (okk ? "" : " !");
        }
        std::cout << std::endl;
    }
    aoclsparse_destroy(&A);
    aoclsparse_destroy(&C);

    return (ok ? 0 : 6);
}

Note

aoclsparse_syrk assumes that the input CSR matrix has sorted column indices in each row. If not, call aoclsparse_order_mat() before calling aoclsparse_syrk.

Note

aoclsparse_syrk currently does not support aoclsparse_operation_transpose for complex A.

Parameters

opA – [in] Matrix \(A\) operation type.
A – [in] Sorted sparse CSR matrix \(A\).
*C – [out] Pointer to the new sparse CSR symmetric/Hermitian matrix \(C\). Only upper triangle of the result matrix is computed. The column indices of the output matrix in CSR format might be unsorted. The matrix should be freed by aoclsparse_destroy() when no longer needed.

Return values

aoclsparse_status_success – The operation completed successfully.
aoclsparse_status_invalid_pointer – A, C is invalid.
aoclsparse_status_wrong_type – A and its operation type do not match.
aoclsparse_status_not_implemented – The input matrix is not in the CSR format or opA is aoclsparse_operation_transpose and A has complex values.
aoclsparse_status_invalid_value – The value of opA is invalid.
aoclsparse_status_unsorted_input – Input matrices are not sorted.
aoclsparse_status_memory_error – Memory allocation failure.

aoclsparse_?syrkd()#

aoclsparse_status aoclsparse_ssyrkd(const aoclsparse_operation opA, const aoclsparse_matrix A, const float alpha, const float beta, float *C, const aoclsparse_order orderC, const aoclsparse_int ldc)#

aoclsparse_status aoclsparse_dsyrkd(const aoclsparse_operation opA, const aoclsparse_matrix A, const double alpha, const double beta, double *C, const aoclsparse_order orderC, const aoclsparse_int ldc)#

aoclsparse_status aoclsparse_csyrkd(const aoclsparse_operation opA, const aoclsparse_matrix A, const aoclsparse_float_complex alpha, const aoclsparse_float_complex beta, aoclsparse_float_complex *C, const aoclsparse_order orderC, const aoclsparse_int ldc)#

aoclsparse_status aoclsparse_zsyrkd(const aoclsparse_operation opA, const aoclsparse_matrix A, const aoclsparse_double_complex alpha, const aoclsparse_double_complex beta, aoclsparse_double_complex *C, const aoclsparse_order orderC, const aoclsparse_int ldc)#

Multiplication of a sparse matrix and its transpose (or conjugate transpose) for all data types.

aoclsparse_syrkd multiplies a sparse matrix with its transpose (or conjugate transpose) in CSR storage format. The result is stored in a dense format, such that

\[ C := \alpha \cdot A \cdot op(A) + \beta \cdot C \]

if \(opA\) is aoclsparse_operation_none.

Otherwise,

\[ C := \alpha \cdot op(A) \cdot A + \beta \cdot C \]

where \(A\) is a \(m \times n\) sparse matrix, opA is one of aoclsparse_operation_none, aoclsparse_operation_transpose (for real matrices) or aoclsparse_operation_conjugate_transpose (for complex matrices). The output matrix \(C\) is a dense symmetric (or Hermitian) matrix stored as an upper triangular matrix.

Example (tests/examples/sample_dsyrkd.cpp)

/* ************************************************************************
 * Copyright (c) 2024 Advanced Micro Devices, Inc.
 *
 * Permission is hereby granted, free of charge, to any person obtaining a copy
 * of this software and associated documentation files (the "Software"), to deal
 * in the Software without restriction, including without limitation the rights
 * to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
 * copies of the Software, and to permit persons to whom the Software is
 * furnished to do so, subject to the following conditions:
 *
 * The above copyright notice and this permission notice shall be included in
 * all copies or substantial portions of the Software.
 *
 * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
 * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
 * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
 * AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
 * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
 * OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
 * THE SOFTWARE.
 *
 * ************************************************************************ */

#include "aoclsparse.h"

#include <complex>
#include <cstring>
#include <iomanip>
#include <iostream>
#include <limits>
#include <vector>

// An example to illustrate the usage of aoclsparse_syrk
int main(void)
{
    std::cout << "-------------------------------" << std::endl
              << "----- syrkd sample program -----" << std::endl
              << "-------------------------------" << std::endl
              << std::endl;

    aoclsparse_matrix     A;
    aoclsparse_status     status;
    aoclsparse_index_base base = aoclsparse_index_base_zero;
    aoclsparse_mat_descr  descr;
    double                alpha  = 2.1;
    double                beta   = 1.3;
    aoclsparse_int        ldc    = 10;
    aoclsparse_operation  op     = aoclsparse_operation_none;
    aoclsparse_order      layout = aoclsparse_order_row;
    aoclsparse_int        m_C;
    double               *C     = nullptr;
    double               *C_exp = nullptr;

    // By default aoclsparse_create_mat_descr sets aoclsparse_matrix_type to aoclsparse_matrix_type_general
    status = aoclsparse_create_mat_descr(&descr);
    if(status != aoclsparse_status_success)
        return status;

    // Matrix sizes
    aoclsparse_int m = 4, k = 3;
    aoclsparse_int nnz = 7;
    // Matrix A
    // [ 0.  -1.2   2.3]
    // [ 4.6  0.    0. ]
    // [ 0.   3.0  -8.1 ]
    // [ 0.3   0.  -5.1 ]
    aoclsparse_int row_ptr[] = {0, 2, 3, 5, 7};
    aoclsparse_int col_ind[] = {1, 2, 0, 1, 2, 0, 2};
    double         csr_val[] = {-1.2, 2.3, 4.6, 3.0, -8.1, 0.3, -5.1};

    status = aoclsparse_create_dcsr(&A, base, m, k, nnz, row_ptr, col_ind, csr_val);
    if(status != aoclsparse_status_success)
    {
        return status;
    }

    // output matrix dimension
    m_C = 4;

    C     = (double *)malloc(sizeof(double) * ldc * m_C);
    C_exp = (double *)malloc(sizeof(double) * ldc * m_C);

    // set the large C matrix to -1, only a window of size 4x4 (upper triangular part)
    // is updated with the syrkd call
    for(aoclsparse_int i = 0; i < ldc * m_C; i++)
    {
        C[i]     = -1;
        C_exp[i] = -1;
    }

    // Expected value of C which is C_exp
    // [ 12.8330    -1.3000   -47.9830   -25.9330,  -1, -1, ....]
    // [ -1         43.1360    -1.3000     1.5980,  -1, -1, ....]
    // [ -1         -1         155.3810    85.4510, -1, -1, ....]
    // [ -1         -1         -1          53.5100, -1, -1, ....]
    // [ -1         -1         -1          -1,      -1, -1, ....]
    // ...
    C_exp[0]  = 12.8330;
    C_exp[1]  = -1.3000;
    C_exp[2]  = -47.9830;
    C_exp[3]  = -25.9330;
    C_exp[11] = 43.1360;
    C_exp[12] = -1.3000;
    C_exp[13] = 1.5980;
    C_exp[22] = 155.3810;
    C_exp[23] = 85.4510;
    C_exp[33] = 53.5100;

    status = aoclsparse_dsyrkd(op, A, alpha, beta, C, layout, ldc);
    if(status != aoclsparse_status_success)
    {
        return status;
    }

    bool oki, ok = true;
    //Initializing precision tolerance range for double
    const double tol = 1e-03;
    std::cout << "Expected output\n"
              << "-------------------------------\n";
    for(aoclsparse_int i = 0; i < m_C; i++)
    {
        for(aoclsparse_int j = 0; j < ldc; j++)
        {
            std::cout << std::setw(10) << C_exp[i * ldc + j] << std::setw(3) << "";
        }
        std::cout << std::endl;
    }

    std::cout << "Actual output\n"
              << "-------------------------------\n";
    for(aoclsparse_int i = 0; i < m_C; i++)
    {
        for(aoclsparse_int j = 0; j < ldc; j++)
        {
            oki = ((std::abs(C[i * ldc + j] - C_exp[i * ldc + j]) <= tol));
            ok &= oki;
            std::cout << std::setw(10) << C[i * ldc + j] << std::setw(3) << (oki ? "" : " !");
        }
        std::cout << std::endl;
    }
    aoclsparse_destroy_mat_descr(descr);
    aoclsparse_destroy(&A);
    free(C);
    free(C_exp);

    return (ok ? 0 : 6);
}

Note

aoclsparse_syrkd assumes that the input CSR matrix has sorted column indices in each row. If not, call aoclsparse_order_mat() before calling aoclsparse_syrkd.

Note

For complex type, only the real parts of \(\alpha\) and \(\beta\) are taken into account to preserve Hermitian \(C\).

Note

aoclsparse_syrkd currently does not support aoclsparse_operation_transpose for complex A.

Parameters

opA – [in] Matrix \(A\) operation type.
A – [in] Sorted sparse CSR matrix \(A\).
alpha – [in] Scalar \(\alpha\).
beta – [in] Scalar \(\beta\).
C – [inout] Output dense matrix. Only upper triangular part of the matrix is processed during the computation, the strictly lower triangle is not modified.
orderC – [in] Storage format of the output dense matrix, \(C\). It can be aoclsparse_order_row or aoclsparse_order_column.
ldc – [in] Leading dimension of \(C\).

Return values

aoclsparse_status_success – The operation completed successfully.
aoclsparse_status_invalid_pointer – A, C is invalid.
aoclsparse_status_wrong_type – A and its operation type do not match.
aoclsparse_status_not_implemented – The input matrix is not in the CSR format or opA is aoclsparse_operation_transpose and A has complex values.
aoclsparse_status_invalid_value – The value of opA, orderC or ldc is invalid.
aoclsparse_status_unsorted_input – Input matrix is not sorted.
aoclsparse_status_memory_error – Memory allocation failure.

aoclsparse_?sypr()#

aoclsparse_status aoclsparse_sypr(aoclsparse_operation opA, const aoclsparse_matrix A, const aoclsparse_matrix B, const aoclsparse_mat_descr descrB, aoclsparse_matrix *C, const aoclsparse_request request)#

Symmetric product of three sparse matrices for real and complex datatypes stored as a sparse matrix.

aoclsparse_sypr multiplies three sparse matrices in CSR storage format. The result is returned in a newly allocated symmetric or Hermitian sparse matrix stored as an upper triangle in CSR format.

If \(opA\) is aoclsparse_operation_none,

\[ C = A \cdot B \cdot A^T, \]

\[ C = A \cdot B \cdot A^H, \]

for real or complex input matrices, respectively, where \(A\) is a \(m \times n\) general matrix , \(B\) is a \(n \times n\) symmetric (for real data types) or Hermitian (for complex data types) matrix, resulting in a symmetric or Hermitian \(m \times m\) matrix \(C\).

Otherwise,

\[ C = op(A) \cdot B \cdot A, \]

with

\[\begin{split} op(A) = \left\{ \begin{array}{ll} A^T, & \text{if } {\bf\mathsf{opA}} = \text{aoclsparse}\_\text{operation}\_\text{transpose} \\ A^H, & \text{if } {\bf\mathsf{opA}} = \text{aoclsparse}\_\text{operation}\_\text{conjugate}\_\text{transpose} \end{array} \right. \end{split}\]

where \(A\) is a \(m \times n\) matrix and \(B\) is a \(m \times m\) symmetric (or Hermitian) matrix, resulting in a \(n \times n\) symmetric (or Hermitian) matrix \(C\).

Depending on request, aoclsparse_sypr might compute the result in a single stage (aoclsparse_stage_full_computation) or in two stages. Then the first stage (aoclsparse_stage_nnz_count) allocates memory for the new output matrix \(C\) and computes its number of non-zeros and their structure which is followed by the second stage (aoclsparse_stage_finalize) to compute the column indices and values of all elements. The second stage can be invoked multiple times (either after aoclsparse_stage_full_computation or aoclsparse_stage_nnz_count) to recompute the numerical values of \(C\) on assumption that the sparsity structure of the input matrices remained unchanged and only the values of the non-zero elements were modified (e.g., by a call to aoclsparse_supdate_values() and variants).

Example (tests/examples/sample_zsypr.cpp)

/* ************************************************************************
 * Copyright (c) 2024 Advanced Micro Devices, Inc. All rights reserved.
 *
 * Permission is hereby granted, free of charge, to any person obtaining a copy
 * of this software and associated documentation files (the "Software"), to deal
 * in the Software without restriction, including without limitation the rights
 * to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
 * copies of the Software, and to permit persons to whom the Software is
 * furnished to do so, subject to the following conditions:
 *
 * The above copyright notice and this permission notice shall be included in
 * all copies or substantial portions of the Software.
 *
 * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
 * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
 * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
 * AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
 * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
 * OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
 * THE SOFTWARE.
 *
 * ************************************************************************ */

#include "aoclsparse.h"

#include <complex>
#include <iomanip>
#include <iostream>

/* Computes triple symmetric product of complex sparse matrices
 * C:=opA(A)*B*A, where opA = conjugate transpose
*/

// Comment this out for single stage computation
//#define TWO_STAGE_COMPUTATION

int main(void)
{
    std::cout << "-------------------------------" << std::endl
              << "----- SYPR sample program -----" << std::endl
              << "-------------------------------" << std::endl
              << std::endl;

    aoclsparse_status     status;
    aoclsparse_int        nnz_C;
    aoclsparse_request    request;
    aoclsparse_index_base base = aoclsparse_index_base_zero;
    aoclsparse_operation  opA  = aoclsparse_operation_conjugate_transpose;

    // Print aoclsparse version
    std::cout << aoclsparse_get_version() << std::endl;

    // Initialise matrix descriptor and csr matrix structure of inputs A and B
    aoclsparse_mat_descr descrB;
    aoclsparse_matrix    csrA;
    aoclsparse_matrix    csrB;

    // Create matrix descriptor of input matrices
    status = aoclsparse_create_mat_descr(&descrB);
    if(status != aoclsparse_status_success)
    {
        std::cerr << "Error returned from aoclsparse_create_mat_descr, status = " << status << "."
                  << std::endl;
        return 1;
    }

    aoclsparse_set_mat_type(descrB, aoclsparse_matrix_type_hermitian);
    aoclsparse_set_mat_fill_mode(descrB, aoclsparse_fill_mode_lower);
    aoclsparse_set_mat_index_base(descrB, aoclsparse_index_base_zero);

    // Matrix sizes
    aoclsparse_int m_a = 5, n_a = 5, m_b = 5, n_b = 5;
    aoclsparse_int nnz_A = 10, nnz_B = 9;
    // Matrix A
    aoclsparse_int            row_ptr_A[] = {0, 1, 2, 5, 9, 10};
    aoclsparse_int            col_ind_A[] = {0, 0, 1, 2, 4, 0, 1, 2, 3, 4};
    aoclsparse_double_complex val_A[]     = {{-0.86, 0.45},
                                             {-2.62, -0.44},
                                             {-0.87, 0.13},
                                             {-0.66, -1.09},
                                             {0.05, -2.37},
                                             {-1.48, -0.42},
                                             {-0.58, -0.70},
                                             {0.31, -0.96},
                                             {-0.88, -2.37},
                                             {-1.23, 0.21}};
    status = aoclsparse_create_zcsr(&csrA, base, m_a, n_a, nnz_A, row_ptr_A, col_ind_A, val_A);
    if(status != aoclsparse_status_success)
    {
        std::cerr << "Error returned from aoclsparse_create_zcsr, status = " << status << "."
                  << std::endl;
        return 2;
    }

    // Hermitian Matrix B
    aoclsparse_int            row_ptr_B[] = {0, 1, 2, 4, 6, 9};
    aoclsparse_int            col_ind_B[] = {0, 1, 0, 2, 1, 3, 1, 2, 4};
    aoclsparse_double_complex val_B[]     = {{-1.59, 0},
                                             {0.46, 0},
                                             {0.07, -0.51},
                                             {-1.52, 0},
                                             {0.21, -1.33},
                                             {-1.37, 0},
                                             {1.42, -2.08},
                                             {-2.26, -1.00},
                                             {-1.81, 0}};

    status = aoclsparse_create_zcsr(&csrB, base, m_b, n_b, nnz_B, row_ptr_B, col_ind_B, val_B);
    if(status != aoclsparse_status_success)
    {
        std::cerr << "Error returned from aoclsparse_create_zcsr, status = " << status << "."
                  << std::endl;
        return 3;
    }

    aoclsparse_matrix          csrC          = NULL;
    aoclsparse_int            *csr_row_ptr_C = NULL;
    aoclsparse_int            *csr_col_ind_C = NULL;
    aoclsparse_double_complex *csr_val_C     = NULL;
    aoclsparse_int             C_M, C_N;

    // expected output matrix
    aoclsparse_int            C_M_exp = 5, C_N_exp = 5, nnz_C_exp = 14;
    aoclsparse_int            csr_row_ptr_C_exp[] = {0, 5, 9, 12, 13, 14};
    aoclsparse_int            csr_col_ind_C_exp[] = {0, 1, 2, 3, 4, 1, 2, 3, 4, 2, 4, 3, 3, 4};
    aoclsparse_double_complex csr_val_C_exp[]     = {{-0.982439, 0.000000},
                                                     {-3.380974, 2.107664},
                                                     {-4.156461, -1.960739},
                                                     {-10.188348, 1.376974},
                                                     {5.609324, 4.536282},
                                                     {-2.308344, 0.000000},
                                                     {-1.331714, -2.631938},
                                                     {-2.972078, -1.039282},
                                                     {-1.922892, -1.975280},
                                                     {-3.862273, 0.000000},
                                                     {-3.714510, 1.465694},
                                                     {-2.743288, 2.163915},
                                                     {-8.756081, -0.000000},
                                                     {-3.022874, 0.000000}};

#ifdef TWO_STAGE_COMPUTATION
    std::cout << "Invoking aoclsparse_sypr with aoclsparse_stage_nnz_count...\n";
    // aoclsparse_stage_nnz_count: Only rowIndex array of the CSR matrix
    // is computed internally.
    request = aoclsparse_stage_nnz_count;
    status  = aoclsparse_sypr(opA, csrA, csrB, descrB, &csrC, request);
    if(status != aoclsparse_status_success)
    {
        std::cerr << "Error returned from aoclsparse_sypr, status = " << status << "." << std::endl;
        return 4;
    }

    std::cout << "Invoking aoclsparse_sypr with aoclsparse_stage_finalize...\n";
    // aoclsparse_stage_finalize: Finalize computation of remaining
    // output arrays (column indices and values of output matrix entries).
    // Has to be called only after aoclsparse_sypr call with
    // aoclsparse_stage_nnz_count parameter.
    request = aoclsparse_stage_finalize;
    status  = aoclsparse_sypr(opA, csrA, csrB, descrB, &csrC, request);
    if(status != aoclsparse_status_success)
    {
        std::cerr << "Error returned from aoclsparse_sypr, status = " << status << "." << std::endl;
        return 5;
    }

#else // SINGLE STAGE
    std::cout << "Invoking aoclsparse_sypr with aoclsparse_stage_full_computation...\n";
    // aoclsparse_stage_full_computation: Whole computation is performed in
    // a single step.
    request = aoclsparse_stage_full_computation;
    status  = aoclsparse_sypr(opA, csrA, csrB, descrB, &csrC, request);
    if(status != aoclsparse_status_success)
    {
        std::cerr << "Error returned from aoclsparse_sypr, status = " << status << "." << std::endl;
        return 6;
    }

#endif

    aoclsparse_export_zcsr(
        csrC, &base, &C_M, &C_N, &nnz_C, &csr_row_ptr_C, &csr_col_ind_C, &csr_val_C);
    // Check and print the result
    std::cout << std::fixed;
    std::cout.precision(1);
    bool oka, okb, okc, oki, okj, okk, ok = true;
    std::cout << std::endl
              << "Output Matrix C: " << std::endl
              << std::setw(11) << "C_M" << std::setw(3) << "" << std::setw(11) << "expected"
              << std::setw(2) << "" << std::setw(11) << "C_N" << std::setw(3) << "" << std::setw(11)
              << "expected" << std::setw(2) << "" << std::setw(11) << "nnz_C" << std::setw(3) << ""
              << std::setw(11) << "expected" << std::endl;
    oka = C_M == C_M_exp;
    ok &= oka;
    std::cout << std::setw(11) << C_M << std::setw(3) << "" << std::setw(11) << C_M_exp
              << std::setw(2) << (oka ? "" : " !");
    okb = C_N == C_N_exp;
    ok &= okb;
    std::cout << std::setw(11) << C_N << std::setw(3) << "" << std::setw(11) << C_N_exp
              << std::setw(2) << (okb ? "" : " !");
    okc = nnz_C == nnz_C_exp;
    ok &= okc;
    std::cout << std::setw(11) << nnz_C << std::setw(3) << "" << std::setw(11) << nnz_C_exp
              << std::setw(2) << (okc ? "" : " !");
    std::cout << std::endl;
    std::cout << std::endl;
    std::cout << std::setw(14) << "csr_val_C" << std::setw(16) << "expected" << std::setw(14)
              << "csr_col_ind_C" << std::setw(14) << "expected" << std::setw(14) << "csr_row_ptr_C"
              << std::setw(14) << "expected" << std::endl;
    // Initializing precision tolerance range for double
    const double tol = 1e-06;
    for(aoclsparse_int i = 0; i < nnz_C; i++)
    {
        oki = ((std::abs(csr_val_C[i].real - csr_val_C_exp[i].real) <= tol)
               && (std::abs(csr_val_C[i].imag - csr_val_C_exp[i].imag) <= tol));
        ok &= oki;
        std::cout << "(" << std::setw(5) << csr_val_C[i].real << "," << std::setw(5)
                  << csr_val_C[i].imag << "i) "
                  << " (" << std::setw(5) << csr_val_C_exp[i].real << "," << std::setw(5)
                  << csr_val_C_exp[i].imag << "i) " << std::setw(2) << (oki ? "" : " !");
        okj = csr_col_ind_C[i] == csr_col_ind_C_exp[i];
        ok &= okj;
        std::cout << std::setw(11) << csr_col_ind_C[i] << std::setw(3) << "" << std::setw(11)
                  << csr_col_ind_C_exp[i] << std::setw(2) << (okj ? "" : " !");
        if(i <= C_M)
        {
            okk = csr_row_ptr_C[i] == csr_row_ptr_C_exp[i];
            ok &= okk;
            std::cout << " " << std::setw(11) << csr_row_ptr_C[i] << std::setw(3) << ""
                      << std::setw(11) << csr_row_ptr_C_exp[i] << std::setw(2) << (okk ? "" : " !");
        }
        std::cout << std::endl;
    }

    aoclsparse_destroy_mat_descr(descrB);
    aoclsparse_destroy(&csrA);
    aoclsparse_destroy(&csrB);
    aoclsparse_destroy(&csrC);
    return (ok ? 0 : 7);
}

Note

aoclsparse_sypr supports only matrices in CSR format which have sorted column indices in each row. If the matrices are unsorted, you might want to call aoclsparse_order_mat().

Note

Currently, opA = aoclsparse_operation_transpose is supported only for real data types.

Parameters

opA – [in] matrix \(A\) operation type.
A – [in] sorted sparse CSR matrix \(A\).
B – [in] sorted sparse CSR matrix \(B\) to be interpreted as symmetric (or Hermitian).
descrB – [in] descriptor of the sparse CSR matrix \(B\). aoclsparse_matrix_type must be aoclsparse_matrix_type_symmetric for real matrices or aoclsparse_matrix_type_hermitian for complex matrices. aoclsparse_fill_mode might be either aoclsparse_fill_mode_upper or aoclsparse_fill_mode_lower to process the upper or lower triangular matrix part, respectively.
request – [in] Specifies if the computation takes place in one stage (aoclsparse_stage_full_computation) or in two stages (aoclsparse_stage_nnz_count followed by aoclsparse_stage_finalize).
*C – [inout] Pointer to the new sparse CSR symmetric/Hermitian matrix \(C\) . Only upper triangle of the result matrix is computed. Matrix \(C\) will always have zero-based indexing, irrespective of the zero/one-based indexing of the input matrices \(A\) and \(B\). The column indices of the output matrix in CSR format might be unsorted. If request is aoclsparse_stage_finalize, matrix \(C\) must not be modified by the user since the last call to aoclsparse_sypr, in the other cases is \(C\) treated as an output only. The matrix should be freed by aoclsparse_destroy() when no longer needed.

Return values

aoclsparse_status_success – the operation completed successfully.
aoclsparse_status_invalid_pointer – descrB, A, B or C is invalid.
aoclsparse_status_invalid_size – Matrix dimensions do not match A or B is not square.
aoclsparse_status_invalid_value – Input parameters are invalid, for example, descrB does not match B indexing or B is not symmetric/Hermitian, C has been modified between stages or opA or request is not recognized.
aoclsparse_status_wrong_type – A and B matrix data types do not match.
aoclsparse_status_not_implemented – Input matrix A or B is not in CSR format.
aoclsparse_status_unsorted_input – Input matrices are not sorted.
aoclsparse_status_memory_error – Memory allocation failure.

aoclsparse_?syprd()#

aoclsparse_status aoclsparse_ssyprd(const aoclsparse_operation op, const aoclsparse_matrix A, const float *B, const aoclsparse_order orderB, const aoclsparse_int ldb, const float alpha, const float beta, float *C, const aoclsparse_order orderC, const aoclsparse_int ldc)#

aoclsparse_status aoclsparse_dsyprd(const aoclsparse_operation op, const aoclsparse_matrix A, const double *B, const aoclsparse_order orderB, const aoclsparse_int ldb, const double alpha, const double beta, double *C, const aoclsparse_order orderC, const aoclsparse_int ldc)#

aoclsparse_status aoclsparse_csyprd(const aoclsparse_operation op, const aoclsparse_matrix A, const aoclsparse_float_complex *B, const aoclsparse_order orderB, const aoclsparse_int ldb, const aoclsparse_float_complex alpha, const aoclsparse_float_complex beta, aoclsparse_float_complex *C, const aoclsparse_order orderC, const aoclsparse_int ldc)#

aoclsparse_status aoclsparse_zsyprd(const aoclsparse_operation op, const aoclsparse_matrix A, const aoclsparse_double_complex *B, const aoclsparse_order orderB, const aoclsparse_int ldb, const aoclsparse_double_complex alpha, const aoclsparse_double_complex beta, aoclsparse_double_complex *C, const aoclsparse_order orderC, const aoclsparse_int ldc)#

Performs symmetric triple product of a sparse matrix and a dense matrix and stores the output as a dense matrix.

aoclsparse_?syprd performs product of a scalar \(\alpha\), with the symmetric triple product of a sparse \(m \times k\) matrix \(A\), defined in CSR format, with a \(k \times k\) symmetric dense (or Hermitian) matrix \(B\), and a \(k \times m\) \(op(A)\). Adds the resulting matrix to \(m \times m\) symmetric dense (or Hermitian) matrix \(C\) that is multiplied by a scalar \(\beta\), such that

\[ C := \alpha \cdot A \cdot B \cdot op(A) + \beta \cdot C \]

if \(op\) is aoclsparse_operation_none.

Otherwise,

\[ C := \alpha \cdot op(A) \cdot B \cdot A + \beta \cdot C \]

\[\begin{split} op(A) = \left\{ \begin{array}{ll} A^T, & \text{if } {\bf\mathsf{op}} = \text{aoclsparse}\_\text{operation}\_\text{transpose} \text{ (real matrices)}\\ A^H, & \text{if } {\bf\mathsf{op}} = \text{aoclsparse}\_\text{operation}\_\text{conjugate}\_\text{transpose} \text{ (complex matrices)} \end{array} \right. \end{split}\]

Notes

1. This routine assumes the dense matrices (B and C) are stored in full although the computations happen on the upper triangular portion of the matrices.

2. aoclsparse_operation_transpose is only supported for real matrices.

3. aoclsparse_operation_conjugate_transpose is only supported for complex matrices.

4. Complex dense matrices are assumed to be Hermitian matrices.

Parameters

op – [in] Matrix \(A\) operation type.
A – [in] Sparse CSR matrix \(A\) structure.
B – [in] Array of dimension \(ldb \times ldb\). Only the upper triangular matrix is used for computation.
orderB – [in] aoclsparse_order_row or aoclsparse_order_column for dense matrix B.
ldb – [in] Leading dimension of \(B\), must be at least \(\max{(1, k)}\) ( \(op(A) = A\)) or \(\max{(1, m)}\) ( \(op(A) = A^T\) or \(op(A) = A^H\)).
alpha – [in] Scalar \(\alpha\).
beta – [in] Scalar \(\beta\).
C – [inout] Array of dimension \(ldc \times ldc\). Only upper triangular part of the matrix is processed.
orderC – [in] aoclsparse_order_row or aoclsparse_order_column for dense matrix C.
ldc – [in] Leading dimension of \(C\), must be at least \(\max{(1, m)}\) ( \(op(A) = A\)) or \(\max{(1, k)}\) ( \(op(A) = A^T\) or \(op(A) = A^H\)).

Return values

aoclsparse_status_success – The operation completed successfully.
aoclsparse_invalid_operation – The operation is invalid if the matrix B and C has a different layout ordering.
aoclsparse_status_wrong_type – The data type of the matrices are not matching or invalid.
aoclsparse_status_invalid_size – The value of m, k, nnz, ldb or ldc is invalid.
aoclsparse_status_invalid_pointer – The pointer A, B, or C is invalid.
aoclsparse_status_not_implemented – The values of orderB and orderC are different.

Miscellaneous#

aoclsparse_ilu_?smoother()#

aoclsparse_status aoclsparse_silu_smoother(aoclsparse_operation op, aoclsparse_matrix A, const aoclsparse_mat_descr descr, float **precond_csr_val, const float *approx_inv_diag, float *x, const float *b)#

aoclsparse_status aoclsparse_dilu_smoother(aoclsparse_operation op, aoclsparse_matrix A, const aoclsparse_mat_descr descr, double **precond_csr_val, const double *approx_inv_diag, double *x, const double *b)#

Incomplete LU factorization with zero fill-in, ILU(0).

Performs incomplete LU factorization with zero fill-in on symmetric sparse matrix A of size \(n \times n\). It also performs a solve for x in

\[ L U x = b, \qquad \text{where} \qquad LU\approx A.\]

Matrix A should be numerically of full rank. Currently single and double precision datatypes are supported.

Example (tests/examples/sample_itsol_d_gmres.cpp)

/* ************************************************************************
 * Copyright (c) 2022-2023 Advanced Micro Devices, Inc.
 *
 * Permission is hereby granted, free of charge, to any person obtaining a copy
 * of this software and associated documentation files (the "Software"), to deal
 * in the Software without restriction, including without limitation the rights
 * to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
 * copies of the Software, and to permit persons to whom the Software is
 * furnished to do so, subject to the following conditions:
 *
 * The above copyright notice and this permission notice shall be included in
 * all copies or substantial portions of the Software.
 *
 * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
 * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
 * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
 * AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
 * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
 * OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
 * THE SOFTWARE.
 *
 * ************************************************************************ */

#include "aoclsparse.h"

#include <cfloat>
#include <iomanip>
#include <iostream>
#include <math.h>
#include <vector>

#define PREMATURE_STOP_TOLERANCE 1e-4
const char gmres_sol[] = "gmres_direct_d";

double calculate_l2Norm_solvers(const double *xSol, const double *x, aoclsparse_int n)
{
    double norm = 0.0f;
    //residual = xSol - x
    //xSol obtained by initial spmv
    //x obtained by ILU-GMRES iterations
    for(aoclsparse_int i = 0; i < n; i++)
    {
        double a = xSol[i] - x[i];
        norm += a * a;
    }
    norm = sqrt(norm);
    return norm;
}

aoclsparse_int monit(aoclsparse_int n,
                     const double  *x,
                     const double  *r __attribute__((unused)),
                     double        *rinfo,
                     void          *udata __attribute__((unused)))
{
    int    it  = (int)rinfo[30];
    double tol = PREMATURE_STOP_TOLERANCE;

    std::ios oldState(nullptr);
    oldState.copyfmt(std::cout);

    std::ios_base::fmtflags fmt = std::cout.flags();
    fmt |= std::ios_base::scientific | std::ios_base::right | std::ios_base::showpos;

    if(!(it % 10))
    {
        std::cout << std::setw(5) << std::right << " iter"
                  << " " << std::setw(16) << std::right << "optim";
        for(int i = 0; i < n; i++)
            std::cout << std::setw(8) << std::right << "x[" << i << "]";
        std::cout << std::endl;
    }
    std::cout << std::setw(5) << std::right << (int)rinfo[30] << " " << std::setw(16) << std::right
              << std::scientific << std::setprecision(8) << rinfo[0];
    std::cout << std::setprecision(2) << std::showpos;
    for(int i = 0; i < n; i++)
        std::cout << " " << x[i];
    std::cout << std::endl;
    std::cout << std::resetiosflags(fmt);

    //reset std::cout state
    std::cout.copyfmt(oldState);

    if(rinfo[0] < tol) // check for premature stop
    {
        return 1; // request to interrupt
    }
    return 0;
}

int main()
{
    std::vector<aoclsparse_int> csr_row_ptr;
    std::vector<aoclsparse_int> csr_col_ind;
    std::vector<double>         csr_val;

    int               m, n, nnz;
    aoclsparse_status exit_status = aoclsparse_status_success;

    std::string filename = "cage4.mtx";
    //symmetry = 0; //real unsymmetric, symmetry = 0%, spd = no
    //https://www.cise.ufl.edu/research/sparse/MM/vanHeukelum/cage4.tar.gz
    n   = 9;
    m   = 9;
    nnz = 49;
    csr_row_ptr.assign({0, 5, 10, 15, 20, 26, 32, 38, 44, 49});
    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,
                        8, 1, 2, 4, 5, 7, 8, 2, 3, 4, 6, 7, 8, 0, 3, 5, 6, 7, 8, 4, 5, 6, 7, 8});
    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,
                    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,
                    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,
                    0.05, 0.14, 0.11, 0.55, 0.25, 0.08, 0.08, 0.09, 0.08, 0.17});

    // Create aocl sparse matrix
    aoclsparse_matrix     A;
    aoclsparse_index_base base = aoclsparse_index_base_zero;
    aoclsparse_mat_descr  descr_a;
    aoclsparse_operation  trans = aoclsparse_operation_none;
    aoclsparse_create_dcsr(&A,
                           base,
                           (aoclsparse_int)n,
                           (aoclsparse_int)n,
                           (aoclsparse_int)nnz,
                           csr_row_ptr.data(),
                           csr_col_ind.data(),
                           csr_val.data());
    aoclsparse_create_mat_descr(&descr_a);
    aoclsparse_set_mat_type(descr_a, aoclsparse_matrix_type_symmetric);
    aoclsparse_set_mat_fill_mode(descr_a, aoclsparse_fill_mode_lower);
    aoclsparse_set_mat_index_base(descr_a, base);
    aoclsparse_set_sv_hint(A, trans, descr_a, 100);
    aoclsparse_optimize(A);

    // Initialize initial point x0 and right hand side b
    double *expected_sol = NULL;
    double *x            = NULL;
    double *b            = NULL;
    double  norm         = 0.0;
    double  rinfo[100];
    double  alpha = 1.0, beta = 0.;
    int     rs_iters = 7;
    char    rs_iters_string[16];

    expected_sol = new double[n];
    x            = new double[n];
    b            = new double[n];

    double init_x = 1.0, ref_x = 0.5;
    for(int i = 0; i < n; i++)
    {
        expected_sol[i] = ref_x;
        x[i]            = init_x;
    }

    aoclsparse_dmv(trans, &alpha, A, descr_a, expected_sol, &beta, b);

    double                  tol    = PREMATURE_STOP_TOLERANCE;
    aoclsparse_itsol_handle handle = NULL;
    // create GMRES handle
    aoclsparse_itsol_d_init(&handle);

    exit_status = aoclsparse_itsol_option_set(handle, "iterative method", "GMRES");
    if(exit_status != aoclsparse_status_success)
        printf("Warning an iterative method option could not be set, exit status = %d\n",
               exit_status);

    exit_status = aoclsparse_itsol_option_set(handle, "gmres preconditioner", "ILU0");
    if(exit_status != aoclsparse_status_success)
        printf("Warning gmres preconditioner option could not be set, exit status = %d\n",
               exit_status);

    exit_status = aoclsparse_itsol_option_set(handle, "gmres iteration limit", "50");
    if(exit_status != aoclsparse_status_success)
        printf("Warning gmres iteration limit option could not be set, exit status = %d\n",
               exit_status);

    sprintf(rs_iters_string, "%d", (int)rs_iters);
    exit_status = aoclsparse_itsol_option_set(handle, "gmres restart iterations", rs_iters_string);
    if(exit_status != aoclsparse_status_success)
        printf("Warning gmres restart iterations option could not be set, exit status = %d\n",
               exit_status);

    aoclsparse_itsol_handle_prn_options(handle);

    // Call GMRES solver
    aoclsparse_status status;

    status = aoclsparse_itsol_d_solve(handle, n, A, descr_a, b, x, rinfo, NULL, monit, &n);
    if(status == aoclsparse_status_success || status == aoclsparse_status_user_stop
       || aoclsparse_status_maxit == status)
    {
        norm = calculate_l2Norm_solvers(expected_sol, x, n);
        printf("input = %s\n", filename.c_str());
        printf("solver = %s\n", gmres_sol);
        printf("no of rows = %d\n", (int)m);
        printf("no of cols = %d\n", (int)n);
        printf("no of nnz = %d\n", (int)nnz);
        printf("monitoring tolerance = %e\n", tol);
        printf("restart iterations = %d\n", (int)rs_iters);
        printf("residual achieved = %e\n", rinfo[0]);
        printf("total iterations = %d\n", (int)rinfo[30]);
        printf("l2 Norm = %e\n", norm);
    }
    else
    {
        std::cout << "Something unexpected happened!, Status = " << status << std::endl;
    }

    delete[] expected_sol;
    delete[] x;
    delete[] b;
    aoclsparse_itsol_destroy(&handle);
    aoclsparse_destroy_mat_descr(descr_a);
    aoclsparse_destroy(&A);
    printf("\n");
    fflush(stdout);

    return 0;
}

Parameters

op – [in] matrix A operation type. Transpose not supported in this release.
A – [in] sparse symmetric matrix handle. Currently ILU functionality is supported only for CSR matrix format.
descr – [in] descriptor of the sparse matrix handle A. Currently, only aoclsparse_matrix_type_symmetric is supported.
precond_csr_val – [out] pointer that contains L and U factors after ILU factorization operation. A is not overwritten with the factors.
approx_inv_diag – [in] Reserved for future use.
x – [out] array of n elements containing the solution to solving aproximatly \(Ax=b\).
b – [in] Right-hand-side of the linear system of equations \(Ax = b\).

Return values

aoclsparse_status_success – the operation completed successfully.
aoclsparse_status_invalid_size – input parameters contain an invalid value.
aoclsparse_status_invalid_pointer – descr, A is invalid.
aoclsparse_status_not_implemented – aoclsparse_matrix_type is not aoclsparse_matrix_type_symmetric or input matrix A is not in CSR format

Sparse BLAS level 1, 2, and 3 functions - Sparse BLAS level 1, 2, and 3 functions - 5.0 English - 63865

AOCL Sparse Library (63865)

Sparse BLAS level 1, 2, and 3 functions#

Level 1#

aoclsparse_?axpyi()#

aoclsparse_?dotci()#

aoclsparse_?dotui()#

aoclsparse_?doti()#

aoclsparse_?sctr()#

sparse_?sctrs()#

aoclsparse_?roti()#

aoclsparse_?gthr()#

aoclsparse_?gthrz()#

aoclsparse_?gthrs()#

Level 2#

aoclsparse_?mv()#

aoclsparse_?trsv()#

aoclsparse_?dotmv()#

aoclsparse_?ellmv()#

aoclsparse_?diamv()#

aoclsparse_?bsrmv()#

aoclsparse_?csrmv()#

aoclsparse_?csrsv()#

Level 3#

aoclsparse_?trsm()#

aoclsparse_sp2m()#

aoclsparse_spmm()#

aoclsparse_?csrmm()#

aoclsparse_?csr2m()#

aoclsparse_?add()#

aoclsparse_?spmmd()#

aoclsparse_?sp2md()#

aoclsparse_syrk()#

aoclsparse_?syrkd()#

aoclsparse_?sypr()#

aoclsparse_?syprd()#

Miscellaneous#

aoclsparse_ilu_?smoother()#