/*! \file */ /* ************************************************************************ * Copyright (C) 2022-2026 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 "rocsparse_enum.hpp" #include "testing.hpp" // // Split such that L is unit and U is non-unit. // template rocsparse_status rocsparse_host_csritsv_ptr_end(rocsparse_fill_mode fill_mode_, rocsparse_diag_type diag_type_, J m_, I nnz_, const I* __restrict__ ptr_, I* __restrict__ ptr_end_, const J* __restrict__ ind_, rocsparse_index_base base_, rocsparse_int* zero_pivot) { zero_pivot[0] = -1; switch(fill_mode_) { case rocsparse_fill_mode_lower: { switch(diag_type_) { case rocsparse_diag_type_unit: { for(J i = 0; i < m_; ++i) { ptr_end_[i] = ptr_[i + 1]; for(I k = ptr_[i] - base_; k < ptr_[i + 1] - base_; ++k) { const J j = ind_[k] - base_; if(j >= i) { ptr_end_[i] = k + base_; break; } } } break; } case rocsparse_diag_type_non_unit: { zero_pivot[0] = std::numeric_limits::max(); J count_symbolic_diag = 0; for(J i = 0; i < m_; ++i) { ptr_end_[i] = ptr_[i + 1]; bool mark = false; for(I k = ptr_[i] - base_; k < ptr_[i + 1] - base_; ++k) { const J j = ind_[k] - base_; if(j == i) { mark = true; ptr_end_[i] = k + 1 + base_; break; } } if(!mark) { zero_pivot[0] = std::min(zero_pivot[0], i + base_); ++count_symbolic_diag; } } if(zero_pivot[0] == std::numeric_limits::max()) { zero_pivot[0] = -1; } if(count_symbolic_diag > 0) { return rocsparse_status_success; } break; } } break; } case rocsparse_fill_mode_upper: { switch(diag_type_) { case rocsparse_diag_type_unit: { for(J i = 0; i < m_; ++i) { ptr_end_[i] = ptr_[i + 1]; for(I k = ptr_[i] - base_; k < ptr_[i + 1] - base_; ++k) { const J j = ind_[k] - base_; if(j > i) { ptr_end_[i] = k + base_; break; } } } break; } case rocsparse_diag_type_non_unit: { zero_pivot[0] = std::numeric_limits::max(); J count_symbolic_diag = 0; for(J i = 0; i < m_; ++i) { bool mark = false; ptr_end_[i] = ptr_[i + 1]; for(I k = ptr_[i] - base_; k < ptr_[i + 1] - base_; ++k) { const J j = ind_[k] - base_; if(j == i) { ptr_end_[i] = k + base_; mark = true; break; } } if(!mark) { zero_pivot[0] = std::min(zero_pivot[0], i + base_); ; ++count_symbolic_diag; } } if(zero_pivot[0] == std::numeric_limits::max()) { zero_pivot[0] = -1; } if(count_symbolic_diag > 0) { return rocsparse_status_success; } break; } } break; } } return rocsparse_status_success; } template rocsparse_status rocsparse_host_csritsv_buffer_size_template(rocsparse_handle handle, rocsparse_operation trans, J m, I nnz, const rocsparse_mat_descr descr, const T* csr_val, const I* csr_row_ptr, const J* csr_col_ind, rocsparse_mat_info info, size_t* buffer_size) { buffer_size[0] = ((sizeof(I) * m - 1) / 256 + 1) * 256 + sizeof(T) * m; const rocsparse_diag_type diag_type = rocsparse_get_mat_diag_type(descr); if(diag_type == rocsparse_diag_type_non_unit) { buffer_size[0] += sizeof(T) * m; } return rocsparse_status_success; } // // Non optimal host implementation. // template rocsparse_status rocsparse_host_csritsv_solve(rocsparse_handle handle, rocsparse_int* host_nmaxiter, const floating_data_t* host_tol, floating_data_t* host_history, rocsparse_operation trans, J m, I nnz, const T* alpha, const rocsparse_mat_descr descr, const T* csr_val, const I* csr_row_ptr, const J* csr_col_ind, const T* x, T* y, void* temp_buffer, rocsparse_int* zero_pivot) { zero_pivot[0] = -1; static constexpr bool verbose = false; const rocsparse_fill_mode fill_mode = rocsparse_get_mat_fill_mode(descr); const rocsparse_diag_type diag_type = rocsparse_get_mat_diag_type(descr); const rocsparse_matrix_type mat_type = rocsparse_get_mat_type(descr); if(verbose) { std::cout << "diag_type_" << diag_type << std::endl; std::cout << "fill_mode_" << fill_mode << std::endl; } if(m == 0 || nnz == 0) { if(nnz == 0 && diag_type == rocsparse_diag_type_unit) { // // copy and scal. // for(J i = 0; i < m; ++i) y[i] = alpha[0] * x[i]; host_nmaxiter[0] = 1; } return rocsparse_status_success; } const I* ptr_end = nullptr; T* y_p = nullptr; T* inv_diag = nullptr; rocsparse_index_base base = rocsparse_get_mat_index_base(descr); if(mat_type == rocsparse_matrix_type_general) { ptr_end = (I*)temp_buffer; rocsparse_status status = rocsparse_host_csritsv_ptr_end(fill_mode, diag_type, m, nnz, csr_row_ptr, (I*)temp_buffer, csr_col_ind, base, zero_pivot); if(status != rocsparse_status_success) { return status; } temp_buffer = (void*)(((char*)temp_buffer) + ((sizeof(I) * m - 1) / 256 + 1) * 256); y_p = (T*)temp_buffer; temp_buffer = (void*)(y_p + m); } else if(mat_type == rocsparse_matrix_type_triangular) { y_p = (T*)temp_buffer; temp_buffer = (void*)(y_p + m); if(fill_mode == rocsparse_fill_mode_lower) { ptr_end = csr_row_ptr + 1; } else { ptr_end = csr_row_ptr; } switch(diag_type) { case rocsparse_diag_type_non_unit: { if(fill_mode == rocsparse_fill_mode_lower) { for(J i = 0; i < m; ++i) { const J j = csr_col_ind[csr_row_ptr[i + 1] - base - 1] - base; if(i != j) { zero_pivot[0] = i + base; break; } } } else { for(J i = 0; i < m; ++i) { const J j = csr_col_ind[csr_row_ptr[i] - base] - base; if(i != j) { zero_pivot[0] = i + base; break; } } } break; } case rocsparse_diag_type_unit: { break; } } } if(zero_pivot[0] != -1) { return rocsparse_status_success; } const I* b = nullptr; const I* e = nullptr; const I* d = nullptr; switch(fill_mode) { case rocsparse_fill_mode_lower: { b = csr_row_ptr; e = ptr_end; break; } case rocsparse_fill_mode_upper: { b = ptr_end; e = csr_row_ptr + 1; break; } } switch(diag_type) { case rocsparse_diag_type_non_unit: { d = ptr_end; break; } case rocsparse_diag_type_unit: { break; } } switch(diag_type) { case rocsparse_diag_type_unit: { break; } case rocsparse_diag_type_non_unit: { inv_diag = (T*)temp_buffer; temp_buffer = (void*)(inv_diag + m); for(J i = 0; i < m; ++i) { I k = (fill_mode == rocsparse_fill_mode_upper) ? (d[i] - base) : (d[i] - base - 1); if(csr_val[k] == static_cast(0)) { zero_pivot[0] = i + base; return rocsparse_status_success; } if(trans == rocsparse_operation_conjugate_transpose) { inv_diag[i] = static_cast(1) / rocsparse_conj(csr_val[k]); } else { inv_diag[i] = static_cast(1) / csr_val[k]; } } break; } } const floating_data_t nrm0 = static_cast>(1); // // Iterative Loop. // for(J iter = 0; iter < host_nmaxiter[0]; ++iter) { // // Copy y to y_p. // for(J i = 0; i < m; ++i) { y_p[i] = y[i]; } floating_data_t mx_residual = static_cast>(0); floating_data_t mx = static_cast>(0); // // Compute y = alpha // switch(trans) { case rocsparse_operation_none: { switch(diag_type) { case rocsparse_diag_type_non_unit: { for(J i = 0; i < m; ++i) { T sum = static_cast(0); if((e[i] > b[i] + 1)) { for(I k = b[i] - base; k < e[i] - base; ++k) { sum += csr_val[k] * y_p[csr_col_ind[k] - base]; } const T h = inv_diag[i] * (alpha[0] * x[i] - sum); mx = std::max(mx, std::abs(h)); mx_residual = std::max(mx_residual, std::abs(alpha[0] * x[i] - sum)); y[i] = y_p[i] + h; } else { y[i] = inv_diag[i] * alpha[0] * x[i]; mx = std::max(mx, std::abs(y[i] - y_p[i])); mx_residual = std::max(mx_residual, std::abs(alpha[0] * x[i] - y_p[i] / inv_diag[i])); } } break; } case rocsparse_diag_type_unit: { for(J i = 0; i < m; ++i) { T sum = static_cast(0); for(I k = b[i] - base; k < e[i] - base; ++k) sum += csr_val[k] * y_p[csr_col_ind[k] - base]; y[i] = alpha[0] * x[i] - sum; const T h = y[i] - y_p[i]; mx = std::max(mx, std::abs(h)); mx_residual = mx; } break; } } break; } case rocsparse_operation_transpose: case rocsparse_operation_conjugate_transpose: { for(J i = 0; i < m; ++i) { y[i] = static_cast(0); } for(J i = 0; i < m; ++i) { // row i, column csr_col_ind[k] // row csr_col_ind[k] for(I k = b[i] - base; k < e[i] - base; ++k) { const J j = csr_col_ind[k] - base; const T a = (trans == rocsparse_operation_conjugate_transpose) ? rocsparse_conj(csr_val[k]) : csr_val[k]; y[j] += a * y_p[i]; } } switch(diag_type) { case rocsparse_diag_type_non_unit: { for(J i = 0; i < m; ++i) { mx_residual = std::max(mx, std::abs(alpha[0] * x[i] - y[i])); const T h = inv_diag[i] * (alpha[0] * x[i] - y[i]); mx = std::max(mx, std::abs(h)); y[i] = h + y_p[i]; } break; } case rocsparse_diag_type_unit: { for(J i = 0; i < m; ++i) { y[i] = (alpha[0] * x[i] - y[i]); const T h = y[i] - y_p[i]; mx = std::max(mx, std::abs(h)); mx_residual = mx; } break; } } break; } } // // y_k+1 = yk + (alpha * x - (id + T) * yk ) // if(verbose) { std::cout << "iter " << iter << ", mx " << mx / nrm0 << ", mx_residual " << mx_residual / nrm0 << std::endl; } if(host_history) { host_history[iter] = mx; } if(host_tol && (mx_residual <= host_tol[0])) { host_nmaxiter[0] = iter + 1; break; } } return rocsparse_status_success; } template void testing_csritsv_bad_arg(const Arguments& arg) { static constexpr bool verbose = false; // Create rocsparse handle rocsparse_local_handle local_handle; // Create matrix descriptor rocsparse_local_mat_descr local_descr; // Create matrix info rocsparse_local_mat_info local_info; const T h_alpha = static_cast(1); rocsparse_handle handle = local_handle; rocsparse_int* host_nmaxiter = (rocsparse_int*)0x4; rocsparse_int host_nfreeiter = 2; const floating_data_t* host_tol = (const floating_data_t*)0x4; floating_data_t* host_history = (floating_data_t*)0x4; rocsparse_operation trans = rocsparse_operation_none; rocsparse_int m = 32; rocsparse_int nnz = 32; const T* alpha_device_host = &h_alpha; const rocsparse_mat_descr descr = local_descr; const T* csr_val = (const T*)0x4; const rocsparse_int* csr_row_ptr = (const rocsparse_int*)0x4; const rocsparse_int* csr_col_ind = (const rocsparse_int*)0x4; rocsparse_mat_info info = local_info; const T* x = (const T*)0x4; T* y = (T*)0x4; rocsparse_solve_policy solve = rocsparse_solve_policy_auto; rocsparse_solve_policy policy = rocsparse_solve_policy_auto; rocsparse_analysis_policy analysis = rocsparse_analysis_policy_force; void* temp_buffer = (void*)0x4; size_t* buffer_size = (size_t*)0x4; #define PARAMS_BUFFER_SIZE \ handle, trans, m, nnz, descr, csr_val, csr_row_ptr, csr_col_ind, info, buffer_size #define PARAMS_ANALYSIS \ handle, trans, m, nnz, descr, csr_val, csr_row_ptr, csr_col_ind, info, analysis, solve, \ temp_buffer #define PARAMS_SOLVE \ handle, host_nmaxiter, host_tol, host_history, trans, m, nnz, alpha_device_host, descr, \ csr_val, csr_row_ptr, csr_col_ind, info, x, y, policy, temp_buffer #define PARAMSX_SOLVE \ handle, host_nmaxiter, host_nfreeiter, host_tol, host_history, trans, m, nnz, \ alpha_device_host, descr, csr_val, csr_row_ptr, csr_col_ind, info, x, y, policy, \ temp_buffer // // // if(verbose) { std::cout << "bad_arg_analysis(rocsparse_csritsv_buffer_size, PARAMS_BUFFER_SIZE)" << std::endl; } bad_arg_analysis(rocsparse_csritsv_buffer_size, PARAMS_BUFFER_SIZE); // // // if(verbose) { std::cout << "bad_arg_analysis(rocsparse_csritsv_analysis, PARAMS_ANALYSIS)" << std::endl; } bad_arg_analysis(rocsparse_csritsv_analysis, PARAMS_ANALYSIS); // // // if(verbose) { std::cout << "bad_arg_analysis(rocsparse_csritsv_solve, PARAMS_SOLVE)" << std::endl; } // // // { static constexpr int nargs_to_exclude = 2; static constexpr int args_to_exclude[nargs_to_exclude] = {2, 3}; select_bad_arg_analysis( rocsparse_csritsv_solve, nargs_to_exclude, args_to_exclude, PARAMS_SOLVE); } // // // { static constexpr int nargs_to_exclude = 2; static constexpr int args_to_exclude[nargs_to_exclude] = {3, 4}; select_bad_arg_analysis( rocsparse_csritsv_solve_ex, nargs_to_exclude, args_to_exclude, PARAMSX_SOLVE); } for(auto matrix_type : rocsparse_matrix_type_t::values) { if(matrix_type != rocsparse_matrix_type_general && matrix_type != rocsparse_matrix_type_triangular) { CHECK_ROCSPARSE_ERROR(rocsparse_set_mat_type(descr, matrix_type)); EXPECT_ROCSPARSE_STATUS(rocsparse_csritsv_buffer_size(PARAMS_BUFFER_SIZE), rocsparse_status_not_implemented); EXPECT_ROCSPARSE_STATUS(rocsparse_csritsv_analysis(PARAMS_ANALYSIS), rocsparse_status_not_implemented); EXPECT_ROCSPARSE_STATUS(rocsparse_csritsv_solve(PARAMS_SOLVE), rocsparse_status_not_implemented); EXPECT_ROCSPARSE_STATUS(rocsparse_csritsv_solve_ex(PARAMSX_SOLVE), rocsparse_status_not_implemented); } } CHECK_ROCSPARSE_ERROR(rocsparse_set_mat_type(descr, rocsparse_matrix_type_general)); CHECK_ROCSPARSE_ERROR(rocsparse_set_mat_storage_mode(descr, rocsparse_storage_mode_unsorted)); EXPECT_ROCSPARSE_STATUS(rocsparse_csritsv_buffer_size(PARAMS_BUFFER_SIZE), rocsparse_status_requires_sorted_storage); EXPECT_ROCSPARSE_STATUS(rocsparse_csritsv_analysis(PARAMS_ANALYSIS), rocsparse_status_requires_sorted_storage); EXPECT_ROCSPARSE_STATUS(rocsparse_csritsv_solve(PARAMS_SOLVE), rocsparse_status_requires_sorted_storage); EXPECT_ROCSPARSE_STATUS(rocsparse_csritsv_solve_ex(PARAMSX_SOLVE), rocsparse_status_requires_sorted_storage); CHECK_ROCSPARSE_ERROR(rocsparse_set_mat_storage_mode(descr, rocsparse_storage_mode_sorted)); #undef PARAMS_BUFFER_SIZE #undef PARAMS_ANALYSIS #undef PARAMS_SOLVE #undef PARAMSX_SOLVE // Test rocsparse_csritsv_zero_pivot() rocsparse_int position; EXPECT_ROCSPARSE_STATUS(rocsparse_csritsv_zero_pivot(nullptr, descr, info, &position), rocsparse_status_invalid_handle); EXPECT_ROCSPARSE_STATUS(rocsparse_csritsv_zero_pivot(handle, descr, nullptr, &position), rocsparse_status_invalid_pointer); EXPECT_ROCSPARSE_STATUS(rocsparse_csritsv_zero_pivot(handle, descr, info, nullptr), rocsparse_status_invalid_pointer); // Test rocsparse_csritsv_clear() EXPECT_ROCSPARSE_STATUS(rocsparse_csritsv_clear(nullptr, descr, info), rocsparse_status_invalid_handle); EXPECT_ROCSPARSE_STATUS(rocsparse_csritsv_clear(handle, nullptr, info), rocsparse_status_invalid_pointer); EXPECT_ROCSPARSE_STATUS(rocsparse_csritsv_clear(handle, descr, nullptr), rocsparse_status_invalid_pointer); } template void testing_csritsv(const Arguments& arg) { static const bool verbose = false; // // Set nmaxiter. // rocsparse_int host_nmaxiter[1] = {arg.nmaxiter}; rocsparse_int host_nfreeiter = arg.nfreeiter; // // Tolerance for the iterative method. // floating_data_t tol_iterative = static_cast>(1.0e-6); if(std::is_same, double>{}) tol_iterative = static_cast>(1.0e-14); floating_data_t host_tol[1] = {tol_iterative}; // // Tolerance for the comparison. // floating_data_t tol_compare = static_cast>(1.0e-5); if(std::is_same, double>{}) tol_compare = static_cast>(1.0e-10); floating_data_t tol_compare_with_direct = static_cast>(1.0e-3); if(std::is_same, double>{}) tol_compare_with_direct = static_cast>(1.0e-10); // // Convergence history. // host_vector> host_history(host_nmaxiter[0]); // // Grab parametrization. // rocsparse_int M = arg.M; rocsparse_int N = arg.N; rocsparse_operation trans = arg.transA; rocsparse_diag_type diag = arg.diag; rocsparse_fill_mode uplo = arg.uplo; rocsparse_analysis_policy apol = arg.apol; rocsparse_solve_policy spol = arg.spol; rocsparse_index_base base = arg.baseA; // // Create host alpha. // host_scalar h_alpha(arg.get_alpha()); // // Create rocsparse handle // rocsparse_local_handle handle; // // Create matrix info. // rocsparse_local_mat_info info; // // Create matrix descriptor // rocsparse_local_mat_descr descr; // // Configure descriptor. // CHECK_ROCSPARSE_ERROR(rocsparse_set_mat_type(descr, arg.matrix_type)); CHECK_ROCSPARSE_ERROR(rocsparse_set_mat_diag_type(descr, diag)); CHECK_ROCSPARSE_ERROR(rocsparse_set_mat_fill_mode(descr, uplo)); CHECK_ROCSPARSE_ERROR(rocsparse_set_mat_index_base(descr, base)); // // Define routine parameters. // #define PARAMS_BUFFER_SIZE(A_) \ handle, trans, A_.m, A_.nnz, descr, A_.val, A_.ptr, A_.ind, info, &buffer_size #define PARAMS_ANALYSIS(A_) \ handle, trans, A_.m, A_.nnz, descr, A_.val, A_.ptr, A_.ind, info, apol, spol, dbuffer #define PARAMS_SOLVE(alpha_, A_, x_, y_) \ handle, host_nmaxiter, host_tol, host_history, trans, A_.m, A_.nnz, alpha_, descr, A_.val, \ A_.ptr, A_.ind, info, x_, y_, spol, dbuffer #define PARAMSX_SOLVE(alpha_, A_, x_, y_) \ handle, host_nmaxiter, host_nfreeiter, host_tol, host_history, trans, A_.m, A_.nnz, alpha_, \ descr, A_.val, A_.ptr, A_.ind, info, x_, y_, spol, dbuffer rocsparse_int host_zero_pivot; #define HOST_PARAMS_SOLVE(alpha_, A_, x_, y_) \ handle, host_nmaxiter, host_tol, host_history, trans, A_.m, A_.nnz, alpha_, descr, A_.val, \ A_.ptr, A_.ind, x_, y_, hbuffer, &host_zero_pivot // // Non-squared matrices are not supported // if(M != N) { if(verbose) { std::cerr << "// rocSPARSE.WARNING clients testing_csritsv, skipping non-squared matrices." << std::endl; } return; } if(M == 0) { size_t buffer_size; rocsparse_int pivot; device_vector dx, dy, dbuffer; device_csr_matrix dA; CHECK_ROCSPARSE_ERROR(rocsparse_set_pointer_mode(handle, rocsparse_pointer_mode_host)); EXPECT_ROCSPARSE_STATUS(rocsparse_csritsv_buffer_size(PARAMS_BUFFER_SIZE(dA)), rocsparse_status_success); EXPECT_ROCSPARSE_STATUS(rocsparse_csritsv_analysis(PARAMS_ANALYSIS(dA)), rocsparse_status_success); EXPECT_ROCSPARSE_STATUS(rocsparse_csritsv_solve(PARAMS_SOLVE(h_alpha, dA, dx, dy)), rocsparse_status_success); EXPECT_ROCSPARSE_STATUS(rocsparse_csritsv_solve_ex(PARAMSX_SOLVE(h_alpha, dA, dx, dy)), rocsparse_status_success); EXPECT_ROCSPARSE_STATUS(rocsparse_csritsv_zero_pivot(handle, descr, info, &pivot), rocsparse_status_success); EXPECT_ROCSPARSE_STATUS(rocsparse_csritsv_clear(handle, descr, info), rocsparse_status_success); return; } // // Declare matrix A. // host_csr_matrix hA; if(diag == rocsparse_diag_type_unit && rocsparse_matrix_type_triangular == arg.matrix_type) { host_csr_matrix hB; static constexpr bool to_int = false; static constexpr bool full_rank = true; rocsparse_matrix_factory matrix_factory(arg, to_int, full_rank); matrix_factory.init_csr(hB, M, N); // // Let's remove the diagonal elements... // int ndiag = 0; for(int i = 0; i < hB.m; ++i) { for(int k = hB.ptr[i] - hB.base; k < hB.ptr[i + 1] - hB.base; ++k) { const int j = hB.ind[k] - hB.base; if(i == j) { ++ndiag; } } } hA.define(hB.m, hB.n, hB.nnz - ndiag, hB.base); if(ndiag > 0) { hA.ptr[0] = hA.base; for(int i = 0; i < hB.m; ++i) { bool find = false; for(int k = hB.ptr[i] - hB.base; k < hB.ptr[i + 1] - hB.base; ++k) { const int j = hB.ind[k] - hB.base; if(i == j) { find = true; break; } } if(find) { hA.ptr[i + 1] = hB.ptr[i + 1] - hB.ptr[i] - 1; } else { hA.ptr[i + 1] = hB.ptr[i + 1] - hB.ptr[i]; } } for(int i = 1; i <= hA.m; ++i) hA.ptr[i] += hA.ptr[i - 1]; for(int i = 0; i < hB.m; ++i) { int s = 0; for(int k = hB.ptr[i] - hB.base; k < hB.ptr[i + 1] - hB.base; ++k) { if((hB.ind[k] - hB.base) != i) { hA.ind[hA.ptr[i] - hA.base + s] = hB.ind[k] - hB.base + hA.base; hA.val[hA.ptr[i] - hA.base + s] = hB.val[k]; ++s; } } } } else { hA.transfer_from(hB); } } else { static constexpr bool to_int = false; static constexpr bool full_rank = true; rocsparse_matrix_factory matrix_factory(arg, to_int, full_rank); matrix_factory.init_csr(hA, M, N); } // // Again, since we import or generate the matrix, Non-squared matrices are not supported. // if(M != N) { if(verbose) { std::cerr << "// rocSPARSE.WARNING clients testing_csritsv, skipping non-squared matrices." << std::endl; } return; } // // Let's normalize ... it helps on horrible matrices. // floating_data_t mx = 0; for(rocsparse_int i = 0; i < hA.nnz; ++i) mx = std::max(mx, std::abs(hA.val[i])); if(mx > 0) { for(rocsparse_int i = 0; i < hA.nnz; ++i) hA.val[i] /= mx; } // // Declare and initialize host X. // host_dense_matrix hx(M, 1); rocsparse_matrix_utils::init(hx); // // Define and transfer A and X from host to device. // device_csr_matrix dA(hA); device_dense_matrix dx(hx); // // Define Y on device. // device_dense_matrix dy(M, 1); host_scalar h_analysis_pivot, h_solve_pivot; if(verbose) { std::cout << "M : " << M << std::endl; std::cout << "N : " << N << std::endl; std::cout << "NNZ : " << hA.nnz << std::endl; } // // Buffer for calculation on device. // void* dbuffer; { size_t buffer_size; CHECK_ROCSPARSE_ERROR(rocsparse_csritsv_buffer_size(PARAMS_BUFFER_SIZE(dA))); if(verbose) { std::cout << "csritsv_buffer_size " << buffer_size << std::endl; } CHECK_HIP_ERROR(rocsparse_hipMalloc(&dbuffer, buffer_size)); } if(arg.unit_check) { if(verbose) { std::cout << "testing_csritsv unit_check pointer mode host ... " << std::endl; } if(verbose) { std::cout << " - compute host iterative" << std::endl; } // // Buffer for calculation on host. // void* hbuffer; { size_t buffer_size; CHECK_ROCSPARSE_ERROR( (rocsparse_host_csritsv_buffer_size_template)(PARAMS_BUFFER_SIZE(hA))); if(verbose) { std::cout << "host_csritsv_buffer_size " << buffer_size << std::endl; } hbuffer = malloc(buffer_size); } // // Calculate on host. // host_dense_matrix hy_iterative(M, 1); for(rocsparse_int i = 0; i < M; ++i) { hy_iterative[i] = static_cast(0); } // // Compute the iterative method on host. // CHECK_ROCSPARSE_ERROR( (rocsparse_host_csritsv_solve)(HOST_PARAMS_SOLVE( h_alpha, hA, hx, hy_iterative))); const bool host_iterative_convergence = (host_nmaxiter[0] < arg.nmaxiter); if(false == host_iterative_convergence) { // // it didn't converge, for some reasons ... // if(verbose) { std::cerr << "host csritsv didn't converge. " << std::endl; } } // // Calculate the direct method on host. // if(verbose) { std::cout << " - compute host direct" << std::endl; } host_dense_matrix hy_direct(M, 1); host_csrsv(trans, hA.m, hA.nnz, *h_alpha, hA.ptr, hA.ind, hA.val, hx, (int64_t)1, hy_direct, diag, uplo, base, h_analysis_pivot, h_solve_pivot); if(verbose) { std::cout << "h_analysis_pivot " << *h_analysis_pivot << std::endl; std::cout << "h_solve_pivot " << *h_solve_pivot << std::endl; } if(*h_analysis_pivot == -1 && *h_solve_pivot == -1) { if(host_iterative_convergence) { if(verbose) { std::cout << " - compare host iterative and direct solutions" << std::endl; } // // Compare direct solution and iterative solution calculated on host // In other words, make sure what we calculate with the iterative method (even on host) makes sense. // hy_direct.near_check(hy_iterative, tol_compare_with_direct); } } else { if(verbose) { std::cout << " - zero pivot detected with direct method" << std::endl; } } if(verbose) { std::cout << " - compute device iterative" << std::endl; } // // Now we calculate on gpu // host_scalar analysis_no_pivot(-1); host_scalar analysis_pivot; host_scalar solve_pivot; CHECK_ROCSPARSE_ERROR(rocsparse_set_pointer_mode(handle, rocsparse_pointer_mode_host)); if(verbose) { std::cout << " - compute device iterative, check if default zero pivot is -1" << std::endl; } // // CHECK IF DEFAULT ZERO PIVOT IS -1 // EXPECT_ROCSPARSE_STATUS(rocsparse_csritsv_zero_pivot(handle, descr, info, analysis_pivot), rocsparse_status_success); analysis_no_pivot.unit_check(analysis_pivot); if(verbose) { std::cout << " - compute device iterative, calling solve before analysis must fail" << std::endl; } // // Call before analysis // host_nmaxiter[0] = arg.nmaxiter; EXPECT_ROCSPARSE_STATUS(rocsparse_csritsv_solve(PARAMS_SOLVE(h_alpha, dA, dx, dy)), rocsparse_status_invalid_pointer); if(verbose) { std::cout << " - compute device iterative, call analysis" << std::endl; } // Call it twice, for code coverage. CHECK_ROCSPARSE_ERROR(rocsparse_csritsv_analysis(PARAMS_ANALYSIS(dA))); if(verbose) { std::cout << " - compute device iterative, call analysis twice (second call for coverage)" << std::endl; } CHECK_ROCSPARSE_ERROR(rocsparse_csritsv_analysis(PARAMS_ANALYSIS(dA))); { auto st = rocsparse_csritsv_zero_pivot(handle, descr, info, analysis_pivot); EXPECT_ROCSPARSE_STATUS(st, (*analysis_pivot != -1) ? rocsparse_status_zero_pivot : rocsparse_status_success); } if(verbose) { std::cout << " - compute device iterative, call solve" << std::endl; } CHECK_HIP_ERROR(hipMemset(dy, 0, sizeof(T) * M)); host_nmaxiter[0] = arg.nmaxiter; CHECK_ROCSPARSE_ERROR(rocsparse_csritsv_solve(PARAMS_SOLVE(h_alpha, dA, dx, dy))); { auto st = rocsparse_csritsv_zero_pivot(handle, descr, info, solve_pivot); EXPECT_ROCSPARSE_STATUS( st, (*solve_pivot != -1) ? rocsparse_status_zero_pivot : rocsparse_status_success); } const bool device_iterative_convergence = (host_nmaxiter[0] < arg.nmaxiter); h_analysis_pivot.unit_check(analysis_pivot); h_solve_pivot.unit_check(solve_pivot); if(*h_analysis_pivot == -1 && *h_solve_pivot == -1) { if(verbose) { std::cout << " - compute device iterative, compare solutions between host and device." << std::endl; } if(host_iterative_convergence) { hy_iterative.near_check(dy, tol_compare); } } host_nmaxiter[0] = arg.nmaxiter; CHECK_HIP_ERROR(hipMemset(dy, 0, sizeof(T) * M)); CHECK_ROCSPARSE_ERROR(rocsparse_csritsv_solve_ex(PARAMSX_SOLVE(h_alpha, dA, dx, dy))); { auto st = rocsparse_csritsv_zero_pivot(handle, descr, info, solve_pivot); CHECK_HIP_ERROR(hipDeviceSynchronize()); EXPECT_ROCSPARSE_STATUS( st, (*solve_pivot != -1) ? rocsparse_status_zero_pivot : rocsparse_status_success); } const bool device_iterative_convergence_ex = (host_nmaxiter[0] < arg.nmaxiter); if(*h_analysis_pivot == -1 && *h_solve_pivot == -1) { if(verbose) { std::cout << " - compute device iterative_ex, compare solutions between host and device." << std::endl; } if(host_iterative_convergence) { hy_iterative.near_check(dy, tol_compare); } } // // ============================================================================ // Regression test for rocsparse_copy_mat_info on a csritsv info. // // For a general matrix, csritsv analysis allocates device buffers (the ptr_end // submatrix offsets and the zero-pivot position) that are owned by the // rocsparse_csritsv_info and freed in its destructor. Copying a mat_info used to // be broken in two ways that this test guards against: // // 1. ptr_end was copied shallowly, so a copied rocsparse_mat_info shared the // SAME device allocation as its source. Destroying the source freed that // buffer out from under the copy (use-after-free), and destroying both // objects double-freed it. The fix deep-copies the owned buffer so each // object owns its allocation. This test reproduces the exact scenario: // copy, destroy the source, keep using the copy, then destroy the copy. // // 2. The pivot/position metadata was copied using the destination's (invalid, // default) index type instead of the source's, corrupting the copied info // so that solving with it produced a wrong result. The solve-and-compare // below deterministically fails unless that metadata is copied correctly. // ============================================================================ // if(arg.matrix_type == rocsparse_matrix_type_general && *h_analysis_pivot == -1 && *h_solve_pivot == -1 && host_iterative_convergence) { rocsparse_pointer_mode pointer_mode_save; CHECK_ROCSPARSE_ERROR(rocsparse_get_pointer_mode(handle, &pointer_mode_save)); CHECK_ROCSPARSE_ERROR(rocsparse_set_pointer_mode(handle, rocsparse_pointer_mode_host)); rocsparse_mat_info src_info = nullptr; rocsparse_mat_info copy_info = nullptr; CHECK_ROCSPARSE_ERROR(rocsparse_create_mat_info(&src_info)); CHECK_ROCSPARSE_ERROR(rocsparse_create_mat_info(©_info)); // // Analyze the source info; this allocates the owned ptr_end device buffer. // CHECK_ROCSPARSE_ERROR(rocsparse_csritsv_analysis(handle, trans, dA.m, dA.nnz, descr, dA.val, dA.ptr, dA.ind, src_info, apol, spol, dbuffer)); // // Copy the analyzed source info into a freshly created destination info. // CHECK_ROCSPARSE_ERROR(rocsparse_copy_mat_info(copy_info, src_info)); // // Destroy the source. With the shallow-copy bug this frees the device buffer // still referenced by copy_info (use-after-free); with the deep copy copy_info // keeps its own buffer. The copy must remain fully usable afterwards. // CHECK_ROCSPARSE_ERROR(rocsparse_destroy_mat_info(src_info)); // // Solve with the (now sole) copy and compare against the reference iterative // solution. This fails if the pivot/position metadata was not copied with the // correct index type, and exercises the copy's ptr_end buffer after the source // has been destroyed. // CHECK_HIP_ERROR(hipMemset(dy, 0, sizeof(T) * M)); host_nmaxiter[0] = arg.nmaxiter; CHECK_ROCSPARSE_ERROR(rocsparse_csritsv_solve(handle, host_nmaxiter, host_tol, host_history, trans, dA.m, dA.nnz, h_alpha, descr, dA.val, dA.ptr, dA.ind, copy_info, dx, dy, spol, dbuffer)); hy_iterative.near_check(dy, tol_compare); // // Destroying the copy frees its ptr_end buffer. With the deep copy this is the // buffer's only owner and the free is clean; with the shallow-copy bug this is // a second free of the already-freed source buffer (double-free). // CHECK_ROCSPARSE_ERROR(rocsparse_destroy_mat_info(copy_info)); CHECK_HIP_ERROR(hipDeviceSynchronize()); CHECK_ROCSPARSE_ERROR(rocsparse_set_pointer_mode(handle, pointer_mode_save)); } if(device_iterative_convergence != host_iterative_convergence) { if(verbose) { if(trans == rocsparse_operation_none) { std::cout << " - WARNING host and device iterative convergence differs on " "NonTranspose case ... it happens rarely but it does." << std::endl; } else { std::cout << " - WARNING host and device iterative convergence differs, it happens " "with transpose cases since csrmv uses atomics." << std::endl; } } } if(device_iterative_convergence_ex != host_iterative_convergence) { if(verbose) { if(trans == rocsparse_operation_none) { std::cout << " - WARNING host and device iterative_ex convergence differs on " "NonTranspose case ... it happens rarely but it does." << std::endl; } else { std::cout << " - WARNING host and device iterative_ex convergence differs, it " "happens " "with transpose cases since csrmv uses atomics." << std::endl; } } } if(verbose) { std::cout << " - compute device iterative, check pivot." << std::endl; std::cout << "h_analysis_pivot " << *h_analysis_pivot << std::endl; std::cout << "h_solve_pivot " << *h_solve_pivot << std::endl; } // // RESET MAT INFO. // info.reset(); if(verbose) { std::cout << "testing_csritsv unit_check pointer mode device ... " << std::endl; } if(verbose) { std::cout << " - compute device iterative" << std::endl; } device_scalar d_analysis_pivot; device_scalar d_solve_pivot; device_scalar d_alpha(h_alpha); // // POINTER MODE DEVICE // CHECK_ROCSPARSE_ERROR(rocsparse_set_pointer_mode(handle, rocsparse_pointer_mode_device)); // // CHECK IF DEFAULT ZERO PIVOT IS -1 // EXPECT_ROCSPARSE_STATUS(rocsparse_csritsv_zero_pivot(handle, descr, info, d_analysis_pivot), rocsparse_status_success); analysis_no_pivot.unit_check(d_analysis_pivot); // // Call before analysis // host_nmaxiter[0] = arg.nmaxiter; EXPECT_ROCSPARSE_STATUS(rocsparse_csritsv_solve(PARAMS_SOLVE(h_alpha, dA, dx, dy)), rocsparse_status_invalid_pointer); // // Call it twice. // CHECK_ROCSPARSE_ERROR(rocsparse_csritsv_analysis(PARAMS_ANALYSIS(dA))); CHECK_ROCSPARSE_ERROR(rocsparse_csritsv_analysis(PARAMS_ANALYSIS(dA))); EXPECT_ROCSPARSE_STATUS(rocsparse_csritsv_zero_pivot(handle, descr, info, d_analysis_pivot), (*h_analysis_pivot != -1) ? rocsparse_status_zero_pivot : rocsparse_status_success); CHECK_HIP_ERROR(hipMemset(dy, 0, sizeof(T) * M)); host_nmaxiter[0] = arg.nmaxiter; CHECK_ROCSPARSE_ERROR(rocsparse_csritsv_solve(PARAMS_SOLVE(d_alpha, dA, dx, dy))); EXPECT_ROCSPARSE_STATUS(rocsparse_csritsv_zero_pivot(handle, descr, info, d_solve_pivot), (*h_solve_pivot != -1) ? rocsparse_status_zero_pivot : rocsparse_status_success); h_analysis_pivot.unit_check(d_analysis_pivot); h_solve_pivot.unit_check(d_solve_pivot); if(*h_analysis_pivot == -1 && *h_solve_pivot == -1) { if(ROCSPARSE_REPRODUCIBILITY) { rocsparse_reproducibility::save("Y", dy); } if(host_iterative_convergence) { hy_iterative.near_check(dy, tol_compare); } } else { if(ROCSPARSE_REPRODUCIBILITY) { rocsparse_reproducibility::save( "analysis_pivot", d_analysis_pivot, "solve_pivot", d_solve_pivot); } } CHECK_HIP_ERROR(hipMemset(dy, 0, sizeof(T) * M)); host_nmaxiter[0] = arg.nmaxiter; CHECK_ROCSPARSE_ERROR(rocsparse_csritsv_solve_ex(PARAMSX_SOLVE(d_alpha, dA, dx, dy))); EXPECT_ROCSPARSE_STATUS(rocsparse_csritsv_zero_pivot(handle, descr, info, d_solve_pivot), (*h_solve_pivot != -1) ? rocsparse_status_zero_pivot : rocsparse_status_success); h_analysis_pivot.unit_check(d_analysis_pivot); h_solve_pivot.unit_check(d_solve_pivot); if(*h_analysis_pivot == -1 && *h_solve_pivot == -1) { if(host_iterative_convergence) { hy_iterative.near_check(dy, tol_compare); } } free(hbuffer); if(verbose) { std::cout << "testing_csritsv unit_check done. " << std::endl; } } if(arg.timing) { host_nfreeiter = arg.nfreeiter; int number_cold_calls = 2; int number_hot_calls = arg.iters; CHECK_ROCSPARSE_ERROR(rocsparse_set_pointer_mode(handle, rocsparse_pointer_mode_host)); // Warm up for(int iter = 0; iter < number_cold_calls; ++iter) { CHECK_ROCSPARSE_ERROR(rocsparse_csritsv_analysis(PARAMS_ANALYSIS(dA))); CHECK_ROCSPARSE_ERROR( rocsparse_csritsv_zero_pivot(handle, descr, info, h_analysis_pivot)); CHECK_ROCSPARSE_ERROR( rocsparse_csritsv_solve_ex(PARAMSX_SOLVE(h_alpha, dA, dx, dy))); CHECK_ROCSPARSE_ERROR(rocsparse_csritsv_zero_pivot(handle, descr, info, h_solve_pivot)); CHECK_ROCSPARSE_ERROR(rocsparse_csritsv_clear(handle, descr, info)); } double gpu_analysis_time_used = get_time_us(); CHECK_ROCSPARSE_ERROR(rocsparse_csritsv_analysis(PARAMS_ANALYSIS(dA))); gpu_analysis_time_used = get_time_us() - gpu_analysis_time_used; double gpu_solve_time_used = 0; // Performance run for(int iter = 0; iter < number_hot_calls; ++iter) { CHECK_HIP_ERROR(hipMemset(dy, 0, sizeof(T) * M)); host_nmaxiter[0] = arg.nmaxiter; double gpu_solve_time_used_1 = get_time_us(); CHECK_ROCSPARSE_ERROR( rocsparse_csritsv_solve_ex(PARAMSX_SOLVE(h_alpha, dA, dx, dy))); gpu_solve_time_used_1 = (get_time_us() - gpu_solve_time_used_1); gpu_solve_time_used += gpu_solve_time_used_1; } gpu_solve_time_used /= number_hot_calls; double gflop_count = csrsv_gflop_count(M, dA.nnz, diag); double gbyte_count = csrsv_gbyte_count(M, dA.nnz); double gpu_gflops = get_gpu_gflops(gpu_solve_time_used, gflop_count); double gpu_gbyte = get_gpu_gbyte(gpu_solve_time_used, gbyte_count); display_timing_info(display_key_t::M, M, display_key_t::nnz, dA.nnz, display_key_t::alpha, *h_alpha, display_key_t::nfreeiter, host_nfreeiter, display_key_t::nmaxiter, host_nmaxiter[0], display_key_t::pivot, std::min(*h_analysis_pivot, *h_solve_pivot), display_key_t::trans, rocsparse_operation2string(trans), display_key_t::diag_type, rocsparse_diagtype2string(diag), display_key_t::fill_mode, rocsparse_fillmode2string(uplo), display_key_t::analysis_policy, rocsparse_analysis2string(apol), display_key_t::solve_policy, rocsparse_solve2string(spol), display_key_t::gflops, gpu_gflops, display_key_t::bandwidth, gpu_gbyte, display_key_t::analysis_time_ms, get_gpu_time_msec(gpu_analysis_time_used), display_key_t::time_ms, get_gpu_time_msec(gpu_solve_time_used)); } // Clear csritsv meta data CHECK_ROCSPARSE_ERROR(rocsparse_csritsv_clear(handle, descr, info)); // Free buffer CHECK_HIP_ERROR(rocsparse_hipFree(dbuffer)); } #define INSTANTIATE(TYPE) \ template void testing_csritsv_bad_arg(const Arguments& arg); \ template void testing_csritsv(const Arguments& arg) INSTANTIATE(float); INSTANTIATE(double); INSTANTIATE(rocsparse_float_complex); INSTANTIATE(rocsparse_double_complex); void testing_csritsv_extra(const Arguments& arg) {}