Scientific Computing Linear Algebra

1. Scientific ComputingLinear Algebra Dr. Guy Tel-Zur Sunset in Caruaru by Jaime JaimeJunior. publicdomainpictures.net Version 10-04-11, 14:00

2. MHJ Chapter 4 – Linear Algebra • In this talk we deal with basic matrix operations • Such as the solution of linear equations, calculate the inverse of a matrix, its determinant etc. • Here we focus in particular on so-called director elimination methods, which are in principle determined through a finite number of arithmetic operations. • Iterative methods will be discussed in connection with eigenvalue problems in MHJ chapter 12. • This chapter serves also the purpose of introducing important programming details such as handling memory allocation for matrices and the usage of the libraries which follow these lectures.

3. Libraries • LAPACK based on: • EISPACK – for solving symmetric, un-symmetric and generalized eigenvalue problems • LINPACK - linear equations and least square problems • BLAS: Basic Linear Algebra Subprogram • Levels I, II and III • Nelib: http://www.netlib.org solutions Building blocks

4. LAPACK - Linear Algebra PACKage LAPACK is written in Fortran90 and provides routines for solving systems of simultaneous linear equations, least-squares solutions of linear systems of equations, eigenvalue problems, and singular value problems. The associated matrix factorizations (LU, Cholesky, QR, SVD…) are also provided. Dense and banded matrices are handled, but not general sparse matrices. In all areas, similar functionality is provided for real and complex matrices, in both single and double precision.

5. EISPACK EISPACK is a collection of Fortran subroutines that compute the eigenvalues and eigenvectors of nine classes of matrices: complex general, complex Hermitian, real general, real symmetric, real symmetric banded, real symmetric tridiagonal, special real tridiagonal, generalized real, and generalized real symmetric matices. In addition, two routines are included that use singular value decomposition to solve certain least-squares problems.

6. LINPACK LINPACK is a collection of Fortran subroutines that analyze and solve linear equations and linear least-squares problems. The package solves linear systems whose matrices are general, banded, symmetric indefinite, symmetric positive definite, triangular, and tridiagonal square. In addition, the package computes the QR and singular value decompositions of rectangular matrices and applies them to least-squares problems. LINPACK uses column-oriented algorithms to increase efficiency by preserving locality of reference.

7. LINPACK and EISPACK are based on BLAS I • LAPACK is based on BLAS III

8. BLAS • The BLAS (Basic Linear Algebra Subprograms) are routines that provide standard building blocks for performing basic vector and matrix operations. • The Level 1 BLAS perform scalar, vector and vector-vector operations. • The Level 2 BLAS perform matrix-vector operations • The Level 3 BLAS perform matrix-matrix operations. • Because the BLAS are efficient, portable, and widely available, they are commonly used in the development of high quality linear algebra software, LAPACK for example.

9. LAPACK / CLAPACK / ScaLAPACK for Windows http://icl.cs.utk.edu/lapack-for-windows/

10. For all matrix/problem structures:Ex: LAPACK Table of Contents • BD – bidiagonal • GB – general banded • GE – general • GG – general , pair • GT – tridiagonal • HB – Hermitian banded • HE – Hermitian • HG – upper Hessenberg, pair • HP – Hermitian, packed • HS – upper Hessenberg • OR – (real) orthogonal • OP – (real) orthogonal, packed • PB – positive definite, banded • PO – positive definite • PP – positive definite, packed • PT – positive definite, tridiagonal • SB – symmetric, banded • SP – symmetric, packed • ST – symmetric, tridiagonal • SY – symmetric • TB – triangular, banded • TG – triangular, pair • TB – triangular, banded • TP – triangular, packed • TR – triangular • TZ – trapezoidal • UN – unitary • UP – unitary packed Source: CS267 Lecture 13

11. An Example - DGESV • Ax=b • DGESV computes the solution to a real system of linear equations A * X = B, where A is an N-by-N matrix and X and B are N-by-NRHS matrices. • The LU decomposition with partial pivoting and row interchanges is used to factor A asA = P * L * U,where P is a permutation matrix, L is unit lower triangular, and U is upper triangular. The factored form of A is then used to solve the system of equations A * X = B.

