Numerical Methods, MAT 350
Autumn 2001

Linear Equations and Eigenvalue problems

Due Tuesday November 20 at 5:00 PM

This homework will emphasize the numerical solutions of linear systems and eigenvalue problems. This exercise is not intended to be conceptually difficult, rather it is intended to give you practical experience in using standard subroutines. Thus, if you get stuck with a programming problem, don't waste too much time: find someone to help you.

For this assignment you will use routines from the BLAS and LAPACK libraries. We will use the DOUBLE PRECISION version of the routines. In this case, the subroutines will all start with the letter D. Use GAMS to figure out which subroutines you need and how to call them.

Using Computers

It is recommended that you use the UNIX computers to complete this assignment. The BLAS and LAPACK libraries are already installed on these computers. If you use Microsoft Windows, you will have to download and properly install the needed libraries. Experience has shown that this produces severe headache and eventual despair.

If you are using FORTRAN, you need "link" your program to the libraries when you compile:

g77 myprogram.f -llapack -lblas

If you are using C/C++, you can still link to these libraries. You can construct prototypes based on CLAPACK library routines in GAMS. In that case, you will need the include statement #include <g2c.h> at the beginning of your program. Compile and link your program using:

gcc myprogram.c -llapack -lblas -lg2c
If you are using C++, you can still link to these C routines. In the file with your C++ program you need to declare the prototypes of any C routines in the following manner:
extern "C" {
/* Subroutine */ int dgetrf_(long int *m, long int *n, double *a, 
   long int *lda, long int *ipiv, long int *info);

// You can add various other C prototypes here

More information on mixing C and C++ can be found at the C++ FAQ.

Working with matrices

Working with matrices is a little tricky, so an example may be helpful. Consider the following matrix:
  10  9  8  7
   2  3  4  5
   6  2 12 14
   3 11 13  1
Here is some code from a C program which finds the LU factorization of this matrix:
#define NMAX 10  /* NMAX is the maximum size of the matrices */

/* In C, we store the matrix in one long array.
   Note that we go down the columns of the matrix.
   Also, the physical size (lda=4) is equal to the logical 
   size (n=4) since the columns of the matrix are right next 
   to each other in memory.  */

  long int n=4, lda=4;
  double a[NMAX*NMAX]={10,2,6,3, 9,3,2,11, 8,4,12,13, 7,5,14,1},

/* Call the LAPACK routine which performs the LU factorization. */


/*  Print out the LU factorization and the pivots. */

  printf("The pivots are:\n");
  for(i=0; i<n; i++)
    printf(" %li\n",ipiv[i]);

  printf("\nThe factorization is:\n");
  for(i=0; i<n; i++){
    for(j=0; j<n; j++)
      /* Note that lda is the physical size of the matrix */
      printf(" %lg",a[i+j*lda]);
In Fortran, there are several ways you can define the matrix. For instance, you can use a data statement as shown in the following code taken from a Fortran program:
      PARAMETER (lda=4,n=4)
      INTEGER ipiv(lda)
      DOUBLE PRECISION a(lda,lda)

C Here, it is essential that the physical size (lda=4)
C is equal to the logical size (n=4) since the numbers
C are being read in as an array.
C Note that we go down the columns of the matrix.

      DATA a/10,2,6,3,  9,3,2,11, 
     *       8,4,12,13, 7,5,14,1/

C Call the LAPACK routine which performs the LU factorization. */

      CALL DGETRF( n, n, a, nmax, ipiv, info )

The assignment

In this assignment, you will find yourself modifying code you have written previously. This is a good skill to develop. Sometimes it is useful to "comment out" the old stuff so that you (and I) can see what was done previously.

For this assignment, you are required to hand in a floppy disk or send me an E-mail containing the programs you wrote. However, program outputs, hand calculations, et cetera can be on paper.

  1. Write a short program wich prints out the 4 by 4 matrix that you wrote down in your last assignment.
  2. Write a program to that computes the LU factorization of this matrix
  3. Find another subroutine which solves a real general system of linear equations A x = b using the LU factorization computed by DGETRF.
  4. Include this routine in your program and use it to solve the set of equations you wrote down in your last assignment. Does your answer agree with your homework?
  5. Consider the following 4 by 4 symmetric matrix:
      1  -2   3  -4
     -2   5  -8  11
      3  -8  14 -20
     -4  11 -20  30
    Find its Cholesky factorization by hand.
  6. Find a subroutine which computes the Cholesky factorization of a real symmetric positive definite matrix.
  7. Write a program which uses this subroutine. You should:
  8. Find a subroutine which computes all eigenvalues and, optionally, eigenvectors of a real symmetric matrix.
  9. Write a program which uses this subroutine to calculate the eigenvalues (only) of a 6 by 6 Hilbert matrix.
  10. In the previous homework and in class, you learned about some properties of determinants.
  11. As your final task, look at the eigenvalues of a 25 by 25 Hilbert matrix. What happens to the smallest eigenvalues?