一、程序填空题

程序填空第1题
1、用OpenMP Fortran API并行化下面的程序：

! This program implements a simple molecular dynamics simulation,!
   using the velocity Verlet time integration scheme. The particles
!   interact with a central pair potential.
!
! Author:   Bill Magro, Kuck and Associates, Inc. (KAI), 1998
!
! Parallelism is implemented via OpenMP directives.
! THIS PROGRAM USES THE FORTRAN90 RANDOM_NUMBER FUNCTION AND ARRAY
!   SYNTAX
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
      program md
      implicit none
      ! simulation parameters
      integer ndim       ! dimensionality of the physical space
      integer nparts     ! number of particles
      integer nsteps     ! number of time steps in the simulation
      parameter(ndim=3,nparts=500,nsteps=1000)
      real*8 mass        ! mass of the particles
      real*8 dt          ! time step
      real*8 box(ndim)   ! dimensions of the simulation box
      parameter(mass=1.0,dt=1.0e-4)
      ! simulation variables
      real*8 position(ndim,nparts)
      real*8 velocity(ndim,nparts)
      real*8 force(ndim,nparts)
      real*8 accel(ndim,nparts)
      real*8 potential, kinetic, E0
      integer i
      box(1:ndim) = 10.
      ! set initial positions, velocities, and accelerations
      call initialize(nparts,ndim,box,position,velocity,accel)
      print *,"initialize successful"
      ! compute the forces and energies
      call compute(nparts,ndim,box,position,velocity,mass,force,potential,kinetic)
      print *,"compute successful"
      E0 = potential + kinetic
      ! This is the main time stepping loop
      do i=1,nsteps
          call compute(nparts,ndim,box,position,velocity,mass,force,potential,kinetic)
          write(*,*) potential, kinetic,(potential + kinetic - E0)/E0
          call update(nparts,ndim,position,velocity,force,accel,mass,dt)
      enddo
      end
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! Compute the forces and energies, given positions, masses,
! and velocities
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
      subroutine compute(np,nd,box,pos,vel,mass,f,pot,kin)
      implicit none
      integer np
      integer nd
      real*8  box(nd)
      real*8  pos(nd,np)
      real*8  vel(nd,np)
      real*8  f(nd,np)
      real*8  mass
      real*8  pot
      real*8  kin
      real*8 dotr8
      external dotr8
      real*8 v, dv, x
      integer i, j, k
      real*8  rij(nd)
      real*8  d
      real*8  PI2
      parameter(PI2=3.14159265d0/2.0d0)
      ! statement function for the pair potential and its derivative
      ! This potential is a harmonic well which smoothly saturates to a
      ! maximum value at PI/2.
      v(x) = sin(min(x,PI2))**2.
      dv(x) = 2.*sin(min(x,PI2))*cos(min(x,PI2))
      pot = 0.0
      kin = 0.0
      ! The computation of forces and energies is fully parallel.
     
     
     
     
      do i=1,np
        ! compute potential energy and forces
        f(1:nd,i) = 0.0
        do j=1,np
             if (i .ne. j) then
               call dist(nd,box,pos(1,i),pos(1,j),rij,d)
               ! attribute half of the potential energy to particle 'j'
               pot = pot + 0.5*v(d)
               do k=1,nd
                 f(k,i) = f(k,i) - rij(k)*dv(d)/d
               enddo
             endif
        enddo
        ! compute kinetic energy
        kin = kin + dotr8(nd,vel(1,i),vel(1,i))
      enddo
     
      kin = kin*0.5*mass
      return
      end
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! Initialize the positions, velocities, and accelerations.
! The Fortran90 random_number function is used to choose positions.
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
      subroutine initialize(np,nd,box,pos,vel,acc)
      implicit none
      integer np
      integer nd
      real*8  box(nd)
      real*8  pos(nd,np)
      real*8  vel(nd,np)
      real*8  acc(nd,np)
      integer i, j
      real*8 x
      do i=1,np
        do j=1,nd
          call random_number(x)
          pos(j,i) = box(j)*x
          vel(j,i) = 0.0
          acc(j,i) = 0.0
        enddo
      enddo
      return
      end
