diff --git a/test/computational_routines/geaxpby/geaxpby.f90 b/test/computational_routines/geaxpby/geaxpby.f90
index 3786368b..51a58358 100644
--- a/test/computational_routines/geaxpby/geaxpby.f90
+++ b/test/computational_routines/geaxpby/geaxpby.f90
@@ -4,13 +4,39 @@ program main
 
     implicit none
 
-
     ! MPI variables
     integer(psb_ipk_)       :: my_rank, np
 
     ! Communicator variable
     type(psb_ctxt_type)     :: ctxt
 
+    ! parameters array
+    character(len=64)       :: x(4),y(4)  
+    real(psb_ipk_)          :: alpha(3), beta(3)  
+
+    ! cycle indexes variables
+    integer(psb_ipk_)       :: i,j,k,h,l
+
+
+    ! 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"
+
+    alpha(1) = sone
+    alpha(2) = -sone
+    alpha(3) = szero
+
+    beta(1) = sone
+    beta(2) = -sone
+    beta(3) = szero
+
 
     call psb_init(ctxt)
     call psb_info(ctxt,my_rank,np)
@@ -18,10 +44,25 @@ program main
     if(my_rank == psb_root_) then
         write(psb_out_unit,*) 'Welcome to PSBLAS version: ',psb_version_string_
         write(psb_out_unit,*) 'This is the psb_geaxpby_test sample program'
+
+        call generate_vectors(10000,10000) 
     end if
 
     call psb_barrier(ctxt)
 
+    if(my_rank == psb_root_) then
+        write(psb_out_unit,*) size(x)
+    end if
+
+    do i=1,size(x)
+        do j=1,size(y)
+            do k=1,size(alpha)
+                do h=1,size(beta)
+                    call psb_geaxpby_kernel(x_file=x(i), y_file=y(j), alpha = alpha(k), beta = beta(h), ctxt = ctxt)
+                end do
+            end do
+        end do
+    end do
 
 
 
diff --git a/test/computational_routines/geaxpby/psb_geaxpby_test b/test/computational_routines/geaxpby/psb_geaxpby_test
deleted file mode 100755
index 420ac3ed..00000000
Binary files a/test/computational_routines/geaxpby/psb_geaxpby_test and /dev/null differ
diff --git a/test/computational_routines/geaxpby/psb_geaxpby_test.f90 b/test/computational_routines/geaxpby/psb_geaxpby_test.f90
index ee8e6fdb..30b8e191 100644
--- a/test/computational_routines/geaxpby/psb_geaxpby_test.f90
+++ b/test/computational_routines/geaxpby/psb_geaxpby_test.f90
@@ -47,8 +47,224 @@
 !!          Specified as: An integer value; 0 means no error has been detected.
 !!
 module psb_geaxpby_test
+
     contains
