From 08797b4e99abb5b9599c8f5396a78437fc2e99fc Mon Sep 17 00:00:00 2001 From: Stack-1 Date: Fri, 1 Aug 2025 15:03:54 +0200 Subject: [PATCH] [ADD] Added psb_gedots tests, even if tests using matrix do not work yet --- test/computational_routines/README.md | 2 +- .../gedot/psb_gedot_test.f90 | 7 +- test/computational_routines/gedots/Makefile | 37 + test/computational_routines/gedots/README.md | 49 ++ .../computational_routines/gedots/autotest.sh | 40 ++ .../gedots/psb_gedots_test.f90 | 677 ++++++++++++++++++ .../utils/psb_test_utils.f90 | 13 +- util/psb_c_mat_dist_impl.f90 | 20 +- util/psb_d_mat_dist_impl.f90 | 20 +- util/psb_s_mat_dist_impl.f90 | 20 +- util/psb_z_mat_dist_impl.f90 | 20 +- 11 files changed, 852 insertions(+), 53 deletions(-) create mode 100644 test/computational_routines/gedots/Makefile create mode 100644 test/computational_routines/gedots/README.md create mode 100755 test/computational_routines/gedots/autotest.sh create mode 100644 test/computational_routines/gedots/psb_gedots_test.f90 diff --git a/test/computational_routines/README.md b/test/computational_routines/README.md index e6aec7a0..1947d895 100644 --- a/test/computational_routines/README.md +++ b/test/computational_routines/README.md @@ -32,7 +32,7 @@ In this test suite were considered only computational routines implemented by PS | ------------------------------- | :--------------------------: | ---------------------------------------------------------------------- | :---------------: |:---------------: |:---------------: |:---------------: | |**General Dense Matrix Sum**| `psb_geaxpby`| This subroutine is an interface to the computational kernel for dense matrix sum: $Y \leftarrow \alpha X + \beta Y$ |Yes ✅|Yes ✅|No ❌|No ❌| | **Dot product**|`psb_gedot`|This function computes dot product between two vectors x and y. $dot \leftarrow x^T y$ If x and y are real vectors it computes dot-product as: $dot \leftarrow x^H y$ |Yes ✅|Yes ✅|No ❌|No ❌| -| **Generalized Dot Product** |`psb_gedots`|This subroutine computes a series of dot products among the columns of two dense matrices x and y:$res(i) \leftarrow x(:,i)^T y(:,i)$ If the matrices are complex, then the usual convention applies, i.e. the conjugate transpose of x is used. If x and y are of rank one, then res is a scalar, else it is a rank one array.|No ❌|No ❌|No ❌|No ❌| +| **Generalized Dot Product** |`psb_gedots`|This subroutine computes a series of dot products among the columns of two dense matrices x and y:$res(i) \leftarrow x(:,i)^T y(:,i)$ If the matrices are complex, then the usual convention applies, i.e. the conjugate transpose of x is used. If x and y are of rank one, then res is a scalar, else it is a rank one array.|Yes ✅|Yes ✅|No ❌|No ❌| |**Infinity-Norm of Vector**|`psb_normi`/`psb_geamax`|This function computes the infinity-norm of a vector x. If x is a real vector it computes infinity norm as:$amax \leftarrow max \mid x_i \mid$else if x is a complex vector then it computes the infinity-norm as: $amax \leftarrow max(\mid re(x_i) \mid + \mid im(x_i) \mid)$ |No ❌|No ❌|No ❌|No ❌| |**Generalized Infinity Norm**|`psb_geamaxs`|This subroutine computes a series of infinity norms on the columns of a dense matrix x: $res(i) \leftarrow max_k \mid x(k,i) \mid$ |No ❌|No ❌|No ❌|No ❌| | **1-Norm of Vector**| `psb_norm1` / `psb_geasums`|This function computes the 1-norm of a vector x. If x is a real vector it computes 1-norm as: $asum \leftarrow \mid \mid x_i \mid \mid$ else if x is a complex vector then it computes 1-norm as: $asum \leftarrow \mid \mid re(x) \mid \mid_1 + \mid \mid im(x) \mid \mid_1$ |No ❌|No ❌|No ❌|No ❌| diff --git a/test/computational_routines/gedot/psb_gedot_test.f90 b/test/computational_routines/gedot/psb_gedot_test.f90 index d13f25fe..652a5032 100644 --- a/test/computational_routines/gedot/psb_gedot_test.f90 +++ b/test/computational_routines/gedot/psb_gedot_test.f90 @@ -41,7 +41,7 @@ !! Scope: local !! Type: required !! Intent: in -!! Specified as: an object of type psb desc type. +!! Specified as: an object of type psb_desc_type. !! !! global Descritption: Specifies whether the computation should include the global !! reduction across all processes. @@ -208,16 +208,13 @@ contains type(psb_desc_type) :: desc_a ! communication context - integer(psb_ipk_) :: my_rank, np, info, err_act + integer(psb_ipk_) :: my_rank, np, info ! variables outside PSLBALS data structures real(psb_spk_), allocatable :: x_single_global(:), y_single_global(:) real(psb_dpk_), allocatable :: x_double_global(:), y_double_global(:) integer(psb_ipk_) :: i, nl - ! others - logical :: exists - info = psb_success_ call psb_info(ctxt,my_rank,np) diff --git a/test/computational_routines/gedots/Makefile b/test/computational_routines/gedots/Makefile new file mode 100644 index 00000000..400cb415 --- /dev/null +++ b/test/computational_routines/gedots/Makefile @@ -0,0 +1,37 @@ +INSTALLDIR=../../.. +INCDIR=$(INSTALLDIR)/include/ +MODDIR=$(INSTALLDIR)/modules/ +include $(INCDIR)/Make.inc.psblas +# +# Libraries used +# +LIBDIR = $(INSTALLDIR)/lib/ +PSBLAS_LIB = -L$(LIBDIR) -lpsb_util -lpsb_base +LDLIBS = $(PSBLDLIBS) +EXEDIR = ./runs + +FINCLUDES=$(FMFLAG)$(MODDIR) $(FMFLAG). + +GREEN=\033[0;32m +RED=\033[0;31m +BLUE=\033[0;34m +YELLOW=\033[33m +END_COLOUR=\033[0m + + +all: runsd psb_gedots_test + @printf "$(GREEN)[INFO]\t Compilation success!$(END_COLOUR)\n" + +runsd: + @(if test ! -d $(EXEDIR) ; then mkdir $(EXEDIR); fi) + @printf "$(BLUE)[INFO]\t Build directory $(EXEDIR) correctly initialized$(END_COLOUR)\n" + +psb_gedots_test: + @$(FLINK) $(LOPT) psb_gedots_test.f90 ../utils/psb_test_utils.o -o $(EXEDIR)/psb_gedots_test -I../utils/ -I$(MODDIR) -I. $(PSBLAS_LIB) $(LDLIBS) + @printf "$(BLUE)[INFO]\t Testing files for psb_gedot linked and compiled correctly$(END_COLOUR)\n" + +clean: + @rm -f $(OBJS)\ + *$(.mod) $(EXEDIR)/psb_gedots_test + +.PHONY: all psb_gedots_test clean diff --git a/test/computational_routines/gedots/README.md b/test/computational_routines/gedots/README.md new file mode 100644 index 00000000..26445719 --- /dev/null +++ b/test/computational_routines/gedots/README.md @@ -0,0 +1,49 @@ +# Introduction +This is a directory developed by Luca Pepè Sciarria and Simone Staccone from Tor Vergata University to start to create some unit tests for PSBLAS 3.9, in particular for ```psb_gedot``` routine. + + +## Getting started +Steps to reproduce the tests: + +- Compile the code using ``` make ``` (Optional) +- Launch the script ```./autotest.sh``` or with ```source ./autotest.sh``` if you want to add modules to the .bashrc file permanently. +- Check the output log file ```psb_gedot_test.log``` to collect results and check for errors. In any case a summarization of tests passed should be shown in stdout. + +NOTE: If the code is changed and a new compilation is needed to show the changes, the ```autotest.sh``` script isn't aware of this scenario, therefore it is necessary to manually recompile the code. Moreover, if you are running manually the script and not launching the main ```test.sh``` script, be careful to use the last compiled version of the utils. + +## Test Suite +The ```psb_gedot``` is a function that performs the scalar product between two vectors giving a value as result. The signature of the function is: + +```fortran +psb_gedot(x, y, desc_a, info [,global]) +``` + +The strategy to validate the correctness of the computation is to compare single precision result and double precision result in the test cases in which the test should not give an error. In this way it is possible to have a correctness check of the computation comparing the two results considering a number of significant digits which is tuned on the single precision computation. + +### Parameters Values +We considered high condition number as 10^6, so we will consider only matrices having a lower estimated condition number, in particular the technique used to estimate it is the SVD (Singular Value Decomposition) considering the highest and lowest value on the diagonal (?). +**x** vectors are located in the vectors/ directory. They are generated randomly using the same seed and then saved on different files based on their characteristics. The size of the vector is choosen accordingly to the size of the matrix column space considered for the single test instance. +|Variable|File Name|Coefficients|Coefficients Description| +|:-:|:-:|:-:|:-:| +|$x_1$|x1.txt|$x_i> 0, \forall i$|Positive coefficients| +|$x_2$|x2.txt|$x_i < 0, \forall i$|Negative coefficients +|$x_3$|x3.txt|$x_i \ne 0, \forall i$|Random coefficients +|$x_4$|x4.txt|$x_i = 0, \forall i$|Null coefficients + +**y** vectors are located in the vectors/ directory. They are generated randomly using the same seed and then saved on different files based on their characteristics. The size of the vector is choosen accordingly to the size of the matrix rows space considered for the single test instance. +|Variable|File Name|Coefficients|Coefficients Description| +|:-:|:-:|:-:|:-:| +|$y_1$|y1.txt|$y_i> 0, \forall i$|Positive coefficients| +|$y_2$|y2.txt|$y_i < 0, \forall i$|Negative coefficients +|$y_3$|y3.txt|$y_i \ne 0, \forall i$|Random coefficients +|$y_4$|y4.txt|$y_i = 0, \forall i$|Null coefficients + + + +## TODO +List of things still to add in the test: +- Add computation with broken descriptor and catch the errore result (Use EXCECUTE_COMMAND_LINE from a fortran program and check the exit codes) +- Test using complex data ($dot \leftarrow x^H \cdot y$) +- Test also GPU excecution +- Try multiple distributions +- Add test using matrix (the first attempt was unsuccessful) \ No newline at end of file diff --git a/test/computational_routines/gedots/autotest.sh b/test/computational_routines/gedots/autotest.sh new file mode 100755 index 00000000..6e9dd9e7 --- /dev/null +++ b/test/computational_routines/gedots/autotest.sh @@ -0,0 +1,40 @@ +#!/bin/bash + +# Variables definition +num_procs=$(nproc) + +# Define color codes +GREEN="\033[0;32m" +RED="\033[0;31m" +BLUE="\033[0;34m" +YELLOW="\033[33m" +RESET="\033[0m" + + +# Check if the executable ELF file exists +if [ ! -f "./runs/psb_gedots_test" ]; then + echo -e "${YELLOW}[WARNING] Executable not found. Running make...${RESET}" + make + if [ ! -f "./runs/psb_geaxpbys_test" ]; then + echo -e "${RED}[ERROR] Failed to create executable. Check make command.${RESET}" + fi +else + echo -e "${BLUE}[INFO]\t The executable already exists. Skipping the make process.${RESET}" +fi + + +# Excecute tests and save results +echo -e "${BLUE}[INFO]\t Running the PSBLAS psb_gedots test...${RESET}" +echo "" +echo -e "${BLUE}[INFO]\t Starting single process computation${RESET}" +mpirun -np 1 ./runs/psb_gedots_test +echo -e "${BLUE}[INFO]\t Single process computation terminated correctly${RESET}" +echo "" +echo -e "${BLUE}[INFO]\t Starting $num_procs processes computation${RESET}" +mpirun -np $num_procs ./runs/psb_gedots_test +echo -e "${BLUE}[INFO]\t Multiple processes computation terminated correctly${RESET}" + +rm -f results/* +rm -f -r results/ + +echo -e "${GREEN}[INFO]\t PSBLAS psb_gedots test succesfully completed.${RESET}" \ No newline at end of file diff --git a/test/computational_routines/gedots/psb_gedots_test.f90 b/test/computational_routines/gedots/psb_gedots_test.f90 new file mode 100644 index 00000000..3d67d6a2 --- /dev/null +++ b/test/computational_routines/gedots/psb_gedots_test.f90 @@ -0,0 +1,677 @@ +!> Test program for y = x^T * y or y = x^H * y psb_gedot routine +!! Check the README.md to see all details about the tests. +!! +!! Authors: Luca Pepé Sciarria, Staccone Simone (Tor Vergata University) +!! +!! psb_gedot(x, y, desc_a, info [,global]) +!! +!! Type: Synchronous. +!! +!! ====================================== +!! | Data type | Precision | +!! ====================================== +!! | psb_spk_ | Short Precision Real | +!! | psb_dpk_ | Long Precision Real | +!! | psb_cpk_ | Short Precision Complex| +!! | psb_zpk_ | Long Precision Complex | +!! ====================================== +!! Table 1: Data types +!! +!! ROUTINE PARAMETERS +!! +!! Input: +!! +!! x Description: the local portion of global dense matrix x. +!! Scope: local +!! Type: required +!! Intent: in +!! Specified as: a rank one or two array or an object of type psb_T_vect_type +!! containing numbers of type specified in Table 1. The rank of x must be +!! the same of y. +!! +!! y Description: the local portion of the global dense matrix y. +!! Scope: local +!! Type: required +!! Intent: inout +!! Specified as: a rank one or two array or an object of type psb_T_vect_type +!! containing numbers of the type indicated in Table 1. The rank of y must +!! be the same of x. +!! +!! desc_a Description: contains data structures for communications. +!! Scope: local +!! Type: required +!! Intent: in +!! Specified as: an object of type psb_desc_type. +!! +!! Output: +!! +!! res Description: is the dot product of vectors x and y +!! Scope: global +!! Intent: out +!! Specified as: a number or a rank-one array of the data type indicated in Table 1 +!! +!! info Description: Error code. +!! Scope: local +!! Type: required +!! Intent: out +!! Specified as: An integer value; 0 means no error has been detected. +!! +program main + use psb_base_mod + use psb_util_mod + use psb_test_utils + + implicit none + + ! Communicator variable + type(psb_ctxt_type) :: ctxt + + ! parameters array + character(len=64) :: x(4),y(4) + integer(psb_ipk_) :: arr_size + integer(psb_ipk_) :: tests_number, count + + ! cycle indexes variables + integer(psb_ipk_) :: i,j,k,h,l + integer(psb_ipk_) :: info, unit + + ! results + real(psb_spk_) :: result_single + real(psb_dpk_) :: result_double + real(psb_spk_), allocatable :: result_single_vector(:) + real(psb_dpk_), allocatable :: result_double_vector(:) + type(psb_test_info) :: test_info + + + + ! Initialize parameters + x(1) = "vectors/x1.mtx" + x(2) = "vectors/x2.mtx" + x(3) = "vectors/x3.mtx" + x(4) = "vectors/x4.mtx" + + y(1) = "vectors/y1.mtx" + y(2) = "vectors/y2.mtx" + y(3) = "vectors/y3.mtx" + y(4) = "vectors/y4.mtx" + + arr_size = 100000 + + !! Initialize test metadata + test_info%total_tests = size(x) * size(y) + test_info%threshold_type = GAMMA + test_info%threshold = 1.0D-06 + test_info%kernel_name = "psb_gedots" + + call psb_test_init(test_info) + + if(test_info%my_rank == psb_root_) then + psb_out_unit = test_info%output_unit + call psb_test_generate_input_vectors(arr_size) + end if + + call psb_bcast(test_info%ctxt,test_info%output_unit) + call psb_barrier(test_info%ctxt) + + + if(test_info%my_rank == psb_root_) write(*,'(A)') "[INFO] Starting test excecution ..." + + ! Iterate over test parameters + do i=1,size(x) + do j=1,size(y) + call psb_gedots_real_kernel(x(i), y(j), arr_size, test_info%ctxt,result_single, result_double) + + if(test_info%my_rank == psb_root_) then + if(test_info%np > 1) then + ! If the program is being run on multiple processes, we need to + ! check the result on the root process with the one computed only using + ! a single process + call psb_test_process_check(result_single, test_info) + else + call psb_test_single_double_scalar_check(result_single,result_double,test_info, arr_size) + + ! If the program is being run on a single process, we can save the result directly + call psb_test_save_result(result_single, test_info) + end if + test_info%current_test = test_info%current_test + 1 + end if + + call psb_barrier(test_info%ctxt) + end do + end do + + + ! Test using matrices + x(1) = "../matrix/1138_bus.mtx" + x(2) = "../matrix/crystm03.mtx" + x(3) = "../matrix/qc2534.mtx" + x(4) = "../matrix/rdb5000.mtx" + + y(1) = "../matrix/1138_bus.mtx" + y(2) = "../matrix/crystm03.mtx" + y(3) = "../matrix/qc2534.mtx" + y(4) = "../matrix/rdb5000.mtx" + + + + if(test_info%my_rank == psb_root_) then + allocate(result_single_vector(arr_size)) + allocate(result_double_vector(arr_size)) + end if + + ! Iterate over test parameters + !! do i=1,size(x) + !! do j=1,size(y) + !! call psb_gedots_real_matrix_kernel(x(i), y(j), arr_size, test_info%ctxt,result_single_vector, result_double_vector) +!! + !! if(test_info%my_rank == psb_root_) then + !! if(test_info%np > 1) then + !! ! If the program is being run on multiple processes, we need to + !! ! check the result on the root process with the one computed only using + !! ! a single process + !! call psb_test_process_vector_check(result_single_vector, test_info) + !! else + !! call psb_test_single_double_vector_check(result_single_vector,result_double_vector,test_info, arr_size) + !! + !! ! If the program is being run on a single process, we can save the result directly + !! call psb_test_save_vector_result(result_single_vector, test_info) + !! end if + !! test_info%current_test = test_info%current_test + 1 + !! end if +!! + !! call psb_barrier(test_info%ctxt) + !! end do + !! end do + + if(test_info%my_rank == psb_root_) then + deallocate(result_single_vector) + deallocate(result_double_vector) + end if + + call psb_test_exit(test_info) + +contains + + + subroutine psb_gedots_real_matrix_kernel(x_file, y_file, arr_size, ctxt, result_single, result_double) + ! input parameters + character(len = *), intent(in) :: x_file, y_file + integer(psb_ipk_), intent(in) :: arr_size + type(psb_ctxt_type), intent(in) :: ctxt + + ! output parameters + real(psb_spk_), allocatable, intent(inout) :: result_single(:) + real(psb_dpk_), allocatable, intent(inout) :: result_double(:) + + ! matrices + type(psb_lsspmat_type) :: x_single_aux, y_single_aux + type(psb_ldspmat_type) :: x_double_aux, y_double_aux + + ! sparse matrices + type(psb_sspmat_type) :: x_single, y_single + type(psb_dspmat_type) :: x_double, y_double + + ! matrix parameters + integer(psb_ipk_) :: x_rows, x_cols, x_nnz + integer(psb_ipk_) :: y_rows, y_cols, y_nnz + integer(psb_ipk_) :: nl + + ! matrix descriptor data structure + type(psb_desc_type) :: desc_a + + ! communication context + integer(psb_ipk_) :: my_rank, np, info + + + + + ! Allocate dense arrays for the local part only + integer(psb_ipk_) :: num_local_rows_x, num_local_cols_x, num_local_rows_y, num_local_cols_y + real(psb_spk_), allocatable :: x_dense_single(:,:), y_dense_single(:,:) + real(psb_dpk_), allocatable :: x_dense_double(:,:), y_dense_double(:,:) + ! For single precision x + integer(psb_ipk_), allocatable :: row_x(:), col_x(:) + real(psb_spk_), allocatable :: val_x(:) + integer(psb_ipk_) :: nnz_x, idx + type(psb_s_coo_sparse_mat) :: x_coo_single + ! For single precision y + integer(psb_ipk_), allocatable :: row_y(:), col_y(:) + real(psb_spk_), allocatable :: val_y(:) + integer(psb_ipk_) :: nnz_y + type(psb_s_coo_sparse_mat) :: y_coo_single + + ! For double precision x + integer(psb_ipk_), allocatable :: row_xd(:), col_xd(:) + real(psb_dpk_), allocatable :: val_xd(:) + integer(psb_ipk_) :: nnz_xd + type(psb_d_coo_sparse_mat) :: x_coo_d + ! For double precision y + integer(psb_ipk_), allocatable :: row_yd(:), col_yd(:) + real(psb_dpk_), allocatable :: val_yd(:) + integer(psb_ipk_) :: nnz_yd + type(psb_d_coo_sparse_mat) :: y_coo_d + + + + class(psb_s_coo_sparse_mat), pointer :: a_ptr_s + class(psb_d_coo_sparse_mat), pointer :: a_ptr_d + + + + + + + call psb_info(ctxt,my_rank,np) + + if (my_rank < 0) then + ! This should not happen, but just in case + call psb_error(ctxt) + endif + + call mm_mat_read(a=x_single_aux,info=info, iunit=17, filename=x_file) + if(info /= psb_success_) then + write(psb_out_unit,'(A)') "Error reading single precision matrix x" + return + end if + + call mm_mat_read(a=x_double_aux,info=info, iunit=17, filename=x_file) + if(info /= psb_success_) then + write(psb_out_unit,'(A)') "Error reading double precision matrix x" + return + end if + + call mm_mat_read(a=y_single_aux,info=info, iunit=17, filename=y_file) + if(info /= psb_success_) then + write(psb_out_unit,'(A)') "Error reading single precision matrix y" + return + end if + + call mm_mat_read(a=y_double_aux,info=info, iunit=17, filename=y_file) + if(info /= psb_success_) then + write(psb_out_unit,'(A)') "Error reading double precision matrix y" + return + end if + + + x_rows = x_single_aux%get_nrows() + x_cols = x_single_aux%get_ncols() + x_nnz = x_single_aux%get_nzeros() + + y_rows = y_single_aux%get_nrows() + y_cols = y_single_aux%get_ncols() + y_nnz = y_single_aux%get_nzeros() + + + ! Allocate descriptor as if it was a block rows distribution + nl = (arr_size)/np + mod(arr_size,np) + + + ! part_block it's a macro defined in psb_blockpart_mod to identify BLOCK ROWS distribution + call psb_matdist(x_single_aux, x_single, ctxt,desc_a,info,fmt="COO",parts=part_block) + if(info /= psb_success_) then + write(psb_out_unit,'(A)') "Error in psb_matdist for single precision matrix x" + goto 9999 + end if + + call psb_matdist(x_double_aux, x_double, ctxt,desc_a,info,fmt="COO",parts=part_block) + if(info /= psb_success_) then + write(psb_out_unit,'(A)') "Error in psb_matdist for double precision matrix x" + goto 9999 + end if + + call psb_matdist(y_single_aux, y_single, ctxt,desc_a,info,fmt="COO",parts=part_block) + if(info /= psb_success_) then + write(psb_out_unit,'(A)') "Error in psb_matdist for single precision matrix y" + goto 9999 + end if + + call psb_matdist(y_double_aux, y_double, ctxt,desc_a,info,fmt="COO",parts=part_block) + if(info /= psb_success_) then + write(psb_out_unit,'(A)') "Error in psb_matdist for double precision matrix y" + goto 9999 + end if + + + + + + ! --------------------------------------------------------------------- ! + + + + ! Get local dimensions from the descriptor + num_local_rows_x = x_single%get_nrows() + num_local_cols_x = x_single%get_ncols() + num_local_rows_y = y_single%get_nrows() + num_local_cols_y = y_single%get_ncols() + + allocate(x_dense_single(num_local_rows_x, num_local_cols_x)) + allocate(y_dense_single(num_local_rows_y, num_local_cols_y)) + allocate(x_dense_double(num_local_rows_x, num_local_cols_x)) + allocate(y_dense_double(num_local_rows_y, num_local_cols_y)) + + ! Initialize to zero + x_dense_single = 0.0_psb_spk_ + y_dense_single = 0.0_psb_spk_ + x_dense_double = 0.0_psb_dpk_ + y_dense_double = 0.0_psb_dpk_ + + ! Explicitly extract COO values and insert them in the dense arrays + + + + + select type(a_ptr_s => x_single%a) + type is (psb_s_coo_sparse_mat) + x_coo_single = a_ptr_s + class default + write(psb_out_unit,'(A)') "Error: x_single%a is not of type psb_s_coo_sparse_mat" + return + end select + + nnz_x = x_coo_single%get_nzeros() + allocate(row_x(nnz_x), col_x(nnz_x), val_x(nnz_x)) + row_x = x_coo_single%ia + col_x = x_coo_single%ja + val_x = x_coo_single%val + + do idx = 1, nnz_x + x_dense_single(row_x(idx), col_x(idx)) = val_x(idx) + end do + + + + select type(a_ptr_s => y_single%a) + type is (psb_s_coo_sparse_mat) + y_coo_single = a_ptr_s + class default + write(psb_out_unit,'(A)') "Error: x_single%a is not of type psb_s_coo_sparse_mat" + return + end select + + nnz_y = y_coo_single%get_nzeros() + allocate(row_y(nnz_y), col_y(nnz_y), val_y(nnz_y)) + row_y = y_coo_single%ia + col_y = y_coo_single%ja + val_y = y_coo_single%val + do idx = 1, nnz_y + y_dense_single(row_y(idx), col_y(idx)) = val_y(idx) + end do + + + + select type(a_ptr_d => x_double%a) + type is (psb_d_coo_sparse_mat) + x_coo_d = a_ptr_d + class default + write(psb_out_unit,'(A)') "Error: x_single%a is not of type psb_s_coo_sparse_mat" + return + end select + + + nnz_xd = x_coo_d%get_nzeros() + allocate(row_xd(nnz_xd), col_xd(nnz_xd), val_xd(nnz_xd)) + row_xd = x_coo_d%ia + col_xd = x_coo_d%ja + val_xd = x_coo_d%val + + do idx = 1, nnz_xd + x_dense_double(row_xd(idx), col_xd(idx)) = val_xd(idx) + end do + + + + + + select type(a_ptr_d => y_double%a) + type is (psb_d_coo_sparse_mat) + y_coo_d = a_ptr_d + class default + write(psb_out_unit,'(A)') "Error: x_single%a is not of type psb_s_coo_sparse_mat" + return + end select + + + nnz_yd = y_coo_d%get_nzeros() + allocate(row_yd(nnz_yd), col_yd(nnz_yd), val_yd(nnz_yd)) + row_yd = y_coo_d%ia + col_yd = y_coo_d%ja + val_yd = y_coo_d%val + do idx = 1, nnz_yd + y_dense_double(row_yd(idx), col_yd(idx)) = val_yd(idx) + end do + + + + + ! --------------------------------------------------------------------- ! + + ! y = x^T * y + call psb_gedots(result_single,x_dense_single,y_dense_single,desc_a,info) + if(info /= psb_success_) then + write(psb_out_unit,'(A)') "Error in psb_gedots routine in single precision" + goto 9999 + end if +!! + !! ! y = x^T * y + call psb_gedots(result_double,x_dense_double,y_dense_double,desc_a,info) + if(info /= psb_success_) then + write(psb_out_unit,'(A)') "Error in psb_gedots routine in double precision" + goto 9999 + end if + + + + + + + + + 9999 call psb_spfree(x_single, desc_a,info) + if(info /= psb_success_) then + write(psb_out_unit,'(A)') "Error in single precision vector x free routine" + end if + + call psb_spfree(y_single, desc_a,info) + if(info /= psb_success_) then + write(psb_out_unit,'(A)') "Error in single precision vector y free routine" + end if + + call psb_spfree(x_double, desc_a,info) + if(info /= psb_success_) then + write(psb_out_unit,'(A)') "Error in double precision vector x free routine" + end if + + call psb_spfree(y_double, desc_a,info) + if(info /= psb_success_) then + write(psb_out_unit,'(A)') "Error in double precision vector y free routine" + end if + + call psb_cdfree(desc_a,info) + if(info /= psb_success_) then + write(psb_out_unit,'(A)') "Error in matrix descriptor free routine" + end if + + + end subroutine + + + !> @brief Function to excecute psb_gedots in single precision real + !! vector and compare with the same computation in double + !! precision + !! + !! @param x_file file name of the first vector + !! @param y_file file name of the second vector + !! @param arr_size size of the vectors + !! @param ctxt communication context + !! @param result_single result of the single precision computation + !! @param result_double result of the double precision computation + !! + subroutine psb_gedots_real_kernel(x_file, y_file, arr_size, ctxt, result_single, result_double) + ! input parameters + character(len = *), intent(in) :: x_file, y_file + integer(psb_ipk_), intent(in) :: arr_size + type(psb_ctxt_type), intent(in) :: ctxt + + ! output parameters + real(psb_spk_), intent(out) :: result_single + real(psb_dpk_), intent(out) :: result_double + + ! vectors + type(psb_s_vect_type) :: x_single, y_single + type(psb_d_vect_type) :: x_double, y_double + + ! matrix descriptor data structure + type(psb_desc_type) :: desc_a + + ! communication context + integer(psb_ipk_) :: my_rank, np, info, err_act + + ! variables outside PSLBALS data structures + real(psb_spk_), allocatable :: x_single_global(:), y_single_global(:) + real(psb_dpk_), allocatable :: x_double_global(:), y_double_global(:) + integer(psb_ipk_) :: i, nl + + ! others + logical :: exists + + info = psb_success_ + + call psb_info(ctxt,my_rank,np) + + if (my_rank < 0) then + ! This should not happen, but just in case + call psb_error(ctxt) + endif + + ! Generate random array for b using always the same seed + if(my_rank == psb_root_) then + allocate(x_single_global(arr_size)) + allocate(y_single_global(arr_size)) + allocate(x_double_global(arr_size)) + allocate(y_double_global(arr_size)) + + call mm_array_read(x_single_global,info,filename=x_file) + call mm_array_read(y_single_global,info,filename=y_file) + call mm_array_read(x_double_global,info,filename=x_file) + call mm_array_read(y_double_global,info,filename=y_file) + end if + + ! Allocate descriptor as if it was a block rows distribution + nl = (arr_size)/np + mod(arr_size,np) + + call psb_cdall(ctxt, desc_a, info,nl=nl) + if(info /= psb_success_) then + write(psb_out_unit,'(A)') "Error allocating desc_a data structure" + goto 9999 + end if + + call psb_cdasb(desc_a, info) + if(info /= psb_success_) then + write(psb_out_unit,'(A)') "Error assembling desc_a data structure" + goto 9999 + end if + + + call psb_geall(x_single,desc_a,info) + if(info /= psb_success_) then + write(psb_out_unit,'(A)') "Error allocating single precision x data structure" + goto 9999 + end if + + call psb_geall(x_double,desc_a,info) + if(info /= psb_success_) then + write(psb_out_unit,'(A)') "Error allocating double precision x data structure" + goto 9999 + end if + + ! Populate x class using data from x_global vector + call psb_scatter(x_single_global,x_single,desc_a,info,root=psb_root_) + if(info /= psb_success_) then + write(psb_out_unit,'(A)') "Error in psb_scatter to populate single precision x data structure" + goto 9999 + end if + + call psb_scatter(x_double_global,x_double,desc_a,info,root=psb_root_) + if(info /= psb_success_) then + write(psb_out_unit,'(A)') "Error in psb_scatter to populate double precision x data structure" + goto 9999 + end if + + call psb_geall(y_single,desc_a,info) + if(info /= psb_success_) then + write(psb_out_unit,'(A)') "Error allocating single precision y data structure" + goto 9999 + end if + + call psb_geall(y_double,desc_a,info) + if(info /= psb_success_) then + write(psb_out_unit,'(A)') "Error allocating double precision y data structure" + goto 9999 + end if + + ! Populate y class using data from y_global vector + call psb_scatter(y_single_global,y_single,desc_a,info,root=psb_root_) + if(info /= psb_success_) then + write(psb_out_unit,'(A)') "Error in psb_scatter to populate single precision y data structure" + goto 9999 + end if + + call psb_scatter(y_double_global,y_double,desc_a,info,root=psb_root_) + if(info /= psb_success_) then + write(psb_out_unit,'(A)') "Error in psb_scatter to populate double precision y data structure" + goto 9999 + end if + + ! y = x^T * y + call psb_gedots(result_single,x_single%v%v,y_single%v%v,desc_a,info) + if(info /= psb_success_) then + write(psb_out_unit,'(A)') "Error in psb_gedots routine in single precision" + goto 9999 + end if + + + call psb_gedots(result_double,x_double%v%v,y_double%v%v,desc_a,info) + if(info /= psb_success_) then + write(psb_out_unit,'(A)') "Error in psb_gedots routine in double precision" + goto 9999 + end if + + + ! Deallocate + 9999 call psb_gefree(x_single, desc_a,info) + if(info /= psb_success_) then + write(psb_out_unit,'(A)') "Error in single precision vector x free routine" + end if + + call psb_gefree(y_single, desc_a,info) + if(info /= psb_success_) then + write(psb_out_unit,'(A)') "Error in single precision vector y free routine" + end if + + call psb_gefree(x_double, desc_a,info) + if(info /= psb_success_) then + write(psb_out_unit,'(A)') "Error in double precision vector x free routine" + end if + + call psb_gefree(y_double, desc_a,info) + if(info /= psb_success_) then + write(psb_out_unit,'(A)') "Error in double precision vector y free routine" + end if + + + call psb_cdfree(desc_a,info) + if(info /= psb_success_) then + write(psb_out_unit,'(A)') "Error in matrix descriptor free routine" + end if + + if(my_rank == 0) then + deallocate(x_single_global) + deallocate(y_single_global) + deallocate(x_double_global) + deallocate(y_double_global) + end if + + end subroutine + +end program main \ No newline at end of file diff --git a/test/computational_routines/utils/psb_test_utils.f90 b/test/computational_routines/utils/psb_test_utils.f90 index 809d6f64..6b4a5bec 100644 --- a/test/computational_routines/utils/psb_test_utils.f90 +++ b/test/computational_routines/utils/psb_test_utils.f90 @@ -354,10 +354,10 @@ contains else call psb_test_log_failed(test_info, out_string) test_info%failure = test_info%failure + 1 + write(psb_out_unit,'(A,F20.10)') "Threshold used: ", test_info%threshold end if write(psb_out_unit,'(A,F20.10)') "Single precision result: ", result_single write(psb_out_unit,'(A,F20.10)') "Double precision result: ", result_double - write(psb_out_unit,'(A,F20.10)') "Threshold used: ", test_info%threshold end subroutine !> @brief Subroutine to check the results of a single and double precision vector computation. @@ -391,7 +391,7 @@ contains else call psb_test_log_failed(test_info, out_string) test_info%failure = test_info%failure + 1 - write(psb_out_unit,'(A,F20.10)') "Comparison error occurred at index: ", i + write(psb_out_unit,'(A,I0)') "Comparison error occurred at index: ", i write(psb_out_unit,'(A,F20.10)') "Single precision result: ", result_single(i) write(psb_out_unit,'(A,F20.10)') "Double precision result: ", result_double(i) end if @@ -483,12 +483,11 @@ contains else call psb_test_log_failed(test_info, out_string) test_info%failure = test_info%failure + 1 + write(test_info%output_unit, '(A,F20.10)') "Delta: ", abs(saved_result - result_single) end if - write(test_info%output_unit, '(F20.10,F20.10,A,L,A,L)') & - & local_saved - local_single, local_single - local_saved, " ", local_saved - local_single == 0, & - & " ", local_single - local_saved == 0 - write(test_info%output_unit, '(A,F20.10)') "Multi-process result: ", local_single - write(test_info%output_unit, '(A,F20.10)') "Single process result: ", local_saved + + write(test_info%output_unit, '(A,F20.10)') "Multi-process result: ", result_single + write(test_info%output_unit, '(A,F20.10)') "Single process result: ", saved_result end subroutine diff --git a/util/psb_c_mat_dist_impl.f90 b/util/psb_c_mat_dist_impl.f90 index 26d04041..d28ce0ff 100644 --- a/util/psb_c_mat_dist_impl.f90 +++ b/util/psb_c_mat_dist_impl.f90 @@ -152,10 +152,10 @@ subroutine psb_cmatdist(a_glob, a, ctxt, desc_a,& call psb_errpush(info,name,i_err=(/liwork/),a_err='integer') goto 9999 endif - if (iam == root) then - write (*, fmt = *) 'start matdist',root, size(iwork),& - &nrow, ncol, nnzero,nrhs - endif + !! if (iam == root) then + !! write (*, fmt = *) 'start matdist',root, size(iwork),& + !! &nrow, ncol, nnzero,nrhs + !! endif if (use_parts) then call psb_cdall(ctxt,desc_a,info,mg=nrow,parts=parts) else if (use_vg) then @@ -359,7 +359,7 @@ subroutine psb_cmatdist(a_glob, a, ctxt, desc_a,& goto 9999 end if - if (iam == root) write (*, fmt = *) 'end matdist' + !! if (iam == root) write (*, fmt = *) 'end matdist' call psb_erractionrestore(err_act) return @@ -490,10 +490,10 @@ subroutine psb_lcmatdist(a_glob, a, ctxt, desc_a,& call psb_errpush(info,name,l_err=(/liwork/),a_err='integer') goto 9999 endif - if (iam == root) then - write (*, fmt = *) 'start matdist',root, size(iwork),& - &nrow, ncol, nnzero,nrhs - endif + !! if (iam == root) then + !! write (*, fmt = *) 'start matdist',root, size(iwork),& + !! &nrow, ncol, nnzero,nrhs + !! endif if (use_parts) then call psb_cdall(ctxt,desc_a,info,mg=nrow,parts=parts) else if (use_vg) then @@ -700,7 +700,7 @@ subroutine psb_lcmatdist(a_glob, a, ctxt, desc_a,& goto 9999 end if - if (iam == root) write (*, fmt = *) 'end matdist' + !! if (iam == root) write (*, fmt = *) 'end matdist' call psb_erractionrestore(err_act) return diff --git a/util/psb_d_mat_dist_impl.f90 b/util/psb_d_mat_dist_impl.f90 index 44d853e4..286dcaf6 100644 --- a/util/psb_d_mat_dist_impl.f90 +++ b/util/psb_d_mat_dist_impl.f90 @@ -152,10 +152,10 @@ subroutine psb_dmatdist(a_glob, a, ctxt, desc_a,& call psb_errpush(info,name,i_err=(/liwork/),a_err='integer') goto 9999 endif - if (iam == root) then - write (*, fmt = *) 'start matdist',root, size(iwork),& - &nrow, ncol, nnzero,nrhs - endif + !! if (iam == root) then + !! write (*, fmt = *) 'start matdist',root, size(iwork),& + !! &nrow, ncol, nnzero,nrhs + !! endif if (use_parts) then call psb_cdall(ctxt,desc_a,info,mg=nrow,parts=parts) else if (use_vg) then @@ -359,7 +359,7 @@ subroutine psb_dmatdist(a_glob, a, ctxt, desc_a,& goto 9999 end if - if (iam == root) write (*, fmt = *) 'end matdist' + !! if (iam == root) write (*, fmt = *) 'end matdist' call psb_erractionrestore(err_act) return @@ -490,10 +490,10 @@ subroutine psb_ldmatdist(a_glob, a, ctxt, desc_a,& call psb_errpush(info,name,l_err=(/liwork/),a_err='integer') goto 9999 endif - if (iam == root) then - write (*, fmt = *) 'start matdist',root, size(iwork),& - &nrow, ncol, nnzero,nrhs - endif + !! if (iam == root) then + !! write (*, fmt = *) 'start matdist',root, size(iwork),& + !! &nrow, ncol, nnzero,nrhs + !! endif if (use_parts) then call psb_cdall(ctxt,desc_a,info,mg=nrow,parts=parts) else if (use_vg) then @@ -700,7 +700,7 @@ subroutine psb_ldmatdist(a_glob, a, ctxt, desc_a,& goto 9999 end if - if (iam == root) write (*, fmt = *) 'end matdist' + !! if (iam == root) write (*, fmt = *) 'end matdist' call psb_erractionrestore(err_act) return diff --git a/util/psb_s_mat_dist_impl.f90 b/util/psb_s_mat_dist_impl.f90 index 126938bf..39575144 100644 --- a/util/psb_s_mat_dist_impl.f90 +++ b/util/psb_s_mat_dist_impl.f90 @@ -152,10 +152,10 @@ subroutine psb_smatdist(a_glob, a, ctxt, desc_a,& call psb_errpush(info,name,i_err=(/liwork/),a_err='integer') goto 9999 endif - if (iam == root) then - write (*, fmt = *) 'start matdist',root, size(iwork),& - &nrow, ncol, nnzero,nrhs - endif + !! if (iam == root) then + !! write (*, fmt = *) 'start matdist',root, size(iwork),& + !! &nrow, ncol, nnzero,nrhs + !! endif if (use_parts) then call psb_cdall(ctxt,desc_a,info,mg=nrow,parts=parts) else if (use_vg) then @@ -359,7 +359,7 @@ subroutine psb_smatdist(a_glob, a, ctxt, desc_a,& goto 9999 end if - if (iam == root) write (*, fmt = *) 'end matdist' + !! if (iam == root) write (*, fmt = *) 'end matdist' call psb_erractionrestore(err_act) return @@ -490,10 +490,10 @@ subroutine psb_lsmatdist(a_glob, a, ctxt, desc_a,& call psb_errpush(info,name,l_err=(/liwork/),a_err='integer') goto 9999 endif - if (iam == root) then - write (*, fmt = *) 'start matdist',root, size(iwork),& - &nrow, ncol, nnzero,nrhs - endif + !! if (iam == root) then + !! write (*, fmt = *) 'start matdist',root, size(iwork),& + !! &nrow, ncol, nnzero,nrhs + !! endif if (use_parts) then call psb_cdall(ctxt,desc_a,info,mg=nrow,parts=parts) else if (use_vg) then @@ -700,7 +700,7 @@ subroutine psb_lsmatdist(a_glob, a, ctxt, desc_a,& goto 9999 end if - if (iam == root) write (*, fmt = *) 'end matdist' + !! if (iam == root) write (*, fmt = *) 'end matdist' call psb_erractionrestore(err_act) return diff --git a/util/psb_z_mat_dist_impl.f90 b/util/psb_z_mat_dist_impl.f90 index 48abe46d..58e1b6e2 100644 --- a/util/psb_z_mat_dist_impl.f90 +++ b/util/psb_z_mat_dist_impl.f90 @@ -152,10 +152,10 @@ subroutine psb_zmatdist(a_glob, a, ctxt, desc_a,& call psb_errpush(info,name,i_err=(/liwork/),a_err='integer') goto 9999 endif - if (iam == root) then - write (*, fmt = *) 'start matdist',root, size(iwork),& - &nrow, ncol, nnzero,nrhs - endif + !! if (iam == root) then + !! write (*, fmt = *) 'start matdist',root, size(iwork),& + !! &nrow, ncol, nnzero,nrhs + !! endif if (use_parts) then call psb_cdall(ctxt,desc_a,info,mg=nrow,parts=parts) else if (use_vg) then @@ -359,7 +359,7 @@ subroutine psb_zmatdist(a_glob, a, ctxt, desc_a,& goto 9999 end if - if (iam == root) write (*, fmt = *) 'end matdist' + !! if (iam == root) write (*, fmt = *) 'end matdist' call psb_erractionrestore(err_act) return @@ -490,10 +490,10 @@ subroutine psb_lzmatdist(a_glob, a, ctxt, desc_a,& call psb_errpush(info,name,l_err=(/liwork/),a_err='integer') goto 9999 endif - if (iam == root) then - write (*, fmt = *) 'start matdist',root, size(iwork),& - &nrow, ncol, nnzero,nrhs - endif + !! if (iam == root) then + !! write (*, fmt = *) 'start matdist',root, size(iwork),& + !! &nrow, ncol, nnzero,nrhs + !! endif if (use_parts) then call psb_cdall(ctxt,desc_a,info,mg=nrow,parts=parts) else if (use_vg) then @@ -700,7 +700,7 @@ subroutine psb_lzmatdist(a_glob, a, ctxt, desc_a,& goto 9999 end if - if (iam == root) write (*, fmt = *) 'end matdist' + !! if (iam == root) write (*, fmt = *) 'end matdist' call psb_erractionrestore(err_act) return