We can label the coefficients and values of the right hand side with the following notation where the first subscript gives the equation number and the second gives the variable number that the coefficient multiplies.
We want to tuck the values of x, y, and z in a convenient location, so we introduce one more piece of notation. A vector (array) is introduced with x1=x, x2=y, and x3=z , and we get the final disguised but useful equations:
In Fortran we introduce arrays "A" and "B" to store the coefficients and right hand sides of the
equations, and "X" for the solution.
REAL A(3,3), B(3), X(3)
In terms of the Fortran arrays our equations would be:
A(1,1)*X(1) + A(1,2)*X(2) + A(1,3)*X(3) = B(1)
A(2,1)*X(1) + A(2,2)*X(2) + A(2,3)*X(3) = B(2)
A(3,1)*X(1) + A(3,2)*X(2) + A(3,3)*X(3) = B(3)
Please don't confuse the three equations I've just written with Fortran assignment statements.
They are not valid Fortran, just a way of indicating the way to translate from standard equation
notation, to Fortran array notation.
The work of solving the system of equations should be performed in a dedicated subroutine. As
with sorting and interpolation, the best place to get such subroutines is in NETLIB
(http://www.netlib.org). Within NETLIB, linear algebra subroutines are found in a library called
LINPACK (http://www.netlib.org/linpack ) or in a more recent version called LAPACK
(http://www.netlib.org/lapack). A Fortran 90 version of LAPACK is under development.
In Linear Algebra, the notation for my equations is further simplified to Ax=b. Matrix "A" is multiplied by the column vector "x", and the resulting column vector is equal to the column vector "b". This can also be written as:
This notation combined with the original linear equations provides a definition of multiplication of an column vector by a matrix. It can be extended to provide a definition for the general multiplication of a matrix by another matrix. For example lets look at how to get a final matrix A that is the product of 2 matrices L and U (A = LU). For now take the following definitions for L and U
The general definition for the elements resulting from the matrix product is:
where n is the number of rows in u (number of unknowns). Here n=3. In plainer terms the product can be written as a combination of matrix vector products. The first column of A is the result of multiplying L by the first column in U:
The second column of A is the product of L and the second column of U:
The third column of A is the product of L and the third column of U:
Combining the columns gives:
For those of you not totally up to speed on linear algebra you should go through the above exercise
to calculate the product UL. You will discover that UL is not equal to LU. Matrix multiplication
is not commutative, like regular scalar multiplication. We are not guaranteed that the product UL
is equal to the product LU
In Fortran, the above matrix product can be accomplished with the following subroutine.
subroutine matmult(a,b,n,c) dimension a(n,n),b(n,n),c(n,n) do 100 i=1,n do 100 j=1,n c(i,j)=0. do 100 k=1,n 100 c(i,j)=c(i,j)+a(i,k)*b(k,j) return end
Using this type subroutine or loop structure is not a good idea when Fortran 90 is available, because you have a matrix multiply available that almost always will be implemented more efficiently than any compilation of your Fortran. For arrays established by a "REAL A(3,3), L(3,3), U(3,3)" you do the multiply with the following line:
A = MATMUL(L,U)
Time savings associated with things like matrix multiplies can be a very big deal in many applications. If you look at the Fortran matrix multiply above you will see that the operation takes n3 multiplies and the same number of adds. You don't have to get to very large matrices or very frequent use of the matrix multiply before you are looking at significant savings in computer time by using a optimally coded function like MATMUL.