assign6-grade -- e. torres davila

       comments:
                 o the program compiles using the supplied makefile,
                   and runs

                 o but the makefile is not commented, is supposed to
                   be commented

                 o code for assessing danger of division by zero is
                   addressed in a coherent manner; avoids floating
                   point overflows

                 o the program seems to produce results that are
                   consistent w/ those produced by the c-lang
                   templates for this benchmark

                 o presents a clear, coherent user interface

                 o good, tight coding logic, good commenting, sharp
                   looking source code.

         grade: 97/100


! matvec.F90 -- Eduardo Torres Davila, 11/05/19
! Description -- Perform matrix-vector muliplication y = Ax, for a
! n x n matrix A, represented in double precision. We will prompt the
! user to provide n. We will then perform the matrix-vector product 3 ways:
! programmed row-based, programmed column-based, and by way of an external
! library call dgemv(). When calling each one we will use fortran built
! in timing function cpu_time() to get the run time of each
! matrix-vector product. We will then provide the timings for each and
! giving the ratio between row-based/column-based and ratio for
! column-based/dgemv(). We finally find the infinity norm of yB - yR and
! yB - yC to check that the results are accurate.

      program matvec
! Don't allow implicit type assignment for our variables
         implicit none

! Declaring variables
         double precision, allocatable      :: A(:,:), x(:), errvect(:)
         double precision, allocatable      :: yC(:), yR(:), yB(:)
         double precision                   :: rowerr, colerr, infnorm
         integer                            :: n
         integer                            :: ii, istatus
         real                               :: start, finish
         real                               :: t_R, t_C, t_B

! Prompt the user to give the dimensions of matrix A
         print '("Please provide n which is the size of matrix A.")'
         print '("A = n x n matrix.")'
         read (*,*) n

         print *, " "

! Allocate memory for matrix A
         allocate(A(n,n),STAT=istatus)
         if (istatus .NE. 0) then
            print *, "failed to allocate storage for matrix A. exiting."
            stop
         end if

! Allocate memory for vector x
         allocate(x(n),STAT=istatus)
         if (istatus .NE. 0) then
            print *, "failed to allocate storage for vector x. exiting."
            stop
         end if

! Allocate memory for vector yR
         allocate(yR(n),STAT=istatus)
         if (istatus .NE. 0) then
            print *, "failed to allocate storage for vector y. exiting."
            stop
         end if

! Allocate memory for vector yC
         allocate(yC(n),STAT=istatus)
         if (istatus .NE. 0) then
            print *, "failed to allocate storage for vector y. exiting."
            stop
         end if

! Allocate memory for vector yB
         allocate(yB(n),STAT=istatus)
         if (istatus .NE. 0) then
            print *, "failed to allocate storage for vector y. exiting."
            stop
         end if

! Allocate memory for error vector errvect
         allocate(errvect(n),STAT=istatus)
         if (istatus .NE. 0) then
            print *, "failed to allocate storage for errvect. exiting."
            stop
         end if

! Load vector x and A with values
         call loadvalues(n, x, A)

! Calculate vector y = Ax doing it Row based and timing the computation
         call cpu_time(start)
         call matvecR(n, A, x, yR)
         call cpu_time(finish)

! Print the duration of the computation and display time to STDOUT
         t_R = finish - start
         print '("Time for (row-oriented) matvecR() = ",f10.6,"&
            & seconds.")', t_R

! Calculate vector y = Ax doing it Column based and timing the operation
         call cpu_time(start)
         call matvecC(n, A, x, yC)
         call cpu_time(finish)

! Print the duration of the computation and display time to STDOUT
         t_C = finish - start
         print '("Time for (column-oriented) matvecC() = ",f10.6,"&
            & seconds.")', t_C

! Usage of dgemv() --
! first entry - 'n' means y = \a A x + \b y
! second and third - number of rows and num of cols
! 4th entry  - the scalar \a
! 5th entry  - the matrix A
! 6th entry  - first dimension of A
! 7th entry  - our x vector - x
! 8th entry  - our x increment - 1
! 9th entry  - the scalar \b
! 10th entry - our y vector - yB
! 11th entry - our y increment - 1
         call cpu_time(start)
         call dgemv('N', n, n, 1.0d+0, A, n, x, 1, 0.0d+0, yB, 1)
         call cpu_time(finish)

! Print the duration of the computation and display time to STDOUT
         t_B = finish - start
         print '("Time for dgemv() = ",f10.6," seconds.")', t_B

         print *, " "