! Compute the displacement vector (and its norm) between two particles.
      subroutine dist(nd,box,r1,r2,dr,d)
      implicit none
      integer nd
      real*8 box(nd)
      real*8 r1(nd)
      real*8 r2(nd)
      real*8 dr(nd)
      real*8 d
      integer i
      d = 0.0
      do i=1,nd
        dr(i) = r1(i) - r2(i)
        d = d + dr(i)**2.
      enddo
      d = sqrt(d)
      return
      end
! Return the dot product between two vectors of type real*8 and length n
      real*8 function dotr8(n,x,y)
      implicit none
      integer n
      real*8 x(n)
      real*8 y(n)
      integer i
      dotr8 = 0.0
      do i = 1,n
        dotr8 = dotr8 + x(i)*y(i)
      enddo
      return
      end
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! Perform the time integration, using a velocity Verlet algorithm
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
      subroutine update(np,nd,pos,vel,f,a,mass,dt)
      implicit none
      integer np
      integer nd
      real*8  pos(nd,np)
      real*8  vel(nd,np)
      real*8  f(nd,np)
      real*8  a(nd,np)
      real*8  mass
      real*8  dt
      integer i, j
      real*8  rmass
      rmass = 1.0/mass
      ! The time integration is fully parallel
     
     
     
      do i = 1,np
        do j = 1,nd
          pos(j,i) = pos(j,i) + vel(j,i)*dt + 0.5*dt*dt*a(j,i)
          vel(j,i) = vel(j,i) + 0.5*dt*(f(j,i)*rmass + a(j,i))
          a(j,i) = f(j,i)*rmass
        enddo
      enddo
     
      return
      end

2、用OpenMP Fortran API 并行化下面的程序。

      program main
************************************************************
* program to solve a finite difference
* discretization of Helmholtz equation : 
* (d2/dx2)u + (d2/dy2)u - alpha u = f
* using Jacobi iterative method.
*
* Modified: Sanjiv Shah,       Kuck and Associates, Inc. (KAI), 1998
* Author:   Joseph Robicheaux, Kuck and Associates, Inc. (KAI), 1998
*
* Directives are used in this code to achieve paralleism.
* All do loops are parallized with default 'static' scheduling.
*
* Input :  n - grid dimension in x direction
*          m - grid dimension in y direction
*          alpha - Helmholtz constant (always greater than 0.0)
*          tol   - error tolerance for iterative solver
*          relax - Successice over relaxation parameter
*          mits  - Maximum iterations for iterative solver
*
* On output
*       : u(n,m) - Dependent variable (solutions)
*       : f(n,m) - Right hand side function
*************************************************************
      implicit none
      integer n,m,mits
      double precision tol,relax,alpha
      common /idat/ n,m,mits
      common /fdat/tol,alpha,relax
*
* Read info
*
      write(*,*) "Input n,m - grid dimension in x,y direction "
      read(5,*) n,m
      write(*,*) "Input alpha - Helmholts constant "
      read(5,*) alpha
      write(*,*) "Input relax - Successive over-relaxation parameter"
      read(5,*) relax
      write(*,*) "Input tol - error tolerance for iterative solver"
      read(5,*) tol
      write(*,*) "Input mits - Maximum iterations for solver"
      read(5,*) mits
*
* Calls a driver routine
*
      call driver ()
      stop
      end
      subroutine driver ( )
*************************************************************
* Subroutine driver ()
* This is where the arrays are allocated and initialzed.
*
* Working varaibles/arrays
*     dx  - grid spacing in x direction
*     dy  - grid spacing in y direction
*************************************************************
      implicit none
      integer n,m,mits,mtemp
      double precision tol,relax,alpha
      common /idat/ n,m,mits,mtemp
      common /fdat/tol,alpha,relax
      double precision u(n,m),f(n,m),dx,dy
* Initialize data
      call initialize (n,m,alpha,dx,dy,u,f)
* Solve Helmholtz equation
      call jacobi (n,m,dx,dy,alpha,relax,u,f,tol,mits)
* Check error between exact solution
      call  error_check (n,m,alpha,dx,dy,u,f)
      return
      end
      subroutine initialize (n,m,alpha,dx,dy,u,f)