12. Arguments N (input) INTEGER , The number of linear equations, i.e., the order of the matrix A. N >= 0. NRHS(input) INTEGER The number of right hand sides, i.e., the number of columns of the matrix B. NRHS >= 0. A(input/output) DOUBLE PRECISION array, dimension (LDA,N) On entry, the N-by-N coefficient matrix A. On exit, the factors L and U from the factorization A = P*L*U; the unit diagonal elements of L are not stored. LDA(input) INTEGER The leading dimension of the array A. LDA >= max(1,N). IPIV (output) INTEGER array, dimension (N) The pivot indices that define the permutation matrix P; row i of the matrix was interchanged with row IPIV(i). B(input/output) DOUBLE PRECISION array, dimension (LDB,NRHS) On entry, the N-by-NRHS matrix of right hand side matrix B. On exit, if INFO = 0, the N-by-NRHS solution matrix X. LDB(input) INTEGER The leading dimension of the array B. LDB >= max(1,N). INFO (output) INTEGER = 0: successful exit< 0: if INFO = -i, the i-th argument had an illegal value> 0: if INFO = i, U(i,i) is exactly zero. The factorization has been completed, but the factor U is exactly singular, so the solution could not be computed.

13. Working with a Linear Algebra package using Visual Studio Demo Work dir: C:\Users\telzur\Documents\Weizmann\ScientificComputing\SC2011B\Lectures\04\LAPACK\CLAPACK-EXAMPLE\Release

14. Compare with the following MATLAB code Demo clear all close all disp('solves Ax=b') A= [76, 27, 18; 25, 89, 60; 11, 51, 32] b=[10, 7, 43] invA=inv(A) x=b*inv(A) [L,U]=lu(A) inv(U)*inv(L)*b' Code location: C:\Users\telzur\Documents\BGU\Teaching\ComputationalPhysics\2011A\Lectures\04\lpack.m

15. DGESV SUBROUTINE DGESV( N, NRHS, A, LDA, IPIV, B, LDB, INFO ) * * -- LAPACK driver routine (version 3.2) -- * -- LAPACK is a software package provided by Univ. of Tennessee, -- * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- * November 2006 * * .. Scalar Arguments .. INTEGER INFO, LDA, LDB, N, NRHS * .. * .. Array Arguments .. INTEGER IPIV( * ) DOUBLE PRECISION A( LDA, * ), B( LDB, * ) * .. * * Purpose * ======= * * DGESV computes the solution to a real system of linear equations * A * X = B, * where A is an N-by-N matrix and X and B are N-by-NRHS matrices. * * The LU decomposition with partial pivoting and row interchanges is * used to factor A as * A = P * L * U, * where P is a permutation matrix, L is unit lower triangular, and U is * upper triangular. The factored form of A is then used to solve the * system of equations A * X = B. *

16. * Arguments * ========= * * N (input) INTEGER * The number of linear equations, i.e., the order of the * matrix A. N >= 0. * * NRHS (input) INTEGER * The number of right hand sides, i.e., the number of columns * of the matrix B. NRHS >= 0. * * A (input/output) DOUBLE PRECISION array, dimension (LDA,N) * On entry, the N-by-N coefficient matrix A. * On exit, the factors L and U from the factorization * A = P*L*U; the unit diagonal elements of L are not stored. * * LDA (input) INTEGER * The leading dimension of the array A. LDA >= max(1,N). * * IPIV (output) INTEGER array, dimension (N) * The pivot indices that define the permutation matrix P; * row i of the matrix was interchanged with row IPIV(i). * * B (input/output) DOUBLE PRECISION array, dimension (LDB,NRHS) * On entry, the N-by-NRHS matrix of right hand side matrix B. * On exit, if INFO = 0, the N-by-NRHS solution matrix X. * * LDB (input) INTEGER * The leading dimension of the array B. LDB >= max(1,N). * * INFO (output) INTEGER * = 0: successful exit * < 0: if INFO = -i, the i-th argument had an illegal value * > 0: if INFO = i, U(i,i) is exactly zero. The factorization * has been completed, but the factor U is exactly * singular, so the solution could not be computed. * * ===================================================================== * Cont’

