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.
250 lines
7.0 KiB
Fortran
250 lines
7.0 KiB
Fortran
!> Test program for y = AX spsb_pmm routine
|
|
!! Check the README.md to see all details about the tests.
|
|
!!
|
|
!! Author: Luca Pepé Sciarria, Staccone Simone (Tor Vergata University)
|
|
program psb_spmm_test
|
|
use psb_base_mod
|
|
use psb_util_mod
|
|
|
|
implicit none
|
|
|
|
! input parameters
|
|
character(len=256) :: mtx_file, file_name, name
|
|
|
|
! Testing parameters
|
|
character(len=256),dimension(2) :: mtx_files
|
|
|
|
! sparse matrices
|
|
type(psb_sspmat_type) :: a
|
|
type(psb_lsspmat_type) :: aux_a
|
|
|
|
! vectors
|
|
type(psb_s_vect_type) :: x, y
|
|
|
|
! matrix descriptor data structure
|
|
type(psb_desc_type) :: desc_a
|
|
|
|
! communication context
|
|
type(psb_ctxt_type) :: ctxt
|
|
integer(psb_ipk_) :: my_rank, np, info, err_act
|
|
|
|
! matrix parameters
|
|
integer(psb_ipk_) :: m, n, nnz
|
|
integer(psb_ipk_) :: nr, nt ! In BLOCK ROWS distributin, the number of rows
|
|
|
|
|
|
real(psb_spk_), allocatable :: x_global(:), y_global(:)
|
|
integer(psb_ipk_) :: i
|
|
|
|
name = "psb_spmm_test" ! Name of the program to output in case of error
|
|
info = psb_success_
|
|
mtx_file = "matrix/1138_bus.mtx"
|
|
mtx_files(1) = "matrix/1138_bus.mtx"
|
|
|
|
|
|
call psb_init(ctxt)
|
|
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
|
|
|
|
if (my_rank == psb_root_) then
|
|
write(psb_out_unit,*) 'Welcome to PSBLAS version: ',psb_version_string_
|
|
write(psb_out_unit,*) 'This is the ',trim(name),' sample program'
|
|
end if
|
|
|
|
|
|
call mm_mat_read(aux_a,info,filename=mtx_file)
|
|
if(info /= psb_success_) then
|
|
write(psb_out_unit,*) "Error while reading matric ", mtx_file
|
|
goto 9999
|
|
end if
|
|
|
|
|
|
! part_block it's a macro defined in psb_blockpart_mod to identify BLOCK ROWS distribution
|
|
call psb_matdist(aux_a, a, ctxt,desc_a,info,fmt="COO",parts=part_block)
|
|
|
|
|
|
m = aux_a%get_nrows()
|
|
n = aux_a%get_ncols()
|
|
nnz = aux_a%get_nzeros()
|
|
|
|
call psb_bcast(ctxt,m)
|
|
call psb_bcast(ctxt,n)
|
|
call psb_bcast(ctxt,nnz)
|
|
|
|
if(my_rank == psb_root_) then
|
|
write(psb_out_unit,*) "Matrix stats"
|
|
write(psb_out_unit,*) "ROWS:", m
|
|
write(psb_out_unit,*) "COLS:", n
|
|
write(psb_out_unit,*) "NNZ: ", nnz
|
|
end if
|
|
|
|
! Generate random array for b using always the same seed
|
|
if(my_rank == psb_root_) then
|
|
allocate(x_global(n))
|
|
allocate(y_global(n))
|
|
call generate_vectors(n) ! True for x
|
|
call mm_array_read(x_global,info,filename="vectors/x1.mtx")
|
|
call mm_array_read(y_global,info,filename="vectors/y1.mtx")
|
|
end if
|
|
|
|
call psb_geall(x,desc_a,info)
|
|
if(info /= psb_success_) then
|
|
write(psb_out_unit,*) "Error allocating x data structure"
|
|
goto 9999
|
|
end if
|
|
|
|
! Populate x class using data from x_global vector
|
|
call psb_scatter(x_global,x,desc_a,info,root=psb_root_)
|
|
if(info /= psb_success_) then
|
|
write(psb_out_unit,*) "Error in psb_scatter to populate x data structure"
|
|
goto 9999
|
|
end if
|
|
|
|
|
|
call psb_geall(y,desc_a,info)
|
|
if(info /= psb_success_) then
|
|
write(psb_out_unit,*) "Error allocating y data structure"
|
|
goto 9999
|
|
end if
|
|
|
|
! Populate y class using data from y_global vector
|
|
call psb_scatter(y_global,y,desc_a,info,root=psb_root_)
|
|
if(info /= psb_success_) then
|
|
write(psb_out_unit,*) "Error in psb_scatter to populate y data structure"
|
|
goto 9999
|
|
end if
|
|
|
|
|
|
! y = alpha * A * x + betha * y
|
|
call psb_spmm(sone,a,x,sone,y,desc_a,info)
|
|
if(info /= psb_success_) then
|
|
write(psb_out_unit,*) "Error in psb_spmm routine"
|
|
goto 9999
|
|
end if
|
|
|
|
! Save result to output file
|
|
call mm_array_write(y,"Result vector",info,filename="sol_m1_x1_y1.mtx")
|
|
|
|
! Deallocate
|
|
call psb_gefree(x, desc_a,info)
|
|
if(info /= psb_success_) then
|
|
write(psb_out_unit,*) "Error in vector x free routine"
|
|
goto 9999
|
|
end if
|
|
|
|
call psb_gefree(y, desc_a,info)
|
|
if(info /= psb_success_) then
|
|
write(psb_out_unit,*) "Error in vector y free routine"
|
|
goto 9999
|
|
end if
|
|
|
|
call psb_spfree(a, desc_a,info)
|
|
if(info /= psb_success_) then
|
|
write(psb_out_unit,*) "Error in matrix A free routine"
|
|
goto 9999
|
|
end if
|
|
|
|
call psb_cdfree(desc_a,info)
|
|
if(info /= psb_success_) then
|
|
write(psb_out_unit,*) "Error in matrix descriptor free routine"
|
|
goto 9999
|
|
end if
|
|
|
|
deallocate(x_global)
|
|
deallocate(y_global)
|
|
|
|
call psb_exit(ctxt)
|
|
stop
|
|
|
|
|
|
! Error handling
|
|
9999 call psb_error(ctxt)
|
|
|
|
call psb_errpush(info,name)
|
|
call psb_error_handler(ctxt,err_act)
|
|
call psb_exit(ctxt)
|
|
stop
|
|
|
|
|
|
contains
|
|
|
|
!> @brief Function to randomly generate x and y vectors
|
|
!! and save them on multiple files based on their
|
|
!! coefficients values.
|
|
!!
|
|
!! @param n The size of the vector.
|
|
subroutine generate_vectors(n)
|
|
use psb_base_mod
|
|
|
|
implicit none
|
|
|
|
integer(psb_ipk_), intent(in) :: n
|
|
integer(psb_ipk_), allocatable :: x(:), y(:)
|
|
real(psb_spk_), allocatable :: v(:)
|
|
integer(psb_ipk_) :: i
|
|
|
|
allocate(x(n))
|
|
allocate(y(n))
|
|
allocate(v(2*n))
|
|
|
|
|
|
call random_init(repeatable=.true.,image_distinct=.true.)
|
|
call random_number(v)
|
|
|
|
do i = 1,n
|
|
x(i) = int(v(i) * 100)
|
|
end do
|
|
|
|
do i = 1,n
|
|
y(i) = int(v(i+n) * 100)
|
|
end do
|
|
|
|
|
|
! Write only positive in x_1
|
|
call mm_array_write(x,"Positive vector",info,filename="vectors/x1.mtx")
|
|
call mm_array_write(y,"Positive vector",info,filename="vectors/y1.mtx")
|
|
|
|
! Write only negative in x_2
|
|
do i=1,n
|
|
x(i) = -x(i)
|
|
y(i) = -y(i)
|
|
end do
|
|
|
|
call mm_array_write(x,"Negative vector",info,filename="vectors/x2.mtx")
|
|
call mm_array_write(y,"Negative vector",info,filename="vectors/y2.mtx")
|
|
|
|
|
|
! Since numbers are less than one and always positive, we have to generate negative ones subtractiong 50
|
|
do i=1,n
|
|
x(i) = -x(i) ! Make the values positive again
|
|
x(i) = x(i) - 50
|
|
y(i) = -y(i) ! Make the values positive again
|
|
y(i) = y(i) - 50
|
|
end do
|
|
|
|
! Write random in x_3
|
|
call mm_array_write(x,"Random vector",info,filename="vectors/x3.mtx")
|
|
call mm_array_write(y,"Random vector",info,filename="vectors/y3.mtx")
|
|
|
|
! Write zero in x_4
|
|
do i=1,n
|
|
x(i) = 0
|
|
y(i) = 0
|
|
end do
|
|
|
|
call mm_array_write(x,"Null vector",info,filename="vectors/x4.mtx")
|
|
call mm_array_write(y,"Null vector",info,filename="vectors/y4.mtx")
|
|
|
|
deallocate(x)
|
|
deallocate(y)
|
|
deallocate(v)
|
|
|
|
end subroutine
|
|
|
|
end program
|
|
|