Numerical Linear Algebra (Direct Methods) ========================================= Numerical Linear Algebra (Direct Methods) (Mar 1994) Author: Jesus ROJO We are at present developing Packages which, once loaded, allow the development of a great part of the normal numerical routines. It is possible that some Packages may need others and this is the case with the Package we present, NUMERIAL, which carries out the most usual routines corresponding to the Direct Methods for solving linear systems. It requires the use of the NUMERICO Package that accompanies it. Let us first of all present these routines by forming an important Package for educational tasks. NUMERICO allows to use: CORounding[a,prec], which rounds off numbers in the usual way and CONorm1[X,n,prec], CONorm2[X,n,prec], CONormi[X,n,prec], that calculate the most usual vector norms. NUMERIAL allows us to use: ALFactor[A,b,n,prec,lu,pi,pp], ALCholeski[A,b,n,prec,pp], ALFactorwithRefinement[A,b,n,prec,lu,pi,no,tol,pp,pr], ALCholeskiwithRefinement[A,b,n,prec,no,tol,pp,pr], ALFactorwithnRefinements[A,b,n,prec,lu,pi,nref,pp,pr], ALCholeskiwithnRefinements[A,b,n,prec,nref,pp,pr], that is, LU and Choleski decomposition with possible iterative improvement. The initials CO and AL allow us to group the modules by Packages. File-package NUMERICO.M: instructions. -------------------------------------- CORounding[a,prec] gives the real number a rounded off to prec signifiant decimal places. Note: CORounding[3.55,2] gives 3.6 and CORounding[3.54999,2] gives 3.5. CORounding[-3.55,2] gives -3.6 and CORounding[-3.54999,2] gives -3.5. CONorm1[X,n,prec] gives the 1-norm of the vector X of size n rounded off to prec significant decimal places. CONorm2[X,n,prec] gives the euclidean norm (2-norm) of the vector X of size n rounded off to prec significant decimal places. CONormi[X,n,prec] gives the maximum norm (infinity norm) of the vector X of size n rounded off to prec significant decimal places. File-package NUMERIAL.M: instructions. -------------------------------------- ALFactor[A,b,n,prec,lu,pi,pp] decomposes the matrix A (with or without pivoting) in the form L U , and obtains the solution to the system A x = b . The integer n is the size of the square matrix A . The integer prec is the number of significant digits with which one wishes to work the method. The character lu is worth "l" or "u" and indicates the matrix (L or U) that conserves the diagonal; the other possesses a diagonal equal to 1. The string pi is worth "cp" or "sp" and indicates if the method is carried out with or without partial pivoting. The string pp is worth "dp" or "sp" and indicates if the method is carried out with or without partial double precision. The output will be printed on the screen and also in a file called LU.DAT. The matrices L and U , the matrix P.A or the rounded A (both under the name PA), the vector P of permutation of the rows and the solution X of the system remain in the memory. ALCholeski[A,b,n,prec,pp] decomposes the matrix A in the form L tL , and obtains the solution to the system A x = b . The method checks if the matrix A is symmetric and if it is positive definite. The integer n is the size of the square matrix A . The integer prec is the number of significant digits with which one wishes to work the method. The string pp is worth "dp" or "sp" and indicates if the method is carried out with or without partial double precision. The output will be printed on the screen and also in a file called CH.DAT. The matrix L , the rounded matrix A (under the name AR) and the solution X of the system remain in the memory. ALFactorwithRefinement[A,b,n,prec,lu,pi,no,tol,pp,pr] decomposes the matrix A (with or without pivoting) in the form L U , obtains the solution to the system A x = b and improves it through refinement until the relative difference is less than tol. The integer n is the size of the square matrix A . The integer prec is the number of significant digits with which one wishes to work the method. The character lu is worth "l" or "u" and indicates the matrix (L or U) that conserves the diagonal; the other possesses a diagonal equal to 1. The string pi is worth "cp" or "sp" and indicates if the method is carried out with or without partial pivoting. The real tol gives us the criteria for ending the iterative refinement. The string no is worth "n1" , "n2" or "ni" and indicates that the norm to be compared with tol is either the 1-norm, the euclidean norm or the maximum one. The string pp is worth "dp" or "sp" and indicates if the method is carried out with or without partial double precision. The string pr is worth "dr" or "sr" and indicates if the method is carried out with or without double precision in the refinement. The output will be printed on the screen and also in a file called LR.DAT. The matrices L and U , the matrix P.A or the rounded A (both under the name PA), the vector P of permutation of the rows, the solution X of the system and the successive refinements XR remain in the memory. ALCholeskiwithRefinement[A,b,n,prec,no,tol,pp,pr] decomposes the matrix A in the form L tL , obtains the solution to the system A x = b and improves it through refinement until the relative difference is less than tol. The method checks if the matrix A is symmetric and if it is positive definite. The integer n is the size of the square matrix A . The integer prec is the number of significant digits with which one wishes to work the method. The real tol gives us the criteria for ending the iterative refinement. The string no is worth "n1" , "n2" or "ni" and indicates that the norm to be compared with tol is either the 1-norm, the euclidean norm or the maximum one. The string pp is worth "dp" or "sp" and indicates if the method is carried out with or without partial double precision. The string pr is worth "dr" or "sr" and indicates if the method is carried out with or without double precision in the refinement. The output will be printed on the screen and also in a file called CR.DAT. The matrix L , the rounded matrix A (under the name AR), the solution X of the system and the successive refinements XR remain in the memory. ALFactorwithnRefinements[A,b,n,prec,lu,pi,nref,pp,pr] decomposes the matrix A (with or without pivoting) in the form L U , obtains the solution to the system A x = b and improves it through nref refinements. The integer n is the size of the square matrix A . The integer prec is the number of significant digits with which one wishes to work the method. The character lu is worth "l" or "u" and indicates the matrix (L or U) that conserves the diagonal; the other possesses a diagonal equal to 1. The string pi is worth "cp" or "sp" and indicates if the method is carried out with or without partial pivoting. The string pp is worth "dp" or "sp" and indicates if the method is carried out with or without partial double precision. The string pr is worth "dr" or "sr" and indicates if the method is carried out with or without double precision in the refinement. The output will be printed on the screen and also in a file called LN.DAT. The matrices L and U , the matrix P.A or the rounded A (both under the name PA), the vector P of permutation of the rows, the solution X of the system and the successive refinements XR remain in the memory. ALCholeskiwithnRefinements[A,b,n,prec,nref,pp,pr] decomposes the matrix A in the form L tL , obtains the solution to the system A x = b and improves it through nref refinements. The method checks if the matrix A is symmetric and if it is positive definite. The integer n is the size of the square matrix A . The integer prec is the number of significant digits with which one wishes to work the method. The string pp is worth "dp" or "sp" and indicates if the method is carried out with or without partial double precision. The string pr is worth "dr" or "sr" and indicates if the method is carried out with or without double precision in the refinement. The output will be printed on the screen and also in a file called CN.DAT. The matrix L , the rounded matrix A (under the name AR), the solution X of the system and the successive refinements XR remain in the memory. Use of these routines: some examples ----------------------------------- are given in the file NUMSAMPL.M that accompanies this set of Packages. The example uses each of the routines once. For instance, after the execution of ALCholeskiwithnRefinements[AC,b,2,3,3,"sp","sr"] which is the last line of NUMSAMPL.M, operations with the matrices L, AR, X and XR can be done. Thus, to check the magnitude of the residual of the solutin obtained, we can do Residual = b - AR . X , wich gives us {0.006, -0.002} . Example of output file: LN.DAT ------------------------------ For the matrix [ A | b ] = 2 1 1.222 4 2.1 3.666 of size 2 that is [ A | b ] = 2. 1. 1.22 4. 2.1 3.67 with 3 significant decimal places, the decomposition L U turns out to be (with partial double precision) (with diagonal in U and without pivoting) A = 2. 1. 4. 2.1 L = 1. 0 2. 1. U = 2. 1. 0 0.1 and the refined solution is X = {-4.97, 11.2} with refinements XR = -5.55 12.3 -5.69 12.6 -5.33 11.9 -4.97 11.2