17. * .. External Subroutines .. EXTERNAL DGETRF, DGETRS, XERBLA * .. * .. Intrinsic Functions .. INTRINSIC MAX * .. * .. Executable Statements .. * * Test the input parameters. * INFO = 0 IF( N.LT.0 ) THEN INFO = -1 ELSE IF( NRHS.LT.0 ) THEN INFO = -2 ELSE IF( LDA.LT.MAX( 1, N ) ) THEN INFO = -4 ELSE IF( LDB.LT.MAX( 1, N ) ) THEN INFO = -7 END IF IF( INFO.NE.0 ) THEN CALL XERBLA( 'DGESV ', -INFO ) RETURN END IF * * Compute the LU factorization of A. * CALL DGETRF( N, N, A, LDA, IPIV, INFO ) IF( INFO.EQ.0 ) THEN * * Solve the system A*X = B, overwriting B with X. * CALL DGETRS( 'No transpose', N, NRHS, A, LDA, IPIV, B, LDB, $ INFO ) END IF RETURN * * End of DGESV * END Cont’

18. LAPACK Benchmark MFLOPS is calculated as: MFLOPS = (1x10**-6) * [(2/3)*N**3 + 2*N**2] / (CPU seconds) http://www.ats.ucla.edu/clusters/common/software/libraries/lapack_benchmark.htm

19. Cost • Solving A*x=b using GE • Factorize A = L*U using GE (cost = 2/3 n3 flops) • Solve L*y = b for y, using substitution (cost = n2 flops) • Solve U*x = y for x, using substitution (cost = n2 flops) PDGESV = ScaLAPACK Parallel LU (“Linpack Benchmark”),

20. Mathematicalintermezzo

21. http://ocw.mit.edu/courses/mathematics/18-085-computational-science-and-engineering-i-fall-2008/http://ocw.mit.edu/courses/mathematics/18-085-computational-science-and-engineering-i-fall-2008/

22. The Kn Matrices What are the properties of K ? The next slides are from: Computational Science and Engineering, by Gilbert Strang. An excellent recommended book. His course is available online from MIT Opencourseware. Very Recommended!!!!

23. Kn Properties • These matrices are symmetric. • The matrices Kn are sparse. • These matrices are tridiagonal. • The matrices have constant diagonals. • All the matrices K = Kn are invertible. • The symmetric matrices Kn are positive definite. Source: Computational Science and Engineering, by Gilbert Strang

24. In Signal Processing D=Kn/4 is a “High-Pass” Filter. Du picks out the rapidly varying parts of a vector u • K are called Toeplitz Matrix and MATLAB has a function for generating such matrices, e.g. K = toeplitz([2 -1 zeros(l,2)])constructs K4 from row 1 Source: Computational Science and Engineering, by Gilbert Strang

25. More properties • (Pivots) An invertible matrix has n nonzero pivots. • A positive definite symmetric matrix has n positive pivots. • (Eigenvalues) An invertible matrix has n nonzero eigenvalues. • A positive definite symmetric matrix has n positive eigenvalues. (Positive definite means a quadratic form that satisfies: XTAX>0 ) Source: Computational Science and Engineering, by Gilbert Strang

26. More special matrices T=Top, B=Bottom Source: Computational Science and Engineering, by Gilbert Strang

27. Building K,T,B,C in Matlab Demo Make a demo also using Octave Source: Computational Science and Engineering, by Gilbert Strang

29. Demo Matlab Demos K4=toeplitz([2 -1 0 0]) C4=toeplitz([2 -1 0 -1]) inv(K4) …OK inv(C4) …not OK  singular eig(K4) positive >0  positive definite eig(C4) >=0  semi positive definite Both are symmetric: try transpose(K4) Check: [L,U]=lu(T4) all pivots are 1 [L,U]=lu(B4) 0 on the diagonal of U  B isn’t invertibe inv(T4) …OK inv(B4) …not OK  singular

30. Demo demos e=ones(8,1) B.e=0 . = dot product BT=transpose(B)  B is symmetric System to solve: Bu=f uTBT=fT BT.e=B.e=0 Thefore fT.e=0

31. Demo Continued det(K8) = Π(diagonal of U) = 2/1 * 3/2 * 4/3 … * 9/8 = 9

