Reminder. The second and last midterm exam will be held as scheduled next Tuesday, Oct. 30. The format will be similar to the previous midterm. Coverage: Fortran programming.
There are several versions of LAPACK; the newest (Version 3) is fully compatible with Fortran 95 and, arguably, is the most convenient to use. For today's exercise, however, we'll use the Fortran-77 calling interface to LAPACK, because it's still the most widely used. Both LAPACK and BLAS are available on the machines in ECA 221.
Exercise 1. Write a subroutine,
linsolve(a,n,b,info)that solves the system Ax=b for a given n×n matrix A and n-vector b. Here a and b are intent(inout) variables that, on input, define the n×n matrix A and n-vector b, respectively. The input argument n gives the size of the problem. On output, b should hold the solution x; your routine is free to overwrite the contents of a. The output argument info is an integer flag that should be set to zero unless an error is encountered in solving the system; in this case, set info to any nonzero value of your choice.
Put another way, linsolve overwrites b with the quantity A-1b, assuming that A is invertible. It is much more efficient to use a routine like linsolve than to form A-1 explicitly, particularly when the system is large.
The declarations will look like
subroutine linsolve(a,n,b,info)You may add linsolve to the module linsys that you wrote last week. Change D as necessary if you altered the name of the double-precision kind variable in precision. (A use precision statement is not necessary, since one appears at the top of linsys. However, a use precision statement would be necessary inside linsolve were a standalone subroutine outside of a module.)
integer,intent(in):: n ! the size of the problem
real(D),intent(inout):: a(n,n) ! overwritten by DGESV
real(D),intent(inout):: b(n) ! holds RHS in input and solution on output
integer,intent(out):: info ! return code
Your routine will call the LAPACK library routine dgesv, with the Fortran 77 calling sequence (i.e., the first one listed on the manual page). Your routine must provide the necessary workspace for a problem of arbitrary size. Your routine must not print anything; all communication is through its arguments (similarly to last week's trisolve).
Test your code with your choice of inputs.
Exercise 2. Write a subroutine,
multinewt(f,df,n,x,epsilon,info)that implements Newton's method for an arbitrary differentiable vector-valued function f. Newton's method for vector-valued functions looks just like the scalar case above, except that the division by f ' is replaced by Df(xk)-1 f(xk), where Df is the Jacobian matrix of partial derivatives of f, evaluated at xk.
In the Fortran code, f and df are subroutines that define the function and its Jacobian. The input argument n gives the size of the problem (i.e., the number of component functions in f). The argument x on input is a starting guess for the solution; on return, x holds an estimate of the solution when the routine is successful. The output argument info is an integer flag that is set to zero if the routine completes successfully and to any nonzero value of your choice otherwise. The input argument epsilon defines the stopping criterion. Your routine should halt after 100 iterations or when the difference between two successive iterates is less than epsilon, whichever occurs first. You may define the difference between two vectors x and y as
maxk |xk - yk|,that is, the maximum absolute value of the difference between any two corresponding components.
The first two arguments, which require an interface block, define subroutines that evaluate f and Df, respectively. They are invoked by
call f(x,n,fcn)respectively. In each case, x is an input argument that specifies the current point. The output argument fcn is an n-vector whose contents are f(x). The output argument dfcn is an n×n matrix whose contents are Df(x).
call df(x,n,dfcn)
Use linsolve to evaluate Df(xk)-1 f(xk). Do not form Df(xk)-1 explicitly!
Here is a code skeleton to get you started.
For testing purposes, try solving the system
x2y - xy2 = -6
x3y + xy3 = -10
Write a module akin to last week's myfuncs.f90 that defines this function and its Jacobian matrix. Read in the starting point and tolerance from the standard input. There are two solutions. Which one you find depends on the starting guess that you supply.
You will need to link with external libraries to obtain the object code for dgesv and the BLAS. Modify the makefile that you wrote last week so that it includes the line
LIBS = -llapack-3 -lblas-3This tells the compiler drive to look for the LAPACK library and the BLAS library. The -l options must be listed in the order shown, because LAPACK depends on BLAS, and the linker processes libraries in the order in which they are given on the command line (from left to right). No spaces may appear between -l and the library names.
Edit the remainder of your previous makefile appropriately and use it to compile and link your test program.
Note to Macintosh users: Apple includes a full implementation of LAPACK and the BLAS on all OS X machines. To link with it, replace the above line with
LIBS = -framework Accelerate
Exercise 3. Use the same multinewt module that you wrote in Exercise 2 to solve the system
s(y - x) = 0There are three solutions, each of which defines a fixed point for the Lorenz equations, which are an archetype for the chaotic behavior of the atmosphere. Use the parameters s=10, b=8/3, and r=28. I suggest that you write a main program that reads in one 3-vector as a guess, then run your main program with multiple guesses until you locate the three solutions. (You can solve this system by hand to double-check your computed answers.)
x(r - z) - y = 0
xy - bz = 0.
Copyright(c) 2007 by Eric J. Kostelich. All rights reserved.