1 / 53

Computational Physics Linear Algebra

Computational Physics Linear Algebra. Dr. Guy Tel- Zur. Sunset in Caruaru by Jaime JaimeJunior . publicdomainpictures.net. Version 4-11-10, 14:00. MHJ Chapter 4 – Linear Algebra. In this talk we deal with basic matrix operations Such as the solution of linear equations, calculate the

kato
Download Presentation

Computational Physics Linear Algebra

An Image/Link below is provided (as is) to download presentation Download Policy: Content on the Website is provided to you AS IS for your information and personal use and may not be sold / licensed / shared on other websites without getting consent from its author. Content is provided to you AS IS for your information and personal use only. Download presentation by click this link. While downloading, if for some reason you are not able to download a presentation, the publisher may have deleted the file from their server. During download, if you can't get a presentation, the file might be deleted by the publisher.

E N D

Presentation Transcript

  1. Computational PhysicsLinear Algebra Dr. Guy Tel-Zur Sunset in Caruaru by Jaime JaimeJunior. publicdomainpictures.net Version 4-11-10, 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

  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. Working with a Linear Algebra package using Visual Studio Demo Work dir: C:\Users\telzur\Documents\BGU\Teaching\ComputationalPhysics\2011A\Lectures\04\LAPACK\CLAPACK-EXAMPLE\Release>

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

  12. 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. *

  13. * 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’

  14. * .. 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’

  15. Mathematicalintermezzo

  16. 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/

  17. 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!!!!

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

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

  20. 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. Source: Computational Science and Engineering, by Gilbert Strang

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

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

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

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

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

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

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

  28. 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 demonstrate 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

  29. 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]

  30. Using Visual Studio 2010 Make demo in class

  31. Static program execution

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

  33. #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[]

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

  35. /* * 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

  36. Skip to section 4.4

  37. 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,

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

  39. (4.4) Linear Systems: An Example

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

  41. Case #1: Assume first that the function f does not depend on u(x). Then our linear equation reduces to Au = f , (4.8) which is nothing but a simple linear equation with a tridiagonal matrix A. We will solve such a system of equations in subsection 4.4.3.

  42. Case #2 If we assume that our boundary value problem is that of a quantum mechanical particle confined by a harmonic oscillator potential, then our function f takes the form (assuming that all constants m = h_bar = omega = 1 with λ being the eigenvalue. Inserting this into our equation, we define first a new matrix A as להבנה נוספת ראש השקף הבא

  43. We will solve this type of equations in chapter 12.

  44. (4.4.2) LU Decomposition LU Factorization Continue from page 83 at MHJ chap 4 PDF

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

More Related