32. Demo For large size n THIS IS THE CODE TO CREATE K,T,B,C AS SPARSE MATRICES Then Matlab will not operate on all the zeros! function K=SKTBC(type,n,sparse)% SKTBC Create finite difference model matrix.% K=SKTBC(TYPE,N,SPARSE) creates model matrix TYPE of size N-by-N.% TYPE is one of the characters 'K', 'T', 'B', or 'C'.% The command K = SKTBC('K', 100, 1) gives a sparse representation% K=SKTBC uses the defaults TYPE='K', N=10, and SPARSE=false.% Change the 3rd argument from 1 to 0 for dense representation!% If no 3rd argument is given, the default is dense% If no argument at all, KTBC will give 10 by 10 matrix Kif nargin<1, type='K'; endif nargin<2, n=10; ende=ones(n,1);K=spdiags([-e,2*e,-e],-1:1,n,n);switch type  case 'K'  case 'T'   K(1,1)=1;  case 'B'   K(1,1)=1;   K(n,n)=1;  case 'C'   K(1,n)=-1;   K(n,1)=-1;  otherwise    error('Unknown matrix type.');endif nargin<3 | ~sparse  K=full(K);end http://math.mit.edu/classes/18.085/create-sparse.html

34. Source: Numerical recipes in FORTRAN77: the art of scientific computing, page 65.  By William H. Press

35. In Mathematics: Vectors and Matrices Are mapped to Computers as Memory arrays Fixed Memory allocation vs. Dynamic Memory Allocation Compile time vs. Run time In the next slides we will study two programs that address these issues Let’s start with Table 4.2 in the book: Matrix handling program where arrays are defined at compilation time – Next slide Use Visual Studio for the demo solution file under: chap4_static