-    subroutine psb_spmm_kernel()
+    
+    subroutine psb_geaxpby_kernel(x_file, y_file, alpha, beta, ctxt)
+        use psb_base_mod
+        use psb_util_mod
+
+        implicit none 
+
+        ! input parameters
+        character(len = *), intent(in)      :: x_file, y_file
+        real(psb_spk_), intent(in)          :: alpha, beta
+
+        character(len=:), allocatable       :: output_file_name       
+
+        ! vectors
+        type(psb_s_vect_type)           	:: x, y
+
+        ! matrix descriptor data structure
+        type(psb_desc_type)             	:: desc_a
+
+        ! communication context
+        type(psb_ctxt_type), intent(in)     :: ctxt
+        integer(psb_ipk_)               	:: my_rank, np, info, err_act
+
+        ! variables outside PSLBALS data structures
+        real(psb_spk_), allocatable     	:: x_global(:), y_global(:)
+        integer(psb_ipk_)               	:: i
+
+        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
+            write(*,*) "Here"
+            allocate(x_global(10000))
+            allocate(y_global(10000))
+            call mm_array_read(x_global,info,filename=x_file)
+            call mm_array_read(y_global,info,filename=y_file)
+        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 * x + beta * y
+        call psb_geaxpby(alpha,x,beta,y,desc_a,info)
+        if(info /= psb_success_) then
+            write(psb_out_unit,*) "Error in psb_spmm routine"
+            goto 9999
+        end if
+
+        ! Make the root process be the one that saves everything on file
+        if(np == 1) then 
+            output_file_name = "serial/"
+        else 
+            output_file_name = "parallel/"
+        end if
+
+        output_file_name = output_file_name // "sol_" // x_file(9:10) // "_" // y_file(9:10)
+
+        if(alpha == sone) then 
+            output_file_name = output_file_name // "_a1"
+        else if(alpha == -sone) then 
+            output_file_name = output_file_name // "_a2"
+        else if(alpha == szero) then
+            output_file_name = output_file_name // "_a3"
+        end if
+
+        if(beta == sone) then 
+            output_file_name = output_file_name // "_b1.mtx"
+        else if(beta == -sone) then 
+            output_file_name = output_file_name // "_b2.mtx"
+        else if(beta == szero) then
+            output_file_name = output_file_name // "_b3.mtx"
+        end if
+
+
+        ! Save result to output file
+        call mm_array_write(y,"Result vector",info,filename=output_file_name)
+
+        ! 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_cdfree(desc_a,info)
+        if(info /= psb_success_) then
+            write(psb_out_unit,*) "Error in matrix descriptor free routine"
+            goto 9999
+        end if
+
+        if(my_rank == 0) then
+            deallocate(x_global)
+            deallocate(y_global)
+        end if
+
+        return
+
+
+        ! Error handling
+        9999 call psb_error(ctxt)
+
+        call psb_error_handler(ctxt,err_act)
+        stop
+
+    end subroutine
+
+
+
+
+        !> @brief Function to randomly generate x and y vectors 
+    !!        and save them on multiple files based on their
+    !!        coefficients values.
+    !!
+    subroutine generate_vectors(rows, cols)
+        use psb_base_mod
+        use psb_util_mod
+
+        implicit none 
+
+        integer(psb_ipk_), intent(in)               :: rows, cols
+        real(psb_spk_), allocatable                 :: x(:), y(:)
+        integer(psb_ipk_)                           :: i, info
+
+        allocate(x(rows))
+        allocate(y(cols))
+        
+        call random_init(repeatable=.true.,image_distinct=.true.) 
+        call random_number(x)
+        call random_number(y)
+
+        ! 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,rows 
+            x(i) = -x(i)
+        end do  
+
+        do i=1,cols
+            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,rows
+            x(i) = -x(i) ! Make the values positive again  
+            x(i) = x(i) - 0.5
+        end do   
+
+        do i=1,cols
+            y(i) = -y(i) ! Make the values positive again  
+            y(i) = y(i) - 0.5
+        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,rows
+            x(i) = 0
+        end do   
+
+        do i=1,cols
+            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)
+
     end subroutine
 
 end module psb_geaxpby_test
\ No newline at end of file
diff --git a/test/computational_routines/spmm/psb_spmm_test.f90 b/test/computational_routines/spmm/psb_spmm_test.f90
index 35f749e9..938a113b 100644
--- a/test/computational_routines/spmm/psb_spmm_test.f90
+++ b/test/computational_routines/spmm/psb_spmm_test.f90
@@ -36,13 +36,12 @@ module psb_spmm_test
         integer(psb_ipk_)               	:: rows, cols, nnz
         integer(psb_ipk_)               	:: nr, nt ! In BLOCK ROWS distributin, the number of rows 
 
-
+        ! variables outside PSLBALS data structures
         real(psb_spk_), allocatable     	:: x_global(:), y_global(:)
         integer(psb_ipk_)               	:: i
 
         info = psb_success_
 
-
         call psb_info(ctxt,my_rank,np)
 
         if (my_rank < 0) then 
diff --git a/test/computational_routines/spmm/spmm.f90 b/test/computational_routines/spmm/spmm.f90
index 07344a0f..88f2a748 100644
--- a/test/computational_routines/spmm/spmm.f90
+++ b/test/computational_routines/spmm/spmm.f90
@@ -56,10 +56,10 @@ program main
 
 
     !! 1138_bus matrix (sparse)
-    do i=1,4 
-        do j=1,4
-            do k=1,3
-                do h=1,3
+    do i=1,size(x)
+        do j=1,size(y)
+            do k=1,size(alpha)
+                do h=1,size(beta)
                     call psb_spmm_kernel(mtx_file="matrix/1138_bus.mtx",x_file=x(i), y_file=y(j), &
                         & alpha = alpha(k), beta = beta(h), ctxt = ctxt)
                 end do