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},
         b[NMAX]={23,16,17,21};

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

  dgetrf_(&n,&n,a,&lda,ipiv,&info);

/*  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]);
    printf("\n");
  }

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.