36. int main() { intk,m, row = 3, col = 5; intvec[5]; // line a: a standard C++ declaration of a vector intmatr[3][5]; // line b: a standard fixed-size C++ declaration of a matrix for(k = 0; k < col; k++) vec[k] = k; // data into vector[] for(m = 0; < row; m++) { // data into matr[][] for(k = 0; k < col ; k++) matr[m][k] = m + 10 * k; } printf("\n Vector data in main():\n”); // print vector data for(k = 0; k < col; k++) printf("vetor[%d]= %d ",k, vec[k]); printf("\n Matrix data in main():"); for(m = 0; m < row; m++) { printf(“\n”); for(k = 0; k < col; k++) printf("m atr[%d][%d]= %d ",m,k,matr[m][k]); } printf(“\n”); sub_1(row, col, vec, matr); // line c: transfer vec[] and matr[][] addresses to func. sub_1(). return 0; } // End: function main() void sub_1(int row, int col, intvec[], intmatr[][5]) { //line d: a func. def. intk,m; printf("\n Vetor data in sub1():\n"); // print vector data for(k = 0; k < col; k++) printf("vetor[%d]= %d ",k, vec[k]); printf("\n Matrix data in sub1():"); for(m = 0; m < row; m++) { printf(“\n”); for(k = 0; k < col; k++) { printf("matr[%d][%d]= %d ",m, k, matr[m][k]); } } printf(“\n”); } // End: function sub_1() Table 4.2 Demo equiv to int *vec equiv to int (*matr)[5]

37. Using Visual Studio 2010 Make demo in class

38. Static program execution

39. Table 4.3: Matrix handling program with dynamic array allocation. התכנית ב 4.3 לא מתקמפלת ומלאה באגים. בעתיד היא תתווסף כתכנית לדוגמה. בשלב הזה נתייחס אליה כאל פסאודו-קוד מבלי להריצה

40. #include <stdio.h> int main() { int *vec; // line a int **matr; // line b int m, k, row, col, total = 0; printf("\n Read in number of rows= "); // line c scanf("%d",&row); printf("\n Read in number of column= "); scanf("%d", &col); vec = new int [col]; // line d matr = (int **)matrix(row, col, sizeof(int)); // line e for(k = 0; k < col; k++) vec[k] = k; // store data in vector[] for(m = 0; < row; m++) { // store data in array[][] for(k = 0; k < col; k++) matr[m][k] = m + 10 * k; } printf("\n Vetor data in main():\n"); // print vector data for(k = 0; k < col; k++) printf("vetor[%d]= %d ",k,vec[k]); printf("\n Array data in main():"); for(m = 0; m < row; m++) { printf("\n"); for(k = 0; k < col; k++) { printf("m atrix[%d][%d]= %d ",m, k, matr[m][k]); } } printf("\n"); for(m = 0; m < row; m++) { // access the array for(k = 0; k < col; k++) total += matr[m][k]; } printf("\n Total= %d\n",total); sub_1(row, col, vec, matr); free_matrix((void **)matr); // line f delete [] vec; // line g return 0; } // End: function main() we declare a pointer to an integer declares a pointer-to-a-pointer which will contain the address to a pointer of row vectors, each with col integers read in the size of vec[] and matr[][] through the numbers row and col we reserve memory for the vector we use a user-defined function to reserve necessary memory for matrix[row][col] and again matr contains the address to the reserved memory location. The remaining part of the function main() are as in the previous case down to line f. the same procedure is performed for vec[]

41. Continued from previous slide void sub_1(int row, int col, intvec[], int ??matr) // line h { intk,m; printf("\n Vetor data in sub1():\n"); // print vector data for(k = 0; k < col; k++) printf("vetor[%d]= %d ",k, vec[k]); printf("\n Matrix data in sub1():"); for(m = 0; m < row; m++) { printf("\n"); for(k = 0; k < col; k++) { printf("matrix[%d][%d]= %d ",m,k,matr[m][k]); } } printf("\n"); } // End: function sub_1() in line h an important difference from the previous case occurs. First, the vector declaration is the same, but the matr declaration is quite different. The corresponding parameter in the call to sub_1[] in line g is a double pointer. Consequently, matr in line h must be a double pointer

42. /* * The function * void **matrix( ) * reserves dynamic memory for a two-dimensional matrix * using the C++ command new . No initialization of the elements . * Input data : * int row - number of rows * int col - number of columns * intnum_bytes - number of bytes for each element * Returns avoid **pointer to the reserved memory location. */ void **matrix( int row , int col , intnum_bytes ) { int i , num; char **pointer, *ptr; pointer = new(nothrow) char* [ row ] ; if ( !pointer ) { cout << "Exeption handling Memory aloation failed"; cout << "for"<< row << "row addresses!"<< endl; return NULL; } i = ( row * col * num_bytes ) / size of ( char ) ; pointer [ 0 ] = new( nothrow ) char [ i ] ; if ( !pointer [ 0 ] ) { cout << "Exeption handling :Memory allocation failed"; cout << "for address to " << i << " characters!"<< endl ; return NULL; } ptr = pointer [ 0 ] ; num = col * num_bytes ; for ( i = 0 ; i < row ; i ++ , ptr += num ) { pointer [ i ] = ptr ; } return ( void **) pointer ; } // end : function void **matrix ( ) lib.cpp

43. Skip to section 4.4

44. C++ and Fortran features of matrix handling If we store a matrix in the sequence a11 a12 . . . a1n a21 a22 . . . a2n . . . ann This is called “row-major” order (we go along a given row i and pick up all column elements j) C++ stores them by row-major. We cloud also store in column-major order a11 a21 . . . an1 a12 a22 . . . an2 . . . ann. Fortran stores matrices by column-major,

45. 4.4 Linear Systems • Gauss Elimination • Upper/Lower triangular matrix • Solve Linear System • LU algorithm • Cholesky decomposition (a special case) • QR

46. (4.4) Linear Systems: An Example

47. This is the stiffness matrix, K, we met earlier! …And we switched from a Differential Equation to Linear Algebra.

48. LU Decomposition (Factorization) Theorem: if the n-by-n matrix A is non singular, there exists a permutation matrix P (the identity matrix with its rows permuted), a non-singular lower triangular matrix L and a non-singular upper triangular matrix U such that: A=PLU. To solve Ax=b, we solve the equivalent PLUx=b as follows: LUx=P-1b=PTb (permute entries of b) Ux=L-1 (PTb) (forward substitution) x=U-1 (L-1PTb) (backward substitution) For further reading: Applied Numerical Linear Algebra by James Demmel (SIAM)

49. LU Decomposition (Factorization) Reference: MHJ chap 4 (section 4.4.2) det{L}=1

50. Cholesky’s Factorization MHJ Chapter 4, page 87 (2009 edition) One has to check first if A is positive definite!