# Commits

committed 73306ff Merge

• Participants
• Parent commits fe2f4b4, 544b6ff
• Branches coursera

# File codes/f2py/Makefile

+
+.PHONY: clean
+
+%.so : %.f90
+	f2py -m $* -c$<
+
+clean:
+	rm -f *.so *.pyc
+

# File codes/f2py/fcn1.f90

+! $UWHPSC/codes/f2py/fcn1.f90 + +! A test function for trying out f2py: +!$ f2py -m fcn1 -c fcn1.f90
+!   $python +! >>> import fcn1 +! >>> fcn1.f1(1.) +! 2.7182818284590451 + + +function f1(x) + real(kind=8), intent(in) :: x + real(kind=8) :: f1 + f1 = exp(x) +end function f1 # File codes/f2py/jacobi1.f90 + +!$UWHPSC/codes/f2py/jacobi1
+
+! For use with f2py:
+!   $f2py -m jacobi1 -c jacobi1.f90 +! Then in Python you can import jacobi1 and use it to take a fixed number of +! iterations using Jacobi in one space dimension to solve +! u(i-1) - 2*u(i) + u(i+1) = -dx**2 * f(i). +! Python call: +! >>> u = jacobi1.iterate(u0, iters, f) +! where  +! iters is the number of iterations to take, +! u0(0:n+1) is the initial guess (with desired boundary values  +! in u0(0) and u0(n+1), these are not changed. +! f(0:n+1) contains values of f(x) at each grid point. +!  + +subroutine iterate(u0,iters,f,u,n) + implicit none + integer, intent(in) :: n, iters + real(kind=8), dimension(0:n+1), intent(in) :: u0, f + real(kind=8), dimension(0:n+1), intent(out) :: u + + ! local variables: + real(kind=8), dimension(:), allocatable :: uold + real(kind=8) :: dx + integer :: i, iter + + ! allocate storage for boundary points too: + allocate(uold(0:n+1)) + + ! grid spacing: + dx = 1.d0 / (n+1.d0) + + uold = u0 + + ! initialize boundary values of u: + u(0) = u0(0) + u(n+1) = u0(n+1) + + ! take the iterations requested: + do iter=1,iters + do i=1,n + u(i) = 0.5d0*(uold(i-1) + uold(i+1) + dx**2*f(i)) + enddo + uold = u + enddo + +end subroutine iterate # File codes/f2py/plot_jacobi_iterates.py + +from matplotlib import pyplot as plt +import numpy as np +import time + +import jacobi1 # created from jacobi1.f90 using f2py + +n = 51 + +# grid points: +x = np.linspace(0., 1., n) + +# boundary conditions: +alpha = 20. +beta = 60. + +# initial guess: +u = np.linspace(alpha, beta, n) + +# right hand side: +f = 100. * np.exp(x) + +niters = 0 +plt.clf() +plt.plot(x, u, 'o-') +plt.title("After %s iterations" % niters) +plt.ylim([10., 80.]) + +ans = raw_input("Type return to start... ") + +iters_per_plot = 200 # number of iterations between plots +nplots = 20 # number of plots to produce + +for nn in range(nplots): + u = jacobi1.iterate(u, iters_per_plot, f) + plt.plot(x, u, 'o-') + niters = niters + iters_per_plot + plt.title("After %s iterations" % niters) + plt.ylim([10., 80.]) + plt.draw() + time.sleep(.5) +  + + # File codes/f2py/sub1.f90 + +!$UWHPSC/codes/f2py/sub1.f90
+
+! A test subroutine for trying out f2py:
+!   $f2py -m sub1 -c sub1.f90 +!$ python
+!   >>> import sub1
+!   >>> sub1.mysub(3., 5.)
+!   (8.0, -2.0)
+
+subroutine mysub(a,b,c,d)
+    real (kind=8), intent(in) :: a,b
+    real (kind=8), intent(out) :: c,d
+    c = a+b
+    d = a-b
+end subroutine mysub

# File codes/fortran/zeroin/Makefile

+# $CLASSHG/codes/fortran/zeroin/Makefile + +OBJECTS = zeroin.o myfcn.o testzeroin.o + +FFLAGS = -g + +.PHONY: test plot clean clobber + +test: testzeroin.exe + ./testzeroin.exe + +testzeroin.exe:$(OBJECTS)
+	gfortran $(OBJECTS) -o testzeroin.exe + +%.o : %.f + gfortran$(FFLAGS) -c $<  + +%.o : %.f90 + gfortran$(FFLAGS) -c $<  + +clean: + rm -f *.o *.exe *.mod + # File codes/fortran/zeroin/README.txt + +Test program for zeroin. +A zero of the function f(x) is computed in the interval ax,bx . + +This function came from Netlib: + http://www.netlib.org/go/zeroin.f + # File codes/fortran/zeroin/myfcn.f90 +!$CLASSHG/codes/fortran/zeroin/myfcn.f90
+
+! Provides the function power(x) = x**n - a.
+! A zero of this function will be found using zeroin.
+
+module myfcn
+
+    implicit none
+    real(kind=8) :: a
+    integer :: n
+    save
+
+contains
+
+real(kind=8) function power(x)
+
+    implicit none
+    real(kind=8) :: x
+    power = x**n - a
+
+end function power
+
+end module myfcn

+! $UWHPSC/codes/fortran/roots/testzeroin.f90 + +program main + + use myfcn + ! This module holds n, a, and provides function power(x) + + implicit none + real(kind=8) :: x, x1, x2, y, tol + real(kind=8), external :: zeroin + + print *, "Test routine for computing n'th root of a" + print *, "Input a " + read *, a + print *, "Input n " + read *, n + print *, "Input x1,x2: interval containing zero" + read *, x1, x2 + + tol = 1.d-10 ! tolerance requested to zeroin + + y = zeroin(x1, x2, power, tol) + + print *, 'Estimated value: ', y + print *, 'Correct value: ', a ** (1.d0/n) + +end program main # File codes/fortran/zeroin/zeroin.f +c$UWHPSC/codes/fortran/zeroin/zeroin.f
+c     This function came from Netlib:
+c        http://www.netlib.org/go/zeroin.f
+c     With some modification to remove dependence on d1mach
+c     by computing eps, the machine precision, in this routine.
+c
+c
+c     =================================================
+      function zeroin(ax,bx,f,tol)
+c     =================================================
+      implicit double precision (a-h,o-z)
+      external f
+c
+c      a zero of the function  f(x)  is computed in the interval ax,bx .
+c      (Standard routine from netlib)
+c
+c  input..
+c
+c  ax     left endpoint of initial interval
+c  bx     right endpoint of initial interval
+c  f      function subprogram which evaluates f(x) for any x in
+c         the interval  ax,bx
+c  tol    desired length of the interval of uncertainty of the
+c         final result ( .ge. 0.d0)
+c
+c
+c  output..
+c
+c  zeroin abcissa approximating a zero of  f  in the interval ax,bx
+c
+c
+c      it is assumed  that   f(ax)   and   f(bx)   have  opposite  signs
+c  without  a  check.  zeroin  returns a zero  x  in the given interval
+c  ax,bx  to within a tolerance  4*macheps*dabs(x) + tol, where macheps
+c  is the relative machine precision.
+c
+c      this function subprogram is a slightly  modified  translation  of
+c  the algol 60 procedure  zero  given in  richard brent, algorithms for
+c  minimization without derivatives, prentice - hall, inc. (1973).
+c
+c
+c
+c  compute eps, the relative machine precision
+c
+      eps = 1.d0
+   10 eps = eps/2.d0
+      tol1 = 1.d0 + eps
+      if (tol1 .gt. 1.d0) go to 10
+c
+c initialization
+c
+      a = ax
+      b = bx
+      fa = f(a)
+      fb = f(b)
+c
+c begin step
+c
+   20 c = a
+      fc = fa
+      d = b - a
+      e = d
+   30 if (dabs(fc) .ge. dabs(fb)) go to 40
+      a = b
+      b = c
+      c = a
+      fa = fb
+      fb = fc
+      fc = fa
+c
+c convergence test
+c
+   40 tol1 = 2.d0*eps*dabs(b) + 0.5*tol
+      xm = .5*(c - b)
+      if (dabs(xm) .le. tol1) go to 90
+      if (fb .eq. 0.d0) go to 90
+c
+c is bisection necessary
+c
+      if (dabs(e) .lt. tol1) go to 70
+      if (dabs(fa) .le. dabs(fb)) go to 70
+c
+c is quadratic interpolation possible
+c
+      if (a .ne. c) go to 50
+c
+c linear interpolation
+c
+      s = fb/fa
+      p = 2.d0*xm*s
+      q = 1.d0 - s
+      go to 60
+c
+c inverse quadratic interpolation
+c
+   50 q = fa/fc
+      r = fb/fc
+      s = fb/fa
+      p = s*(2.d0*xm*q*(q - r) - (b - a)*(r - 1.d0))
+      q = (q - 1.d0)*(r - 1.d0)*(s - 1.d0)
+c
+c adjust signs
+c
+   60 if (p .gt. 0.d0) q = -q
+      p = dabs(p)
+c
+c is interpolation acceptable
+c
+      if ((2.d0*p) .ge. (3.d0*xm*q - dabs(tol1*q))) go to 70
+      if (p .ge. dabs(0.5*e*q)) go to 70
+      e = d
+      d = p/q
+      go to 80
+c
+c bisection
+c
+   70 d = xm
+      e = d
+c
+c complete step
+c
+   80 a = b
+      fa = fb
+      if (dabs(d) .gt. tol1) b = b + d
+      if (dabs(d) .le. tol1) b = b + dsign(tol1, xm)
+      fb = f(b)
+      if ((fb*(fc/dabs(fc))) .gt. 0.d0) go to 20
+      go to 30
+c
+c done
+c
+   90 zeroin = b
+      return
+      end

# File codes/homework5/Makefile

+
+OBJECTS = functions.o quadrature.o test.o
+FFLAGS = -fopenmp
+
+.PHONY: test clean
+
+test: test.exe
+	./test.exe
+
+test.exe: $(OBJECTS) + gfortran$(FFLAGS) $(OBJECTS) -o test.exe + +%.o : %.f90 + gfortran$(FFLAGS) -c  $<  + +clean: + rm -f *.o *.exe *.mod + # File codes/homework5/README.txt + +Sample code to start with for Homework 5. + +For convenience a Makefile is provided: +$ make test
+
+to compile and run.
+
+A notebook describing Simpson's rule is in the notebook subdirectory, in
+various formats.
+

# File codes/homework5/functions.f90

+
+module functions
+
+    use omp_lib
+    implicit none
+    integer :: fevals(0:7)
+    real(kind=8) :: k
+    save
+
+contains
+
+    real(kind=8) function f(x)
+        implicit none
+        real(kind=8), intent(in) :: x
+        integer thread_num
+
+        ! keep track of number of function evaluations by
+        ! each thread:
+        thread_num = 0   ! serial mode