aoclsparse_?sypr() - 5.2 English - 68552

AOCL API Guide (68552)
Document ID
68552
Release Date
2025-12-29
Version
5.2 English
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, \]
or
\[ 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).

  1/* ************************************************************************
  2 * Copyright (c) 2024 Advanced Micro Devices, Inc. All rights reserved.
  3 *
  4 * Permission is hereby granted, free of charge, to any person obtaining a copy
  5 * of this software and associated documentation files (the "Software"), to deal
  6 * in the Software without restriction, including without limitation the rights
  7 * to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
  8 * copies of the Software, and to permit persons to whom the Software is
  9 * furnished to do so, subject to the following conditions:
 10 *
 11 * The above copyright notice and this permission notice shall be included in
 12 * all copies or substantial portions of the Software.
 13 *
 14 * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
 15 * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
 16 * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
 17 * AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
 18 * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
 19 * OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
 20 * THE SOFTWARE.
 21 *
 22 * ************************************************************************ */
 23
 24#include "aoclsparse.h"
 25
 26#include <complex>
 27#include <iomanip>
 28#include <iostream>
 29
 30/* Computes triple symmetric product of complex sparse matrices
 31 * C:=opA(A)*B*A, where opA = conjugate transpose
 32*/
 33
 34// Comment this out for single stage computation
 35//#define TWO_STAGE_COMPUTATION
 36
 37int main(void)
 38{
 39    std::cout << "-------------------------------" << std::endl
 40              << "----- SYPR sample program -----" << std::endl
 41              << "-------------------------------" << std::endl
 42              << std::endl;
 43
 44    aoclsparse_status     status;
 45    aoclsparse_int        nnz_C;
 46    aoclsparse_request    request;
 47    aoclsparse_index_base base = aoclsparse_index_base_zero;
 48    aoclsparse_operation  opA  = aoclsparse_operation_conjugate_transpose;
 49
 50    // Print aoclsparse version
 51    std::cout << aoclsparse_get_version() << std::endl;
 52
 53    // Initialise matrix descriptor and csr matrix structure of inputs A and B
 54    aoclsparse_mat_descr descrB;
 55    aoclsparse_matrix    csrA;
 56    aoclsparse_matrix    csrB;
 57
 58    // Create matrix descriptor of input matrices
 59    status = aoclsparse_create_mat_descr(&descrB);
 60    if(status != aoclsparse_status_success)
 61    {
 62        std::cerr << "Error returned from aoclsparse_create_mat_descr, status = " << status << "."
 63                  << std::endl;
 64        return 1;
 65    }
 66
 67    aoclsparse_set_mat_type(descrB, aoclsparse_matrix_type_hermitian);
 68    aoclsparse_set_mat_fill_mode(descrB, aoclsparse_fill_mode_lower);
 69    aoclsparse_set_mat_index_base(descrB, aoclsparse_index_base_zero);
 70
 71    // Matrix sizes
 72    aoclsparse_int m_a = 5, n_a = 5, m_b = 5, n_b = 5;
 73    aoclsparse_int nnz_A = 10, nnz_B = 9;
 74    // Matrix A
 75    aoclsparse_int            row_ptr_A[] = {0, 1, 2, 5, 9, 10};
 76    aoclsparse_int            col_ind_A[] = {0, 0, 1, 2, 4, 0, 1, 2, 3, 4};
 77    aoclsparse_double_complex val_A[]     = {{-0.86, 0.45},
 78                                             {-2.62, -0.44},
 79                                             {-0.87, 0.13},
 80                                             {-0.66, -1.09},
 81                                             {0.05, -2.37},
 82                                             {-1.48, -0.42},
 83                                             {-0.58, -0.70},
 84                                             {0.31, -0.96},
 85                                             {-0.88, -2.37},
 86                                             {-1.23, 0.21}};
 87    status = aoclsparse_create_zcsr(&csrA, base, m_a, n_a, nnz_A, row_ptr_A, col_ind_A, val_A);
 88    if(status != aoclsparse_status_success)
 89    {
 90        std::cerr << "Error returned from aoclsparse_create_zcsr, status = " << status << "."
 91                  << std::endl;
 92        return 2;
 93    }
 94
 95    // Hermitian Matrix B
 96    aoclsparse_int            row_ptr_B[] = {0, 1, 2, 4, 6, 9};
 97    aoclsparse_int            col_ind_B[] = {0, 1, 0, 2, 1, 3, 1, 2, 4};
 98    aoclsparse_double_complex val_B[]     = {{-1.59, 0},
 99                                             {0.46, 0},
100                                             {0.07, -0.51},
101                                             {-1.52, 0},
102                                             {0.21, -1.33},
103                                             {-1.37, 0},
104                                             {1.42, -2.08},
105                                             {-2.26, -1.00},
106                                             {-1.81, 0}};
107
108    status = aoclsparse_create_zcsr(&csrB, base, m_b, n_b, nnz_B, row_ptr_B, col_ind_B, val_B);
109    if(status != aoclsparse_status_success)
110    {
111        std::cerr << "Error returned from aoclsparse_create_zcsr, status = " << status << "."
112                  << std::endl;
113        return 3;
114    }
115
116    aoclsparse_matrix          csrC          = NULL;
117    aoclsparse_int            *csr_row_ptr_C = NULL;
118    aoclsparse_int            *csr_col_ind_C = NULL;
119    aoclsparse_double_complex *csr_val_C     = NULL;
120    aoclsparse_int             C_M, C_N;
121
122    // expected output matrix
123    aoclsparse_int            C_M_exp = 5, C_N_exp = 5, nnz_C_exp = 14;
124    aoclsparse_int            csr_row_ptr_C_exp[] = {0, 5, 9, 12, 13, 14};
125    aoclsparse_int            csr_col_ind_C_exp[] = {0, 1, 2, 3, 4, 1, 2, 3, 4, 2, 4, 3, 3, 4};
126    aoclsparse_double_complex csr_val_C_exp[]     = {{-0.982439, 0.000000},
127                                                     {-3.380974, 2.107664},
128                                                     {-4.156461, -1.960739},
129                                                     {-10.188348, 1.376974},
130                                                     {5.609324, 4.536282},
131                                                     {-2.308344, 0.000000},
132                                                     {-1.331714, -2.631938},
133                                                     {-2.972078, -1.039282},
134                                                     {-1.922892, -1.975280},
135                                                     {-3.862273, 0.000000},
136                                                     {-3.714510, 1.465694},
137                                                     {-2.743288, 2.163915},
138                                                     {-8.756081, -0.000000},
139                                                     {-3.022874, 0.000000}};
140
141#ifdef TWO_STAGE_COMPUTATION
142    std::cout << "Invoking aoclsparse_sypr with aoclsparse_stage_nnz_count...\n";
143    // aoclsparse_stage_nnz_count: Only rowIndex array of the CSR matrix
144    // is computed internally.
145    request = aoclsparse_stage_nnz_count;
146    status  = aoclsparse_sypr(opA, csrA, csrB, descrB, &csrC, request);
147    if(status != aoclsparse_status_success)
148    {
149        std::cerr << "Error returned from aoclsparse_sypr, status = " << status << "." << std::endl;
150        return 4;
151    }
152
153    std::cout << "Invoking aoclsparse_sypr with aoclsparse_stage_finalize...\n";
154    // aoclsparse_stage_finalize: Finalize computation of remaining
155    // output arrays (column indices and values of output matrix entries).
156    // Has to be called only after aoclsparse_sypr call with
157    // aoclsparse_stage_nnz_count parameter.
158    request = aoclsparse_stage_finalize;
159    status  = aoclsparse_sypr(opA, csrA, csrB, descrB, &csrC, request);
160    if(status != aoclsparse_status_success)
161    {
162        std::cerr << "Error returned from aoclsparse_sypr, status = " << status << "." << std::endl;
163        return 5;
164    }
165
166#else // SINGLE STAGE
167    std::cout << "Invoking aoclsparse_sypr with aoclsparse_stage_full_computation...\n";
168    // aoclsparse_stage_full_computation: Whole computation is performed in
169    // a single step.
170    request = aoclsparse_stage_full_computation;
171    status  = aoclsparse_sypr(opA, csrA, csrB, descrB, &csrC, request);
172    if(status != aoclsparse_status_success)
173    {
174        std::cerr << "Error returned from aoclsparse_sypr, status = " << status << "." << std::endl;
175        return 6;
176    }
177
178#endif
179
180    aoclsparse_export_zcsr(
181        csrC, &base, &C_M, &C_N, &nnz_C, &csr_row_ptr_C, &csr_col_ind_C, &csr_val_C);
182    // Check and print the result
183    std::cout << std::fixed;
184    std::cout.precision(1);
185    bool oka, okb, okc, oki, okj, okk, ok = true;
186    std::cout << std::endl
187              << "Output Matrix C: " << std::endl
188              << std::setw(11) << "C_M" << std::setw(3) << "" << std::setw(11) << "expected"
189              << std::setw(2) << "" << std::setw(11) << "C_N" << std::setw(3) << "" << std::setw(11)
190              << "expected" << std::setw(2) << "" << std::setw(11) << "nnz_C" << std::setw(3) << ""
191              << std::setw(11) << "expected" << std::endl;
192    oka = C_M == C_M_exp;
193    ok &= oka;
194    std::cout << std::setw(11) << C_M << std::setw(3) << "" << std::setw(11) << C_M_exp
195              << std::setw(2) << (oka ? "" : " !");
196    okb = C_N == C_N_exp;
197    ok &= okb;
198    std::cout << std::setw(11) << C_N << std::setw(3) << "" << std::setw(11) << C_N_exp
199              << std::setw(2) << (okb ? "" : " !");
200    okc = nnz_C == nnz_C_exp;
201    ok &= okc;
202    std::cout << std::setw(11) << nnz_C << std::setw(3) << "" << std::setw(11) << nnz_C_exp
203              << std::setw(2) << (okc ? "" : " !");
204    std::cout << std::endl;
205    std::cout << std::endl;
206    std::cout << std::setw(14) << "csr_val_C" << std::setw(16) << "expected" << std::setw(14)
207              << "csr_col_ind_C" << std::setw(14) << "expected" << std::setw(14) << "csr_row_ptr_C"
208              << std::setw(14) << "expected" << std::endl;
209    // Initializing precision tolerance range for double
210    const double tol = 1e-06;
211    for(aoclsparse_int i = 0; i < nnz_C; i++)
212    {
213        oki = ((std::abs(csr_val_C[i].real - csr_val_C_exp[i].real) <= tol)
214               && (std::abs(csr_val_C[i].imag - csr_val_C_exp[i].imag) <= tol));
215        ok &= oki;
216        std::cout << "(" << std::setw(5) << csr_val_C[i].real << "," << std::setw(5)
217                  << csr_val_C[i].imag << "i) "
218                  << " (" << std::setw(5) << csr_val_C_exp[i].real << "," << std::setw(5)
219                  << csr_val_C_exp[i].imag << "i) " << std::setw(2) << (oki ? "" : " !");
220        okj = csr_col_ind_C[i] == csr_col_ind_C_exp[i];
221        ok &= okj;
222        std::cout << std::setw(11) << csr_col_ind_C[i] << std::setw(3) << "" << std::setw(11)
223                  << csr_col_ind_C_exp[i] << std::setw(2) << (okj ? "" : " !");
224        if(i <= C_M)
225        {
226            okk = csr_row_ptr_C[i] == csr_row_ptr_C_exp[i];
227            ok &= okk;
228            std::cout << " " << std::setw(11) << csr_row_ptr_C[i] << std::setw(3) << ""
229                      << std::setw(11) << csr_row_ptr_C_exp[i] << std::setw(2) << (okk ? "" : " !");
230        }
231        std::cout << std::endl;
232    }
233
234    aoclsparse_destroy_mat_descr(descrB);
235    aoclsparse_destroy(&csrA);
236    aoclsparse_destroy(&csrB);
237    aoclsparse_destroy(&csrC);
238    return (ok ? 0 : 7);
239}

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