! Before computing ratio we make sure we are not dividing by zero
         if ( abs(t_C) .LT. 0.000001 ) then
            print '("WARNING: looks like we have a zero timing result&
               & for col-oriented computation,")'
            print '("WARNING: maybe the problem size is too&
               & fine-grained...")'
         else
            print '("Ratio of R timing...")'
            print '("row oriented / col oriented = ",f10.7,".")', t_R/t_C
         end if

         print *, " "

! Before computing ratio we make sure we are not dividing by zero
         if ( abs(t_B) .LT. 0.000001 ) then
            print '("WARNING: looks like we have a zero timing result&
               & for dgemv() computation,")'
            print '("WARNING: maybe the problem size is too &
               &fine-grained...")'
         else
            print '("Ratio of C timing...")'
            print '("col oriented / call dgemv() = ",f10.7,".")', t_C/t_B
         end if

         print *, " "

! Generate error vector for row-based vector and check infinity norm for
! correctness of the computation.
         errvect = yB - yR
         rowerr = infnorm(n, errvect)
         print '("Norm of row-based matvec error = ",f6.4,".")', rowerr

! Generate error vector for column-based vector and check infinity norm
! for correctness of the computation.
         errvect = yB - yC
         colerr = infnorm(n, errvect)
         print '("Norm of column-based matvec error = ",f6.4,".")',&
            &colerr

! deallocate all memory for matrices/vectors created
         deallocate(A)
         deallocate(x)
         deallocate(yR)
         deallocate(yC)
         deallocate(yB)
         deallocate(errvect)
      end program matvec

! loadvalues() -- this routine loads a vector and matrix of dimensions
!                 n with specific values.
!
! param list:
!         n -- integer value -- input value, the number of array elements
!         x -- empty array of size n with double precision elements
!              output value
!         A -- empty 2 dimensional array with double precision elements
!              which is our n x n matrix -- output value

      subroutine loadvalues(n, x, A)
! Don't allow implicit type assignment for our variables
         implicit none

         integer, intent(in)                             :: n
         integer                                         :: ii, jj
         double precision, dimension(n,n), intent(out)   :: A
         double precision, dimension(n), intent(out)     :: x

         do ii=1, n
            x(ii) = 1.0
            do jj=1, n
               A(ii,jj) = (ii+1)*(jj+1)/(float(n))
            end do
         end do

      end subroutine loadvalues

! matvecR() -- row-oriented matrix-vector multiplication routine.
!
! param list:
!         n -- integer value -- input value, the number of array elements
!         A -- 2 dimensional array filled with double precision elements
!              which is our n x n matrix -- input value
!         x -- array with n double precision values -- input value
!         y -- array with n double precision values -- output value

      subroutine matvecR(n, A, x, y)
! Don't allow implicit type assignment for our variables
         implicit none

         integer, intent(in)             :: n
         integer                         :: ii, jj
         double precision, intent(in)    :: A(n,n), x(n)
         double precision, intent(out)   :: y(n)

         do ii=1, n
            y(ii) = a(ii,1) * x(1)
            do jj=2, n
               y(ii) = y(ii) + a(ii,jj) * x(jj)
            end do
         end do

      end subroutine matvecR

! matvecC() -- column-oriented matrix-vector multiplication routine.
!
! param list:
!         n -- integer value -- input value, the number of array elements
!         A -- 2 dimensional array filled with double precision elements
!              which is our n x n matrix -- input value
!         x -- array with n double precision values -- input value
!         y -- array with n double precision values -- output value

      subroutine matvecC(n, A, x, y)
! Don't allow implicit type assignment for our variables
         implicit none

         integer, intent(in)             :: n
         integer                         :: ii, jj
         double precision, intent(in)    :: A(n,n), x(n)
         double precision, intent(out)   :: y(n)

         do ii=1, n
            y(ii) = a(ii,1) * x(1)
         end do

         do jj=2, n
            do ii=1, n
               y(ii) = y(ii) + a(ii,jj) * x(jj)
            end do
         end do

      end subroutine matvecC

! function body for infnorm which is used to find the maximal value of
! the absolute value of the array xx.
!
! param list:
!         n -- integer value -- input value, the length of the array xx
!         xx -- array with n double precision values -- input value
! return value:
!         infnorm -- double precision -- return the absolute value of
!         the maximal double precision value in array xx

      double precision function infnorm( n, xx )
! Don't allow implicit type assignment for our variables
         implicit none

         integer, intent(in)   :: n
         double precision, intent(in)   :: xx(n)

         infnorm = maxval(dabs(xx))
      end function