You cannot select more than 25 topics
Topics must start with a letter or number, can include dashes ('-') and can be up to 35 characters long.
677 lines
23 KiB
Fortran
677 lines
23 KiB
Fortran
!> 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 |