MacKinnon - Computational Physics (523159), страница 7
Текст из файла (страница 7)
Unfortunately this varies from language to language: in C(++) and Pascal the first index variesslowest and the matrix is stored one complete row after another; whereas in FORTRAN the first index variesfastest and the matrix is stored one complete column after another. It is generally most efficient to accessthe matrix elements in the order in which they are stored. Hence the simple matrix algebra, A B C,should be written in C1 as 1 In C(++) it is usual to use double instead of float when doing numerical programming. In most cases single precision is notaccurate enough.Matrix Algebra27const int N =??;int i, j;double A[N][N], B[N][N], C[N][N];for( i = 0; i < N; i++)for( j = 0; j < N; j++)A[i][j] = B[i][j] + C[i][j];or its equivalent using pointers.
In FORTRAN90 this should be writteninteger :: i,jinteger, parameter :: N = ??real, double precision :: A(N,N), B(N,N), C(N,N)do i = 1,Ndo j=1,NA(j,i) = B(j,i) + C(j,i)end doend doNote the different ordering of the loops in these 2 examples. This is intended to optimise the order in whichthe matrix elements are accessed.
It is perhaps worth noting at this stage that in both C++ and FORTRAN90it is possible to define matrix type variables (classes) so that the programs could be reduced tomatrix A(N,N), B(N,N), C(N,N);A = B + C;The time taken to add or subtract 2# matrices is generally proportional to the total number ofmatrix elements, , although this may be modified by parallel architecture.¥3.3.2 Multiplication of MatricesUnlike addition and subtraction of matrices it is difficult to give a general machine independent rule forthe optimum algorithm for the operation. In fact matrix multiplication is a classic example of an operationwhich is very dependent on the details of the architecture of the machine.
We quote here a general purposeexample but it should be noted that this does not necessarily represent the optimum ordering of the loops.const int L = ??, M = ??, N = ??;int i, j, k;double A[L][N], B[L][M], C[M][N], sum;for ( j = 0; j < N; i++)for ( i = 0; i < L; j++){for ( sum = 0, k = 0; k < M; k++)sum += B[i][k] * C[k][j];A[i][j] = sum;}in C.The time taken to multiply 2 modified by parallel processing.¥$matrices is generally proportional to , although this may be3.4 Elliptic Equations — Poisson’s EquationAt this point we digress to discuss Elliptic Equations such as Poisson’s equation.
In general these equationsare subject to boundary conditions at the outer boundary of the range; there are no initial conditions, suchas we would expect for the Wave or Diffusion equations. Hence they cannot be solved by adapting themethods for simple differential equations.Matrix Algebra283.4.1 One Dimension0Y½ jk"]KÑ&\[K 0We start by considering the one dimensional Poisson’s equation(3.2)½ QeT $ f ½ Q ½ QS#T O K 0 Yj Q(3.3)j Q jk"]KÅhK Q 8gO KÑ& .where we defineõ½ "<KÑ&U½ at KÅK and ½I"]KÑ&;½ S#T at KõK S#T ,The boundary conditions are usually of the form½ and ½ S#T are both known the 8 althoughsometimes the condition is on the first derivative. Since equations (3.3) may be written asand 8$ f ½£T ½ K 0 jTU$}½(3.4)½ eTU$ f ½ 0 OO K 0 j $}½ S#T[(3.5)The 2nd derivative may be discretised in the usual way to giveThis may seem trivial but it maintains the convention that all the terms on the left contain unknowns andeverything on the right is known.
It also allows us to rewrite the (3.3) in matrix form as%&(*) %&(*)%&(*)&) &)&)&) &)&)&) &)&)&) &)&)&) &)&)&) &)&)&)&)&)......... ..&) & ..)& ..)&) &)&)(3.6)&) &)&)&) &)&)........ .......'+ ' +'+$f $ f $ f½£T½0½½ÑQ $f $ f $ fwhich is a simple matrix equation of the formAx½ eT½O K1K10ñ0ñjj 0TU$G½1O K10ñj OO K 0 jQO K1K10ñ0ñjj e$}T ½ S#TOb(3.7)in which A is tridiagonal. Such equations can be solved by methods which we shall consider in section 3.5.For the moment it should suffice to note that the tridiagonal form can be solved particularly efficiently andthat functions for this purpose can be found in most libraries of numerical functions.There are several points which are worthy of note.We could only write the equation in this matrix form because the boundary conditions allowed us toeliminate a term from the 1st and last lines.Periodic boundary conditions, such ascan be implemented, but they have theeffect of adding a non–zero element to the top right and bottom left of the matrix, and , sothat the tridiagonal form is lost.½"<K &o ½I"]KÑ&TTIt is sometimes more efficient to solve Poisson’s or Laplace’s equation using Fast Fourier Transformation (FFT).
Again there are efficient library routines available (Numerical Algorithms Group n.d.).This is especially true in machines with vector processors.3.4.2 2 or more DimensionsPoisson’s and Laplace’s equations can be solved in 2 or more dimensions by simple generalisations of theschemes discussed for 1D. However the resulting matrix will not in general be tridiagonal. The discretisedform of the equation takes the form½ ù á QegT ½ ù á QS#T ½ ù ge T\á Q ½ ù S#T^á Q,$ × ½ ù á Q, O K 0 j ù á QÑ[(3.8)Matrix Algebra29/ R: 8" : & : " : f & : [[ [" /: & : " f : & : " f : f & : [[ [ " f :/ & : [[ [" : & : " : f & : [ [ [" : &The 2 dimensional index pairs ,.- may be mapped on to one dimension for the purpose of setting upthe matrix.
A common choice is so–called dictionary order,Alternatively Fourier transformation can be used either for all dimensions or to reduce the problem totridiagonal form.3.5 Systems of Equations and Matrix InversionWe now move on to look for solutions of problems of the form¥where A is ancase where B is an¥matrixand X and B areunit matrix.AX¥B(3.9)matrices. Matrix inversion is simply the special3.5.1 Exact MethodsMost library routines, for example those in the NAG (Numerical Algorithms Group n.d.) or LaPack (LapackNumerical Library n.d.) libraries are based on some variation of Gaussian elimination. The standardprocedure is to call a first function which performs an LU decomposition of the matrix A,ALU(3.10)where L and U are lower and upper triangular matrices respectively, followed by a function which performsthe operationsYX L e e T T BU(3.11)Y(3.12)on each column of B.The procedure is sometimes described as Gaussian Elimination.
A common variation on this procedure is partial pivoting. This is a trick to avoid numerical instability in the Gaussian Elimination (or LUdecomposition) by sorting the rows of A to avoid dividing by small numbers.There are several aspects of this procedure which should be noted:The LU decomposition is usually done in place, so that the matrix A is overwritten by the matricesL and U.The matrix B is often overwritten by the solution X.A well written LU decomposition routine should be able to spot when A is singular and return a flagto tell you so.Often the LU decomposition routine will also return the determinant of A.Conventionally the lower triangular matrix L is defined such that all its diagonal elements are unity.This is what makes it possible to replace A with both L and U.Some libraries provide separate routines for the Gaussian Elimination and the Back–Substitutionsteps.
Often the Back–Substitution must be run separately for each column of B.Routines are provided for a wide range of different types of matrices. The symmetry of the matrix isnot often used however.0The time taken to solve equations in unknowns is generally proportional to for the Gaussian–Elimination (LU–decomposition) step. The Back–Substitution step goes as for each column of B. Asbefore this may be modified by parallelism.Matrix Algebra303.5.2 Iterative MethodsAs an alternative to the above library routines, especially when large sparse matrices are involved, it ispossible to solve the equations iteratively.The Jacobi Method(3.13)ð A Ò0D L Uwhere D is a diagonal matrix(i.e.) and L and U are strict lower and upper triangular21 Ä forforall Ò ). â Wematrices respectively (i.e.now write the Jacobi or Gauss–Jacobi iterativeprocedure to solve our system of equations asQS#T eT $" & QXD uBL U X v(3.14)QQ this procedure requires thewhere the superscripts refer to the iteration number.
Note that in practicestorage of the diagonal matrix, D, and a function to multiply the vector, X by L U.We first divide the matrix A into 3 parts+This algorithm resembles the iterative solution of hyperbolic (see section 2.3) or parabolic (see section 2.5) partial differential equations, and can be analysed in the same spirit. In particular care must betaken that the method is stable.Simple C code to implement this for a 1D Poisson’s equation is given below.int i;const int N = ??; // incomplete codedouble xa[N], xb[N], b[N];while ( ...
) // incomplete code{for ( i = 0; i < N; i++ )xa[i] = (b[i] - xb[i-1] - xb[i+1]) * 0.5;for ( i = 0; i < N; i++ )xb[i] = xa[i];}Note that 2 arrays are required for X and that the matrices, D, L and U don’t appear explicitely.The Gauss–Seidel MethodQQS#TQS=TAny implementation of the Jacobi Method (3.14) above will involve a loop over the matrix elements in eachcolumn of X. Instead of calculating a completely new matrix Xat each iteration and then replacingX with it before the next iteration, as in the above code, it might seem sensible to replace the elements ofX with those of Xas soon as they have been calculated.
![](https://s.studizba.com/z.php?f=/templates/si/i/noavatar.png&w=100&h=100&t=1"){kind=avatar}
