diff --git a/meson.build b/meson.build index f9411c9..541c06b 100644 --- a/meson.build +++ b/meson.build @@ -251,3 +251,42 @@ if build_tests is_parallel : false) endforeach endif + +# C tests +if build_tests + + fortran_tests_folder = 'tests/C' + + foreach test: cutest_c_tests + package = test[0] + precision = test[1] + name = test[2] + file = test[3] + + if precision == 'single' + libgalahad_precision = libcutest_single + args_precision = pp_flag + '-DREAL_32' + endif + if precision == 'double' + libgalahad_precision = libcutest_double + args_precision = pp_flag + endif + if precision == 'quadruple' + libgalahad_precision = libcutest_quadruple + args_precision = pp_flag + '-DREAL_128' + '-DCUTEST_16btye_reals_exist' + endif + + test(name, + executable(name, file, + fortran_args : args_precision, + c_args : args_precision, + link_with : libgalahad_precision, + link_language : 'fortran', + include_directories: libcutest_include, + install : true, + install_dir : fortran_tests_folder), + suite : [package, precision, 'C'], + workdir : join_paths(meson.project_source_root(), 'src', 'test'), + is_parallel : false) + endforeach +endif diff --git a/src/test/meson.build b/src/test/meson.build index d3f28ff..df2ce45 100644 --- a/src/test/meson.build +++ b/src/test/meson.build @@ -9,10 +9,15 @@ cutest_tests += [['cutest', 'single', 'ctest_single' , files('ctest.F90' ['cutest', 'double', 'utest_threaded_double', files('utest_threaded.F90', 'u_elfun_double.f', 'u_group_double.f', 'u_range_double.f')], ['cutest', 'double', 'lqp_test_double' , files('lqptest.F90' , 'q_elfun_double.f', 'q_group_double.f', 'q_range_double.f')]] +cutest_c_tests += [['cutest', 'single', 'utest_c_single', files('utest.c', 'u_elfun_single.f', 'u_group_single.f', 'u_range_single.f')], + ['cutest', 'double', 'utest_c_double', files('utest.c', 'u_elfun_double.f', 'u_group_double.f', 'u_range_double.f')]] + if build_quadruple cutest_tests += [['cutest', 'quadruple', 'ctest_quadruple' , files('ctest.F90' , 'c_elfun_quadruple.f', 'c_group_quadruple.f', 'c_range_quadruple.f')], ['cutest', 'quadruple', 'ctest_threaded_quadruple', files('ctest_threaded.F90', 'c_elfun_quadruple.f', 'c_group_quadruple.f', 'c_range_quadruple.f')], ['cutest', 'quadruple', 'utest_quadruple' , files('utest.F90' , 'u_elfun_quadruple.f', 'u_group_quadruple.f', 'u_range_quadruple.f')], ['cutest', 'quadruple', 'utest_threaded_quadruple', files('utest_threaded.F90', 'u_elfun_quadruple.f', 'u_group_quadruple.f', 'u_range_quadruple.f')], ['cutest', 'quadruple', 'lqp_test_quadruple' , files('lqptest.F90' , 'q_elfun_quadruple.f', 'q_group_quadruple.f', 'q_range_quadruple.f')]] + + cutest_c_tests += [['cutest', 'quadruple', 'utest_c_quadruple', files('utest.c', 'u_elfun_quadruple.f', 'u_group_quadruple.f', 'u_range_quadruple.f')]] endif diff --git a/src/test/utest.F90 b/src/test/utest.F90 index ca93a19..035c193 100644 --- a/src/test/utest.F90 +++ b/src/test/utest.F90 @@ -275,18 +275,10 @@ PROGRAM CUTEST_test_unconstrained_tools ALLOCATE( H_band( 0 : lbandh, n ), stat = alloc_stat ) IF ( alloc_stat /= 0 ) GO TO 990 - goth = .FALSE. - WRITE( out, "( ' Call CUTEST_ubandh with goth = .FALSE.' )" ) - CALL CUTEST_ubandh_r( status, n, X, nsemib, H_band, lbandh, maxsbw ) - IF ( status /= 0 ) GO to 900 - CALL WRITE_H_BAND( out, n, lbandh, H_band, nsemib ) -! CALL WRITE_H_BAND( out, n, lbandh, H_band, nsemib, maxsbw ) - goth = .TRUE. - WRITE( out, "( ' Call CUTEST_ubandh with goth = .TRUE.' )" ) + WRITE( out, "( ' Call CUTEST_ubandh' )" ) CALL CUTEST_ubandh_r( status, n, X, nsemib, H_band, lbandh, maxsbw ) IF ( status /= 0 ) GO to 900 CALL WRITE_H_BAND( out, n, lbandh, H_band, nsemib ) -! CALL WRITE_H_BAND( out, n, lbandh, H_band, nsemib, maxsbw ) ! calls and time report diff --git a/src/test/utest.c b/src/test/utest.c new file mode 100644 index 0000000..940614e --- /dev/null +++ b/src/test/utest.c @@ -0,0 +1,533 @@ +#include +#include +#include + +#include "cutest_routines.h" +#include "cutest.h" + +#define INPUT_FILE "u_OUTSDIF.d" + +// Function prototypes +void write_x(ipc_ n, rpc_ *X, rpc_ *X_l, rpc_ *X_u); +void write_x_type(ipc_ n, ipc_ *X_type); +void write_p_name(char *p_name); +void write_x_names(ipc_ n, char **X_names); +void write_f(rpc_ f); +void write_g(ipc_ n, rpc_ *G); +void write_h_dense(ipc_ n, ipc_ l_h2_1, rpc_ **H2_val); +void write_h_sparsity_pattern(ipc_ H_ne, ipc_ l_h, ipc_ *H_row, ipc_ *H_col); +void write_h_sparse(ipc_ H_ne, ipc_ l_h, rpc_ *H_val, ipc_ *H_row, ipc_ *H_col); +void write_h_element(ipc_ ne, ipc_ lhe_ptr, ipc_ *HE_row_ptr, ipc_ *HE_val_ptr, ipc_ *lhe_row, ipc_ *HE_row, ipc_ *lhe_val, rpc_ *HE_val); +void write_result(ipc_ n, rpc_ *vector, rpc_ *result); +void write_sresult(ipc_ n, ipc_ nnz_vector, ipc_ *INDEX_nz_vector, rpc_ *vector, ipc_ nnz_result, ipc_ *INDEX_nz_result, rpc_ *result); +void write_h_band(ipc_ n, ipc_ lbandh, rpc_ **H_band, ipc_ nsemib); + +int main() { + // Parameters + ipc_ input = 55; + ipc_ out = 6; + ipc_ buffer = 77; + rpc_ zero = 0.0; + rpc_ one = 1.0; + + // Local Variables + ipc_ i, n, HE_nel, HE_val_ne, HE_row_ne, status; + ipc_ l_h2_1, l_h, lhe_ptr, H_ne, lhe_val, lhe_row; + ipc_ nnz_vector, nnz_result, maxsbw, alloc_stat; + ipc_ nsemib, lbandh; + rpc_ f; + ipc_ grad, byrows; + bool goth; + char p_name[10]; + ipc_ *X_type; + ipc_ *H_row, *H_col; + ipc_ *HE_row, *HE_row_ptr; + ipc_ *HE_val_ptr; + ipc_ *INDEX_nz_vector; + rpc_ *INDEX_nz_result; + rpc_ *X, *X_l, *X_u, *G; + rpc_ *H_val, *HE_val; + rpc_ *vector, *result; + rpc_ **H2_val, **H_band; + char **X_names; + rpc_ CPU[4], CALLS[4]; + + // Open the problem data file + FILE *input_file = fopen(INPUT_FILE, "r"); + if (input_file == NULL) { + perror("Error opening input file"); + return 1; + } + + // Allocate basic arrays + printf("Call CUTEST_udimen\n"); + CUTEST_udimen_r(&status, &input, &n); + printf("* n = %d\n", n); + l_h2_1 = n; + + X = malloc(n * sizeof(rpc_)); + X_l = malloc(n * sizeof(rpc_)); + X_u = malloc(n * sizeof(rpc_)); + G = malloc(n * sizeof(rpc_)); + vector = malloc(n * sizeof(rpc_)); + result = malloc(n * sizeof(rpc_)); + X_names = malloc(n * sizeof(char *)); + X_type = malloc(n * sizeof(ipc_)); + INDEX_nz_vector = malloc(n * sizeof(ipc_)); + INDEX_nz_result = malloc(n * sizeof(ipc_)); + if (X == NULL || X_l == NULL || X_u == NULL || G == NULL || + vector == NULL || result == NULL || X_names == NULL || + X_type == NULL || INDEX_nz_vector == NULL || INDEX_nz_result == NULL) { + perror("Error allocating memory"); + return 1; + } + + H2_val = malloc(l_h2_1 * sizeof(ipc_ *)); + for (i = 0; i < l_h2_1; ++i) { + H2_val[i] = malloc(n * sizeof(rpc_)); + } + + if (H2_val == NULL) { + perror("Error allocating memory"); + return 1; + } + + // Set up SIF data + printf("Call CUTEST_usetup\n"); + CUTEST_usetup_r(&status, &input, &out, &buffer, &n, X, X_l, X_u); + if (status != 0) { + printf("error status = %d\n", status); + } + write_x(n, X, X_l, X_u); + + // Obtain variable and problem names + printf("Call CUTEST_unames\n"); + CUTEST_unames_r(&status, &n, p_name, X_names); + if (status != 0) { + printf("error status = %d\n", status); + } + write_p_name(p_name); + write_x_names(n, X_names); + + // Obtain problem name + printf("Call CUTEST_probname\n"); + CUTEST_probname_r(&status, p_name); + if (status != 0) { + printf("error status = %d\n", status); + } + write_p_name(p_name); + + // Obtain variable names + printf("Call CUTEST_varnames\n"); + CUTEST_varnames_r(&status, &n, X_names); + if (status != 0) { + printf("error status = %d\n", status); + } + write_x_names(n, X_names); + + // Obtain variable types + printf("Call CUTEST_uvartype\n"); + CUTEST_uvartype_r(&status, &n, X_type); + if (status != 0) { + printf("error status = %d\n", status); + } + write_x_type(n, X_type); + + // Compute the objective function value + printf("Call CUTEST_ufn\n"); + CUTEST_ufn_r(&status, &n, X, &f); + if (status != 0) { + printf("error status = %d\n", status); + } + write_f(f); + + // Compute the gradient value + printf("Call CUTEST_ugr\n"); + CUTEST_ugr_r(&status, &n, X, G); + if (status != 0) { + printf("error status = %d\n", status); + } + write_g(n, G); + + // Compute the objective function and gradient values + grad = 1; + printf("Call CUTEST_uofg with grad = TRUE\n"); + CUTEST_uofg_r(&status, &n, X, &f, G, grad); + if (status != 0) { + printf("error status = %d\n", status); + } + write_f(f); + write_g(n, G); + + grad = 0; + printf("Call CUTEST_uofg with grad = FALSE\n"); + CUTEST_uofg_r(&status, &n, X, &f, G, grad); + if (status != 0) { + printf("error status = %d\n", status); + } + write_f(f); + + // Compute the dense Hessian value + printf("Call CUTEST_udh\n"); + CUTEST_udh_r(&status, &n, X, &l_h2_1, H2_val); + if (status != 0) { + printf("error status = %d\n", status); + } + write_h_dense(n, l_h2_1, H2_val); + + // Compute the gradient and dense Hessian values + printf("Call CUTEST_ugrdh\n"); + CUTEST_ugrdh_r(&status, &n, X, G, &l_h2_1, H2_val); + if (status != 0) { + printf("error status = %d\n", status); + } + write_g(n, G); + write_h_dense(n, l_h2_1, H2_val); + + // Compute the number of nonzeros in the sparse Hessian + printf("Call CUTEST_udimsh\n"); + CUTEST_udimsh_r(&status, &H_ne); + if (status != 0) { + printf("error status = %d\n", status); + } + printf("* H_ne = %d\n", H_ne); + + l_h = H_ne; + H_val = malloc(l_h * sizeof(rpc_)); + H_row = malloc(l_h * sizeof(int)); + H_col = malloc(l_h * sizeof(int)); + if (H_val == NULL || H_row == NULL || H_col == NULL) { + perror("Error allocating memory"); + return 1; + } + + // Compute the sparsity pattern of the Hessian + printf("Call CUTEST_ushp\n"); + CUTEST_ushp_r(&status, &n, &H_ne, &l_h, H_row, H_col); + if (status != 0) { + printf("error status = %d\n", status); + } + write_h_sparsity_pattern(H_ne, l_h, H_row, H_col); + + // Compute the sparse Hessian value + printf("Call CUTEST_ush\n"); + CUTEST_ush_r(&status, &n, X, &H_ne, &l_h, H_val, H_row, H_col); + if (status != 0) { + printf("error status = %d\n", status); + } + write_h_sparse(H_ne, l_h, H_val, H_row, H_col); + + // Compute the gradient and sparse Hessian values + printf("Call CUTEST_ugrsh\n"); + CUTEST_ugrsh_r(&status, &n, X, G, &H_ne, &l_h, H_val, H_row, H_col); + if (status != 0) { + printf("error status = %d\n", status); + } + write_g(n, G); + write_h_sparse(H_ne, l_h, H_val, H_row, H_col); + + // Compute the number of nonzeros in the element Hessian + printf("Call CUTEST_udimse\n"); + CUTEST_udimse_r(&status, &HE_nel, &HE_val_ne, &HE_row_ne); + if (status != 0) { + printf("error status = %d\n", status); + } + printf("* H_nel = %d, HE_val_ne = %d, HE_row_ne = %d\n", HE_nel, HE_val_ne, HE_row_ne); + + lhe_ptr = HE_nel + 1; + lhe_val = HE_val_ne; + lhe_row = HE_row_ne; + HE_row_ptr = malloc(lhe_ptr * sizeof(int)); + HE_val_ptr = malloc(lhe_ptr * sizeof(int)); + HE_row = malloc(lhe_row * sizeof(int)); + HE_val = malloc(lhe_val * sizeof(rpc_)); + if (HE_row_ptr == NULL || HE_val_ptr == NULL || HE_row == NULL || HE_val == NULL) { + perror("Error allocating memory"); + return 1; + } + + byrows = false; + printf("Call CUTEST_ueh with byrows = .FALSE.\n"); + CUTEST_ueh_r(&status, &n, X, HE_nel, lhe_ptr, HE_row_ptr, HE_val_ptr, lhe_row, HE_row, lhe_val, HE_val, byrows); + if (status != 0) { + printf("error status = %d\n", status); + } + write_h_element(HE_nel, lhe_ptr, HE_row_ptr, HE_val_ptr, lhe_row, HE_row, lhe_val, HE_val); + + byrows = true; + printf("Call CUTEST_ueh with byrows = .TRUE.\n"); + CUTEST_ueh_r(&status, &n, X, &HE_nel, lhe_ptr, HE_row_ptr, HE_val_ptr, lhe_row, HE_row, lhe_val, HE_val, byrows); + if (status != 0) { + printf("error status = %d\n", status); + } + write_h_element(HE_nel, lhe_ptr, HE_row_ptr, HE_val_ptr, lhe_row, HE_row, lhe_val, HE_val); + + // Compute gradient and element Hessian values + byrows = false; + printf("Call CUTEST_ugreh with byrows = .FALSE.\n"); + CUTEST_ugreh_r(&status, &n, X, G, &HE_nel, lhe_ptr, HE_row_ptr, HE_val_ptr, lhe_row, HE_row, lhe_val, HE_val, byrows); + if (status != 0) { + printf("error status = %d\n", status); + } + write_g(n, G); + write_h_element(HE_nel, lhe_ptr, HE_row_ptr, HE_val_ptr, lhe_row, HE_row, lhe_val, HE_val); + + // byrows = true; + printf("Call CUTEST_ugreh with byrows = .TRUE.\n"); + CUTEST_ugreh_r(&status, n, X, G, HE_nel, lhe_ptr, HE_row_ptr, HE_val_ptr, lhe_row, HE_row, lhe_val, HE_val, byrows); + if (status != 0) { + printf("error status = %d\n", status); + } + write_g(n, G); + write_h_element(HE_nel, lhe_ptr, HE_row_ptr, HE_val_ptr, lhe_row, HE_row, lhe_val, HE_val); + + // Compute Hessian-vector product + vector[0] = 1.0; + for (int i = 1; i < n; i++) vector[i] = 0.0; + goth = false; + + printf("Call CUTEST_uhprod with goth = .FALSE.\n"); + CUTEST_uhprod_r(&status, &n, &goth, X, vector, result); + if (status != 0) { + printf("error status = %d\n", status); + } + write_result(n, vector, result); + + goth = true; + printf("Call CUTEST_uhprod with goth = .TRUE.\n"); + CUTEST_uhprod_r(&status, &n, &goth, X, vector, result); + if (status != 0) { + printf("error status = %d\n", status); + } + write_result(n, vector, result); + + // Compute sparse Hessian-vector product + nnz_vector = 1; + INDEX_nz_vector[0] = 1; + goth = false; + + printf("Call CUTEST_ushprod with goth = .FALSE.\n"); + CUTEST_ushprod_r(&status, &n, &goth, X, nnz_vector, INDEX_nz_vector, vector, &nnz_vector, INDEX_nz_vector, result); + if (status != 0) { + printf("error status = %d\n", status); + } + write_sresult(n, nnz_vector, INDEX_nz_vector, vector, nnz_vector, INDEX_nz_vector, result); + + goth = true; + printf("Call CUTEST_ushprod with goth = .TRUE.\n"); + CUTEST_ushprod_r(&status, &n, &goth, X, nnz_vector, INDEX_nz_vector, vector, &nnz_vector, INDEX_nz_vector, result); + if (status != 0) { + printf("error status = %d\n", status); + } + write_sresult(n, nnz_vector, INDEX_nz_vector, vector, nnz_vector, INDEX_nz_vector, result); + + printf("Call CUTEST_ubandh\n"); + CUTEST_ubandh_r(&status, &n, X, &nsemib, H_band, &lbandh, &maxsbw); + if (status != 0) { + printf("error status = %d\n", status); + } + write_h_band(n, lbandh, H_band, nsemib); + + // Calls and time report + printf("CALL CUTEST_ureport\n"); + CUTEST_ureport_r(&status, CALLS, CPU); + printf("CALLS(1-4) = %d %d %d %d\n", CALLS[0], CALLS[1], CALLS[2], CALLS[3]); + printf("CPU(1-4) = %.2f %.2f %.2f %.2f\n", CPU[0], CPU[1], CPU[2], CPU[3]); + + // Terminal exit + printf("Call CUTEST_uterminate\n"); + CUTEST_uterminate_r(&status); + if (status != 0) { + printf("error status = %d\n", status); + } + + // One more setup + printf("Call CUTEST_usetup\n"); + CUTEST_usetup_r(&status, &stdin, &out, &buffer, &n, X, X_l, X_u); + if (status != 0) { + printf("error status = %d\n", status); + } + + // Terminal exit again + printf("Call CUTEST_uterminate\n"); + CUTEST_uterminate_r(&status); + if (status != 0) { + printf("error status = %d\n", status); + } + + // Cleanup and close files + free(X); + free(X_l); + free(X_u); + free(G); + free(vector); + free(result); + free(X_names); + free(X_type); + free(INDEX_nz_vector); + free(INDEX_nz_result); + free(H2_val); + free(H_val); + free(H_row); + free(H_col); + free(HE_row_ptr); + free(HE_val_ptr); + free(HE_row); + free(HE_val); + + fclose(input_file); + + return 0; +} + +void write_x(ipc_ n, rpc_ *X, rpc_ *X_l, rpc_ *X_u) { + printf(" * i X_l X X_u\n"); + for (ipc_ i = 0; i < n; i++) { + printf(" * %7d %12.4f %12.4f %12.4f\n", i + 1, X_l[i], X[i], X_u[i]); + } +} + +void write_x_type(ipc_ n, ipc_ *X_type) { + printf(" * i X_type\n"); + for (ipc_ i = 0; i < n; i++) { + printf(" * %7d %10d\n", i + 1, X_type[i]); + } +} + +void write_p_name(char *p_name) { + printf(" * p_name = %s\n", p_name); +} + +void write_x_names(ipc_ n, char **X_names) { + printf(" * i X_name\n"); + for (ipc_ i = 0; i < n; i++) { + printf(" * %7d %10s\n", i + 1, X_names[i]); + } +} + +void write_f(rpc_ f) { + printf(" * f = %12.4f\n", f); +} + +void write_g(ipc_ n, rpc_ *G) { + printf(" * i G\n"); + for (ipc_ i = 0; i < n; i++) { + printf(" * %7d %12.4f\n", i + 1, G[i]); + } +} + +void write_h_dense(ipc_ n, ipc_ l_h2_1, rpc_ **H2_val) { + printf(" * H(dense)\n"); + for (ipc_ j = 0; j < n; j += 4) { + if (j + 3 < n) { + printf(" * i j%8d%8d%8d%8d\n", j, j + 1, j + 2, j + 3); + } else if (j + 2 < n) { + printf(" * i j%8d%8d%8d\n", j, j + 1, j + 2); + } else if (j + 1 < n) { + printf(" * i j%8d%8d\n", j, j + 1); + } else { + printf(" * i j%8d\n", j); + } + for (ipc_ i = 0; i < l_h2_1; i++) { + if (j + 3 < n) { + printf(" * %7d %8.4f %8.4f %8.4f %8.4f\n", i + 1, H2_val[i][j], H2_val[i][j + 1], H2_val[i][j + 2], H2_val[i][j + 3]); + } else if (j + 2 < n) { + printf(" * %7d %8.4f %8.4f %8.4f\n", i + 1, H2_val[i][j], H2_val[i][j + 1], H2_val[i][j + 2]); + } else if (j + 1 < n) { + printf(" * %7d %8.4f %8.4f\n", i + 1, H2_val[i][j], H2_val[i][j + 1]); + } else { + printf(" * %7d %8.4f\n", i + 1, H2_val[i][j]); + } + } + } +} + +void write_h_sparsity_pattern(ipc_ H_ne, ipc_ l_h, ipc_ *H_row, ipc_ *H_col) { + printf(" * H(sparse)\n"); + printf(" * row col\n"); + for (ipc_ i = 0; i < H_ne; i += 4) { + if (i + 3 < H_ne) { + printf(" * %7d %7d %7d %7d\n", H_row[i], H_col[i], H_row[i + 1], H_col[i + 1]); + printf(" * %7d %7d %7d %7d\n", H_row[i + 2], H_col[i + 2], H_row[i + 3], H_col[i + 3]); + } else if (i + 2 < H_ne) { + printf(" * %7d %7d %7d %7d\n", H_row[i], H_col[i], H_row[i + 1], H_col[i + 1]); + printf(" * %7d %7d\n", H_row[i + 2], H_col[i + 2]); + } else if (i + 1 < H_ne) { + printf(" * %7d %7d %7d %7d\n", H_row[i], H_col[i], H_row[i + 1], H_col[i + 1]); + } else { + printf(" * %7d %7d\n", H_row[i], H_col[i]); + } + } +} + +void write_h_sparse(ipc_ H_ne, ipc_ l_h, rpc_ *H_val, ipc_ *H_row, ipc_ *H_col) { + printf(" * H(sparse)\n"); + printf(" * row col val\n"); + for (ipc_ i = 0; i < H_ne; i += 2) { + if (i + 1 < H_ne) { + printf(" * %7d %7d %12.4f %7d %7d %12.4f\n", H_row[i], H_col[i], H_val[i], H_row[i + 1], H_col[i + 1], H_val[i + 1]); + } else { + printf(" * %7d %7d %12.4f\n", H_row[i], H_col[i], H_val[i]); + } + } +} + +void write_h_element(ipc_ ne, ipc_ lhe_ptr, ipc_ *HE_row_ptr, ipc_ *HE_val_ptr, ipc_ *lhe_row, ipc_ *HE_row, ipc_ *lhe_val, rpc_ *HE_val) { + printf(" * H(element)\n"); + for (ipc_ i = 0; i < ne; i++) { + if (HE_row_ptr[i + 1] > HE_row_ptr[i]) { + printf(" * element %d\n", i + 1); + printf(" * indices\n"); + for (ipc_ j = HE_row_ptr[i]; j < HE_row_ptr[i + 1]; j++) { + printf(" * %12d", HE_row[j]); + } + printf("\n * values\n"); + for (ipc_ j = HE_val_ptr[i]; j < HE_val_ptr[i + 1]; j++) { + printf(" * %12.4f", HE_val[j]); + } + printf("\n"); + } else { + printf(" * no indices\n"); + } + } +} + +void write_result(ipc_ n, rpc_ *vector, rpc_ *result) { + printf(" * i vector result\n"); + for (ipc_ i = 0; i < n; i++) { + printf(" * %7d %12.4f %12.4f\n", i + 1, vector[i], result[i]); + } +} + +void write_sresult(ipc_ n, ipc_ nnz_vector, ipc_ *INDEX_nz_vector, rpc_ *vector, ipc_ nnz_result, ipc_ *INDEX_nz_result, rpc_ *result) { + printf(" * i vector\n"); + for (ipc_ j = 0; j < nnz_vector; j++) { + ipc_ i = INDEX_nz_vector[j]; + printf(" * %7d %12.4f\n", i + 1, vector[i]); + } + printf(" * i result\n"); + for (ipc_ j = 0; j < nnz_result; j++) { + ipc_ i = INDEX_nz_result[j]; + printf(" * %7d %12.4f\n", i + 1, result[i]); + } +} + +void write_h_band(ipc_ n, ipc_ lbandh, rpc_ **H_band, ipc_ nsemib) { + printf(" * H(band)\n"); + printf(" * i band\n"); + for (ipc_ j = 0; j <= nsemib; j++) { + printf(" * %7d", j); + } + printf("\n"); + for (ipc_ i = 0; i < n; i++) { + printf(" * %7d", i); + for (ipc_ j = 0; j <= nsemib; j++) { + if (i - j >= 0 && i - j < n) { + printf(" %12.4f", H_band[i - j][j]); + } else { + printf(" %12s", "----"); + } + } + printf("\n"); + } +}