Autumn 2001

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.

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 -lg2cIf 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.

10 9 8 7 2 3 4 5 6 2 12 14 3 11 13 1Here 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 )

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.

- Write a short program wich prints out the 4 by 4 matrix that you wrote down in your last assignment.
- Write a program to that computes the LU factorization of this matrix
- Find an LAPACK subroutine which computes an LU factorization of a real general matrix, using partial pivoting with row interchanges.
- Print out the factorized matrix along with the array containing the pivots.
- Does the factorization look anything like the factorization you had in your homework? (It might not, because of the pivoting.)

- Find another subroutine which solves a real general system of
linear equations
**A**=**b**using the LU factorization computed by DGETRF.- Hint: the subroutine DGESV is
*not*what you are looking for. - In our case, we only want to solve one set of equations,
so
**b**and**x**will only have one column.

- Hint: the subroutine DGESV is
- 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?
- 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. - Find a subroutine which computes the Cholesky factorization of a real symmetric positive definite matrix.
- Write a program which uses this subroutine. You should:
- Check your program by using it on the example you did by hand above. (Print out your results.)
- Modify your program (this would be a good place for you
to "comment out" your old code), so that
it computes the Cholesky factorization of an
*n*by*n*Hilbert matrix. - Find the Cholesky factorization for a 6 by 6 Hilbert matrix. Print out your results.

- Find a subroutine which computes all eigenvalues and, optionally, eigenvectors of a real symmetric matrix.
- Write a program which uses this subroutine to calculate the eigenvalues (only) of a 6 by 6 Hilbert matrix.
- In the previous homework
and in class, you learned
about some properties of determinants.
- Use the eigenvalues you calculated above to find the determinant of a 6 by 6 Hilbert matrix.
- Now, use the results of the Cholesky factorization to find the determinant of a 6 by 6 Hilbert matrix.
- Do your two answers agree?

- As your final task, look at the eigenvalues of a 25 by 25 Hilbert matrix. What happens to the smallest eigenvalues?