Advertisement
Guest User

fortran matrix vector_multiplication

a guest
Jan 3rd, 2014
239
0
Never
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
  1. program mvmul_test
  2.     implicit none;
  3.     integer, parameter :: m=2000,n=1000;
  4.     integer :: i,j;
  5.     integer, parameter :: iterations = 15;
  6.     integer, parameter :: kfloat = kind(0.D0);
  7.     real(kind=kfloat), dimension(1:n,1:m) :: A;
  8.     real(kind=kfloat), dimension(1:m)     :: x;
  9.     real(kind=kfloat), dimension(1:n)     :: y1, y2;
  10.     real(kind=kfloat) :: t1, t2, t_naive, t_optimized;
  11.    
  12.     do i = 1,m
  13.         x(i) = 1.0
  14.         do j = 1,n
  15.             A(i,j) = (-1.0)**(i-j);
  16.         end do
  17.     end do
  18.    
  19.     call cpu_time(t1);
  20.     do i = 1,iterations
  21.     call naive(A, m, n, x, y1);
  22.     end do
  23.     call cpu_time(t2);
  24.     t_naive = (t2-t1)/iterations;
  25.    
  26.     call cpu_time(t1);
  27.     do i = 1,iterations
  28.     call optimized(A, m, n, x, y2);
  29.     end do
  30.     call cpu_time(t2);
  31.     t_optimized = t2-t1/iterations;
  32.    
  33.     if( all(y1 == y2) ) then
  34.         write(*,*) "Vectors match! Calculations are OK.";
  35.     endif
  36.    
  37.     write(*,*) "Optimized time:", t_optimized
  38.     write(*,*) "Naive time:", t_naive;
  39.     write(*,*) "Ratio (t_optimized/t_naive):", t_optimized/t_naive
  40.    
  41.    
  42.    
  43.     !-------------------------------------------------------------------------
  44.     contains
  45.    
  46.     subroutine naive(A, m, n, x, y)
  47.         implicit none;
  48.         real(kind=kfloat), dimension(1:n,1:m), intent(in) :: A;
  49.         real(kind=kfloat), dimension(1:m), intent(in) :: x;
  50.         real(kind=kfloat), dimension(1:n), intent(out) :: y;
  51.         integer, intent(in) :: m, n;
  52.         integer :: i,j;
  53.        
  54.         do i = 1,n
  55.             y(i) = 0.0;
  56.         end do
  57.        
  58.         do i = 1,n
  59.             do j = 1,m
  60.                 y(i) = y(i) + A(i,j)*x(j);
  61.             end do
  62.         end do
  63.        
  64.        
  65.     end subroutine
  66.    
  67.     subroutine optimized(A, m, n, x, y)
  68.         implicit none;
  69.         real(kind=kfloat), dimension(1:n,1:m), intent(in) :: A;
  70.         real(kind=kfloat), dimension(1:m), intent(in) :: x;
  71.         real(kind=kfloat), dimension(1:n), intent(out) :: y;
  72.         integer, intent(in) :: m, n;
  73.         integer :: i,j;
  74.        
  75.         do i = 1,n
  76.             y(i) = 0.0;
  77.         end do
  78.        
  79.         do j = 1,m
  80.             do i = 1,n
  81.                 y(i) = y(i) + A(i,j)*x(j);
  82.             end do
  83.         end do
  84.        
  85.        
  86.     end subroutine
  87. end program
Advertisement
Add Comment
Please, Sign In to add comment
Advertisement