******************************************************
* Initializes data
* Assumes exact solution is u(x,y) = (1-x^2)*(1-y^2)
*
******************************************************
      implicit none
      integer n,m
      double precision u(n,m),f(n,m),dx,dy,alpha
      integer i,j, xx,yy
      double precision PI
      parameter (PI=3.1415926)
      dx = 2.0 / (n-1)
      dy = 2.0 / (m-1)
* Initilize initial condition and RHS
     
      do j = 1,m
         do i = 1,n
            xx = -1.0 + dx * dble(i-1)        ! -1 < x < 1
            yy = -1.0 + dy * dble(j-1)        ! -1 < y < 1
            u(i,j) = 0.0
            f(i,j) = -alpha *(1.0-xx*xx)*(1.0-yy*yy) - 2.0*(1.0-xx*xx)-2.0*(1.0-yy*yy)
         enddo
      enddo
     
      return
      end
      subroutine jacobi (n,m,dx,dy,alpha,omega,u,f,tol,maxit)
******************************************************************
* Subroutine HelmholtzJ
* Solves poisson equation on rectangular grid assuming :
* (1) Uniform discretization in each direction, and
* (2) Dirichlect boundary conditions
*
* Jacobi method is used in this routine
*
* Input : n,m   Number of grid points in the X/Y directions
*         dx,dy Grid spacing in the X/Y directions
*         alpha Helmholtz eqn. coefficient
*         omega Relaxation factor
*         f(n,m) Right hand side function
*         u(n,m) Dependent variable/Solution
*         tol    Tolerance for iterative solver
*         maxit  Maximum number of iterations
*
* Output : u(n,m) - Solution
*****************************************************************
      implicit none
      integer n,m,maxit
      double precision dx,dy,f(n,m),u(n,m),alpha, tol,omega
*
* Local variables
*
      integer i,j,k,k_local
      double precision error,resid,rsum,ax,ay,b
      double precision error_local, uold(n,m)
      real ta,tb,tc,td,te,ta1,ta2,tb1,tb2,tc1,tc2,td1,td2
      real te1,te2
      real second
      external second
*
* Initialize coefficients
      ax = 1.0/(dx*dx) ! X-direction coef
      ay = 1.0/(dy*dy) ! Y-direction coef
      b  = -2.0/(dx*dx)-2.0/(dy*dy) - alpha ! Central coeff 
      error = 10.0 * tol
      k = 1
      do while (k.le.maxit .and. error.gt. tol)
         error = 0.0   
* Copy new solution into old
     
     
         do j=1,m
            do i=1,n
               uold(i,j) = u(i,j)
            enddo
         enddo
* Compute stencil, residual, & update
     
         do j = 2,m-1
            do i = 2,n-1
*     Evaluate residual
               resid = (ax*(uold(i-1,j) + uold(i+1,j))& +  ay*(uold(i,j-1) + uold(i,j+1))& + b * uold(i,j) - f(i,j))/b
* Update solution
               u(i,j) = uold(i,j) - omega * resid
* Accumulate residual error
               error = error + resid*resid
            end do
         enddo
     
     
* Error check
         k = k + 1
         error = sqrt(error)/dble(n*m)
*
      enddo                     ! End iteration loop
*
      print *, 'Total Number of Iterations ', k
      print *, 'Residual                   ', error
      return
      end
      subroutine error_check (n,m,alpha,dx,dy,u,f)
      implicit none
************************************************************
* Checks error between numerical and exact solution
*
************************************************************
      integer n,m
      double precision u(n,m),f(n,m),dx,dy,alpha
      integer i,j
      double precision xx,yy,temp,error
      dx = 2.0 / (n-1)
      dy = 2.0 / (m-1)
      error = 0.0
     
      do j = 1,m
         do i = 1,n
            xx = -1.0d0 + dx * dble(i-1)
            yy = -1.0d0 + dy * dble(j-1)
            temp  = u(i,j) - (1.0-xx*xx)*(1.0-yy*yy)
            error = error + temp*temp
         enddo
      enddo
      error = sqrt(error)/dble(n*m)
      print *, 'Solution Error : ',error
      return

      end