(* :Title:StrangLinearAlgebra *) (* :Author: Allan hayes, hay@le.ac.uk hay@haystack.demon.co.uk (NeXTmail) *) (* :Summary: StrangLinearAlgebra is a package of functions intended for use with the textbook, Introduction to Linear Algebra, Gilbert Strang, Wellesley-Cambridge Press,Box 812060, Wellesley, MA 02181 1993." *) (* :Context:StrangLinearAlgebra` *) (* :Package Version: 0.1*) (* :Copyright: Copyright 1993, 1994 Allan Hayes. Permission is hereby granted to modify and/or make copies of this file for any purpose other than direct profit, or as part of a commercial product, provided this copyright notice is left intact. Sale, other than for the cost of media, is prohibited. Permission is hereby granted to reproduce part or all of this file, provided that the source is acknowledged. *) (* :History: MATLAB code by Cleve Moler (1993)*) (* :Keywords: Linear Algebra*) (* :Source:*) (* :Warning: Extends the definition of Transpose as follows, For a vector v Transpose[v] = v . *) (* :Mathematica Version: 2.2 *) (* :Requirement: *) (* :Limitation: The default setting of the option ZeroTest is ZeroTest ->(Chop[N[Together[#, Trig->True] == 0]]&) *) (* :Discussion: Design decisions Names: names are all lower case, as in the text; built-in Mathematica functions are therefore recognizable as beginning with an upper case letter. Input and Output:these are almost all of the same pattern as in the Matlab version referred to in the text ; any differences are pointed out in the usage notes for the individual functions. Faulty Inputs and Failure of Computation: these should give an error message and cause the input to be returned unaltered. Explanation: most of the functions should have the option Explain which when set to True will result in an explanation of the computation being printed. In some cases Explain also has a setting PatternOnly, which will cause the matrices used to be shown with 0, 0.0 and 1, 1.0 to be shown as they are and all other entries to be shown as +. ZeroTest: most of the function have a ZeroTest option like many built-in functions. Simplification of expressions:how much simplification of expressions should be done and how much shown to the user? Sometimes simplification must be used; for example, to decide if an expression is zero. Even here there is a question as to whether an expression should be actually replaced by the simplified form . Some arguments for and against non essential simplification are: for, -- the user should be given understandable output; -- the computation may be speeded up by subsequent steps having simpler input; against, -- simplification can itself take a lot of time. -- the user may want to retain the forms of some expressions and be left to decide what kind of simplification is appropriate . In the present version the use and showing of of simplification is kept to a minimum. It is only used in deciding whether an expression in zero and the expression is replaced by the simplified form only if this is zero. I shall probably use and show more simplification in a later version. COMMENTS ON THIS WILL BE VERY WElCOME. *) BeginPackage[ "StrangLinearAlgebra`", "Graphics`Colors`", "Graphics`Arrow`", "Utilities`FilterOptions`" ]; (*Clear previous definitions made in this context:*) Unprotect["`*"]; ClearAll["`*"]; (*Usage messages:*) StrangLinearAlgebraInfo::usage = " StrangLinearAlgebra is package of functions intended for use with\n the textbook\n\n Introduction to Linear Algebra\n Gilbert Strang\n\n Wellesley-Cambridge Press\n Box 812060, Wellesley, MA 02181\n 1993\n \n The function defined are listed below, fuller information\n may be found under their separate entries.\n \n basis Basis for the column space;\n cofactor The matrix of cofactors;\n cramer Solution by Cramer©s method;\n determ Determinant (uses plu);\n eigen Eigenvalues and eigenvectors;\n eigen2 Two by two eigenvalues and eigenvectors;\n findpiv First pivot (non-zero entry) in submatrix;\n grams Gram-Schmidt orthogonalization;\n invert Matrix inversion by Gauss-Jordan elimination;\n linefit Plot the least squares fit by a line;\n null Nullspace of a matrix;\n permdet Determinant of a permutation;\n plu Pivoting, rectangular, LU factorization;\n projmat Projection matrix onto the column space;\n randperm Random permutation;\n ref Reduced row exchelon form;\n signperm Sign of a permutation;\n slu Simple, square, LU factorization;\n slv Simple linear equation solver;\n solve Particular solution;\n splu Square LU factorization with row exchanges;\n splv The solution to a square, invertible system. "; basis::usage = "basis[m] gives a matrix which has a basis for the column space of the matrix m as its columns;\n\n basis[m, Explain->True] adds some explanation;\n basis[m, ZeroTest->test] uses test to test for objects being zero; the default setting is\n ZeroTest ->(Chop[N[Together[#, Trig->True] == 0]]&)\n\n Examples:\n basis[{{1,2,3,4},{2,1,3,1},{1,1,2,3}}]\n basis[{{1,2,3,4},{2,1,3,1},{1,1,2,3}}, Explain->True] "; cofactor::usage = " cofactor[m] gives the matrix of cofactors of the square matrix; m.\n\n cofactor[m, Explain->True] adds some explanation;\n cofactor[m, Explain->PatternOnly] explains using only the the patterns of the matrices;\n cofactor[m, ZeroTest->test] uses test to test for objects being zero; the default setting is\n ZeroTest ->(Chop[N[Together[#, Trig->True] == 0]]&)\n\n Examples:\n cofactor[{{1,2,3},{2,-1,3},{1,1,2}}]\n cofactor[{{1,2,3},{2,-1,3},{1,1,2}}, Explain->True]\n cofactor[{{1,2,3},{2,-1,3},{1,1,2}}, Explain->PatternOnly] "; cramer::usage = " cramer[m,b] tries to solve m.x = b by Cramer©s method;\n\n cramer[m,b, Explain->True] adds some explanation;\n cramer[m,b, Explain->PatternOnly] explains using only the the patterns of the matrices;\n cramer[m,b, ZeroTest->test] uses test to test for objects being zero; the default setting is\n ZeroTest ->(Chop[N[Together[#, Trig->True] == 0]]&)\n\n Examples:\n cramer[{{1,2,3},{2,-1,3},{1,1,2}}, {1,2,3}]\n cramer[{{1,2,3},{2,-1,3},{1,1,2}}, {1,2,3}, Explain->True]\n\n See also: LinearSolve, Solve, Reduce, NSolve "; determ::usage = " determ[m] gives the matrix of cofactors of the matrix m;\n determ[m, Explain->True] adds some explanation;\n determ[m, Explain->PatternOnly] explains using only the the pattern of the matrix;\n determ[m, ZeroTest->test] uses test to test for objects being zero; the default setting is\n ZeroTest ->(Chop[N[Together[#, Trig->True] == 0]]&)\n\n Example:\n determ[{{1,2,3},{2,-1,3},{1,1,2}}]\n determ[{{1,2,3},{2,-1,3},{1,1,2}},Explain->True]\n determ[{{1,2,3},{2,-1,3},{1,1,2}},Explain->PatternOnly]\n\n See also: Det. "; eigen::usage = " eigen[m], for a square matrix m, gives {tr, det, pol, vals, vecs}, where\n tr is its trace;\n det is its determinant;\n pol is its characteristic polynomial\n vals is its eigenvalues\n vecs is the list of corresponding eigenvectors\n\n eigen[m, Explain->True] adds some explanation;\n eigen[m, ZeroTest->test] uses test to test for objects being zero;\n the default setting is\n ZeroTest ->(Chop[N[Together[#, Trig->True] == 0]]&)\n\n Examples:\n eigen[{{1,0,3},{2,-1,3},{1,1,2}}]\n eigen[{{1,0,3},{2,-1,3},{1,1,2}}, Explain->True]\n\n See also: Eigenvalues, Eigenvalues, Eigensystem. "; eigen2::usage = " eigen[m], for a 2 by 2 matrix m, gives {vals, vecs} where vals is its eigenvalues\n vecs is the list of corresponding eigenvectors.\n\n eigen[m, Explain->True] adds some explanation;\n eigen[m, ZeroTest->test] uses test to test for objects being zero; the default setting is\n ZeroTest ->(Chop[N[Together[#, Trig->True] == 0]]&)\n\n Examples:\n eigen2[{{1,2},{1,1}}]\n eigen2[{{1,2},{1,1}}, Explain->True]\n\n See also: eigen, Eigenvalues, Eigenvalues, Eigensystem. "; findpiv::usage = " findpiv[m, {rmin, rmax}, {cmin,cmax}] gives {pos, newm} where\n pos is the position of non-zero entry found by searching\n column by column in the submatrix of m common to rows\n rmin to rmax and columns cmin to cmax;\n newm is the final form of m after the search (if an entry\n simplifies to zero, it is replaced by zero.\n\n findpiv[m, r, {cmin,cmax}] =\n findpiv[m, {r, number of rows in m}, {cmin,cmax}];\n similarly for findpiv[m, r, c] and findpiv[m, {rmin, rmax}, c]\n\n findpiv[m, Explain->True] adds some explanation;\n findpiv[m, Explain->PatternOnly] explains using only the the pattern of the matrix;\n findpiv[m, ZeroTest->test] uses test to test for objects being zero; the default setting is\n ZeroTest ->(Chop[N[Together[#, Trig->True] == 0]]&)\n\n Examples:\n findpiv[{{1,2,3,4},{0,0,0,2},{0,0,1,1}},2,2]\n findpiv[{{1,2,3,4},{0,0,0,2},{0,0,1,1}},2,2, Explain->True]\n findpiv[{{1,2,3,4},{0,0,0,2},{0,0,1,1}},2,2, Explain->PatternOnly]\n findpiv[{{1,2,3,4},{0,0,0,0},{2,0,0,0}},2,2, Explain->True] "; grams::usage = " grams[m] give a triple {Q, R, ip} where Q and R are matrices\n Q has orthonormal columns that span the column space of m.\n R is upper triangular matrix\n Q.R = m\n ip is the inner product being used - the default is\n ip = Transpose[#1].#2& for real m;\n ip = Conjugate[Transpose[#1]].#2& for complex m. \n\n grams[m, Explain->True] adds some explanation;\n grams[m, Normalized->False] does not normalize the columns in the output, leaving them orthogonal but not necessarily orthonormal;\n grams[m, InnerProduct -> ip] uses the inner product ip; grams[m, ZeroTest->test] uses test to check for objects being zero; the default setting is \n ZeroTest ->(Chop[N[Together[#, Trig->True] == 0]]&)\n \n Examples:\n grams[{{1,2,2,3},{3,6,2,1}}]\n grams[{{1,2,2,3},{3,6,2,1}}, Explain -> True]\n grams[{{1,2,2,3},{3,6,2,1}}, Normalized -> False,Explain -> True]\n grams[{{1,2,2,3},{3,6,2,1}},\n InnerProduct -> (Transpose[#1].{{1,1},{1,5}}.#2&),\n Explain -> True\n ]\n (*The position of the brackets, (...&), in the last one is important.*)\n\n Warning\n This package uses the extension of function Transpose:\n Transpose[v] = v when v is a vector that is automatically loaded by the package.\n\n See also:\n LinearAlgebra`Orthogonalization` "; invert::usage = " invert[m] gives the inverse of the square matrix m.\n\n invert[m, Explain->True] adds some explanation;\n grams[m, ZeroTest->test] uses test to test for objects being zero; the default setting is\n ZeroTest ->(Chop[N[Together[#, Trig->True] == 0]]&)\n\n Examples:\n invert[{{1,-2,2},{0,-1,1},{2,1,0}}]\n invert[{{1,-2,2},{0,-1,1},{2,1,0}}, Explain->True]\n invert[{{1,-2,2},{0,-1,1},{2,1,-1}}, Explain->True]\n Clear[a];invert[{{1,-2,2},{0,-1,1},{2,1,a}}, Explain->True]\n\n See also:Inverse "; linefit::usage = " line[xs,vals,x] calculates and plots the best fit straight line to the data given by values vals at numbers xs; the formula for this is given as a formula in the variable x.\n The output is this formula.\n\n invert[m, Explain->True] adds some explanation.\n Options for Plot may be used to modify the plot.\n\n Examples:\n xs = {0.07, 0.59, 1.33, 1.86, 2.49, 3.12, 3.86, 4.34}\n vals = {0.11, 0.85, 0.77, 2.14, 1.63, 2.27, 3.83, 2.50}\n linefit[xs,vals,x]\n linefit[xs,vals,x, Explain -> True]\n linefit[xs,vals,x,\n Explain -> True, Frame -> True, Background -> Eggshell ]\n\n See also: Fit, Statistics`NonlinearFit`, NumericalMath`PolynomialFit`, NumericalMath`SplineFit`, InterpolatingFunction. "; lsq::usage = " lsq[A,b] calculates the least square solution to the equation\n A.x = b. The output is {x,p,e} where x is th solution, p = A.x\n and e is the error p-b.\n\n lsq[A, Explain->True] adds some explanation.\n\n Example:\n lsq[{{2,1,1},{-1,1,2},{1,2,3}}, {2,3,1}]\n lsq[{{2,1,1},{-1,1,2},{1,2,3}}, {2,3,1}, Explain->True]\n\n See also: Fit, Statistics`NonlinearFit, NumericalMath`PolynomialFit, NumericalMath`SplineFit, InterpolatingFunction. "; null::usage = " null[A] gives a matrix with columns forming a basis for the null space of the matrix A.\n\n null[A, Explain->True] adds some explanation;\n null[A, ZeroTest->test] uses test to test for objects being zero;\n the default setting is\n ZeroTest ->(Chop[N[Together[#, Trig->True] == 0]]&)\n\n Examples:\n null[{{0,1,-1,0,2},{1,2,-1,1,3},{1,1,0,1,0}}]\n null[{{0,1,-1,0,2},{1,2,-1,1,3},{1,1,0,1,0}},Explain->True]\n\n See also: NullSpace (but note that NullSpace give a matrix with ROWS, not columns, forming a basis for the null space of A) "; permdet::usage = " permdet[perm] gives a matrix of the permutation perm.\n\n permdet[m, Explain->True] adds some explanation.\n\n Examples:\n permdet[{1,3,4,2}]\n permdet[{1,3,4,2}, Explain -> True]\n\n See also: signperm, randperm, Sort, Permutations, Signature. "; plu::usage = " plu[A], for the matrix A, finds {P,L,U,pivcol,sign}, where,\n if A is m by n,\n P is an m by m permutation matrix, recording the swaps used;\n L is m by m lower triangular with ones on the diagonal;\n U is m by n upper triangular(the same size as A);\n P.A = L.U;\n pivcol is the list of the positions of the pivot columns in U;\n sign is the sign of the permutation represented by P;\n\n plu[A, Explain->True] adds some explanation;\n plu[A, Explain->PatternOnly] explains using only the the pattern of the matrix;\n plu[A, ZeroTest->test] uses test to test for objects being zero; the default setting is\n ZeroTest ->(Chop[N[Together[#, Trig->True] == 0]]&)\n\n Examples:\n plu[{{1,2,3,4},{2,4,2,1},{-3,1,1,3}}]\n plu[{{1,2,3,4},{2,4,2,1},{-3,1,1,3}}, Explain->True]\n plu[{{1,2,3,4},{2,4,2,1},{-3,1,1,3}}, Explain->PatternOnly]\n\n See also: LUDecomposition, RowReduce. "; projmat::usage = " projmat[A] for an m by n matrix A gives a matrix P such\n that for each vector b of length m, P.b is the projection of b onto the column space of A.\n\n projmat[A, Explain->True] adds some explanation;\n projmat[A, ZeroTest->test] uses test to test for objects being zero; the default setting is\n ZeroTest ->(Chop[N[Together[#, Trig->True] == 0]]&)\n\n Examples:\n projmat[{{-1,-1}, {-1,-1}, {-1,2}}]\n projmat[{{-1,-1}, {-1,-1}, {-1,2}}, Explain->True] "; randperm::usage = " randperm[n], for a non-negative integer gives a random permutation of the integers 1 to n (this is the empty set if n = 0).\n\n randperm[n, Explain->True] adds some explanation.\n\n Examples:\n randperm[10]\n randperm[10, Explain->True]\n\n See also: signperm, permdet, Sort, Permutations, Signature. "; rats::usage = " rats[x] replaces real numbers in the expression x with rational numbers that differ from them by not more than 10^-6.\n\n Example\n rats[1.235]\n\n See also: Rationalize, Chop, Round, Floor, Ceiling, LatticeReduce. "; ref::usage = " ref[A], for a matrix A, gives {R, pivcols} where R is the reduced row echelon form of the matrix A, and pivcols is the list of the positions of the pivot columns of R.\n\n ref[A, Explain->True] adds some explanation;\n ref[A, Explain->PatternOnly] explains using only the the pattern of the matrix;\n ref[A, ZeroTest->test] uses test to test for objects being zero; the default setting is\n ZeroTest ->(Chop[N[Together[#, Trig->True] == 0]]&)\n\n Examples:\n ref[{{0,0,1,1,2},{-2,2,0,0,-2},{-2,2,2,1,2},{0,0,-2,-1,-4}}]\n ref[{{0,0,1,1,2},{-2,2,0,0,-2},{-2,2,2,1,2},{0,0,-2,-1,-4}},\n Explain->True\n ]\n ref[{{0,0,1,1,2},{-2,2,0,0,-2},{-2,2,2,1,2},{0,0,-2,-1,-4}},\n Explain->PatternOnly\n ]\n\n See also: RowReduce. "; signperm::usage = " signperm[perm], where perm is a permutation of an initial segment of the positive integers, gives the sign of perm.\n\n signperm[perm, Explain->True] adds some explanation;\n\n Examples:\n signperm[{3, 1, 4, 5, 2}]\n signperm[{3, 1, 4, 5, 2}, Explain->True]\n\n See also: Signature, detperm,randperm, Permutations, Order, Sort. "; slu::usage = " slu[A] tries to find {L,U}, for a square matrix A, where,\n if A is m by m, then\n L is m by m, lower triangular with ones on the diagonal;\n U is m by m, upper triangular and non-zero on the diagonal;\n A = L.U;\n The algorithm fails if A is singular.\n\n slu[A, Explain->True] adds some explanation;\n slu[A, Explain->PatternOnly] explains using only the the pattern of the matrix;\n slu[A, ZeroTest->test] uses test to test for objects being zero; the default setting is\n ZeroTest ->(Chop[N[Together[#, Trig->True] == 0]]&)\n\n Examples:\n slu[{{1,2,3},{2,3,2},{-3,1,1}}]\n slu[{{1,2,3},{2,3,2},{-3,1,1}}, Explain->True]\n slu[{{1,2,3},{2,3,2},{-3,1,1}}, Explain->PatternOnly]\n slu[{{1,2,3},{2,4,2},{-3,1,1}}, Explain->True]\n Clear[a];slu[{{1,2,3},{2,a,2},{-3,1,1}}, Explain->True]\n\n See also: splu, plu. "; slv::usage = " slv[A,b] tries to solve the equation A.x = b by using the LU factorization of A calculated by the function slu. This method will fail if A is singular.\n If it succeeds the output is the list s such that A.s = b.\n\n slv[A, Explain->True] adds some explanation;\n slv[A, Explain->PatternOnly] explains using only the the pattern of the matrix;\n slv[A, ZeroTest->test] uses test to test for objects being zero; the default setting is\n ZeroTest ->(Chop[N[Together[#, Trig->True] == 0]]&)\n\n Examples:\n slv[{{1,2,3},{2,3,2},{-3,1,1}},{2,1,3}]\n slv[{{1,2,3},{2,3,2},{-3,1,1}},{2,1,3}, Explain->True]\n slv[{{1,2,3},{2,3,2},{-3,1,1}},{2,1,3}, Explain->PatternOnly]\n slv[{{1,2,3},{2,4,2},{-3,1,1}},{2,1,3}, Explain->True]\n Clear[a];slv[{{1,2,3},{2,a,2},{-3,1,1}},{2,1,3}, Explain->True]\n\n See also: splv,solve, LinearSolve, Solve, Reduce,NSolve. "; solve::usage = " solve[A,b] finds a particular solution of the equation A.x = b, if one exists, by using the PLU factorization of A calculated by the function plu.\n If it succeeds the output is a list s such that A.s = b.\n \n solve[A, Explain->True] adds some explanation;\n solve[A, Explain->PatternOnly] explains using only the the pattern of the matrix;\n solve[A, ZeroTest->test] uses test to test for objects being zero; the default setting is\n ZeroTest ->(Chop[N[Together[#, Trig->True] == 0]]&)\n\n Examples:\n solve[{{1,2,3},{2,4,2},{-3,-6,1}},{4,6,-7}]\n solve[{{1,2,3},{2,4,2},{-3,-6,1}},{4,6,-7}, Explain->True]\n solve[{{1,2,3},{2,4,2},{-3,-6,1}},{4,6,-7}, Explain->PatternOnly]\n solve[{{1,2,3},{2,4,2},{-3,-6,1}},{4,6,0}, Explain->True]\n Clear[a];\n solve[{{1,2,3},{2,4,2},{-3,-6,1}},{4,6,a}, Explain->True]\n\n See also; LinearSolve, Solve, Reduce, NSolve. "; splu::usage = " splu[A], for square matrix A, tries to find {P,L,U,sign}, where,\n P is a permutation matrix, recording the swaps used;\n L is square lower triangular with ones on the diagonal;\n U is a upper triangular (the same size as A);\n P.A = L.U;\n sign is the sign of the permutation represented by P;\n The algorithm fails if A is singular\n\n splu[A, Explain->True] adds some explanation;\n splu[A, Explain->PatternOnly] explains using only the the pattern of the matrix;\n splu[A, ZeroTest->test] uses test to test for objects being zero; the default setting is\n ZeroTest ->(Chop[N[Together[#, Trig->True] == 0]]&)\n\n Examples:\n splu[{{1,2,3},{2,4,2},{-3,1,1}}]\n splu[{{1,2,3},{2,4,2},{-3,1,1}}, Explain->True]\n splu[{{1,2,3},{2,4,2},{-3,1,1}}, Explain->PatternOnly]\n splu[{{1,2,3},{2,4,2},{-3,-6,1}}, Explain->True]\n Clear[a];splu[{{1,2,3},{2,4,2},{-3,a,1}}, Explain->True]\n\n See also: plu. "; splv::usage = " splv[A,b], for a square matrix A, tries to solve the equation A.x = b by using the factorization of A calculated by the function splu. This method will fail if A is singular.\n If it succeeds the output is the list s such that A.s = b.\n\n splv[A, Explain->True] adds some explanation;\n splv[A, Explain->PatternOnly] explains using only the the pattern of the matrix;\n splv[A, ZeroTest->test] uses test to test for objects being zero; the default setting is\n ZeroTest ->(Chop[N[Together[#, Trig->True] == 0]]&)\n\n Examples:\n splv[{{1,2,3},{2,4,2},{-3,1,1}},{2,3,1}]\n splv[{{1,2,3},{2,4,2},{-3,1,1}},{2,3,1}, Explain->True]\n splv[{{1,2,3},{2,4,2},{-3,1,1}},{2,3,1}, Explain->PatternOnly]\n splv[{{1,2,3},{2,4,2},{-3,-6,1}},{2,3,1}, Explain->True]\n Clear[a];splv[{{1,2,3},{2,4,2},{-3,a,1}},{2,3,1}, Explain->True]\n\n See also: solve, LinearSolve, Solve, Reduce, NSolve. "; System`Explain::usage = " Explain is an option for many functions in AHPackages;\n Explain -> True causes an explanation of the evaluation of the function to be generated;\n Explain ->PatternOnly causes the pattern of the evaluation to be generated.\n The default is Explain -> False. "; InnerProduct::usage = "InnerProduct is an option for functions concerned with inner products; it allows the user to set the inner product to be used. The default is InnerProduct -> Automatic which sets the inner product to\n Transpose[#1].#2& in the real situation;\n Conjugate[Transpose[#1]].#2& in the complex case.\n The inner product must be given as a pure binary function; one pattern\n is ip = (#1.M.#2)&, where M is a poitive definite matrix. "; Normalized::usage = "Normalized is an option for functions concerned with inner products. Normalized -> True causes the vectors in the answer to be normalized; Normalized -> False leaves the vectors un normalized. "; System`PatternOnly::usage = "PatternOnly is a value for the option Explain that causes the explanation to only show the pattern of the calculation "; StrangLinearAlgebra::usage = "StrangLinearAlgebra is a name for the messages used in this package"; (* Optioms *) Options[basis] = {Explain -> False, ZeroTest ->(Chop[N[Together[#, Trig->True] == 0]]&) }; Options[cofactor] = Options[basis]; Options[cramer] = Options[basis]; Options[determ] = Options[basis]; Options[eigen] = Options[basis]; Options[eigen2] = Options[basis]; Options[findpiv] = Options[basis];; Options[grams] = { Normalized ->True, Explain->False, ZeroTest ->(Chop[N[Together[#, Trig->True] == 0]]&), InnerProduct -> Automatic }; Options[invert] = Options[basis]; Options[linefit] = Prepend[Options[Plot], Explain->False]; Options[lsq] = {Explain->False}; Options[null] = Options[basis]; Options[permdet] = {Explain -> False}; Options[plu] = Options[basis]; Options[projmat] = Options[basis]; Options[randperm] = {Explain -> False}; Options[ref] = Options[basis]; Options[signperm] = {Explain -> False}; Options[slu] = Options[basis]; Options[slv] = Options[basis]; Options[splu] = Options[basis]; Options[splv] = Options[basis]; Options[solve] = Options[basis]; Begin["`Private`"]; (* Messages *) StrangLinearAlgebra::notperm = "Argument `2` of `1` is not a permutation; it does not contain `3` although its length is `4`."; StrangLinearAlgebra::notpermM = "Argument `2` of `1` is not a permutation; its determinant is `3` not -1 or +1."; StrangLinearAlgebra::notposintlst = "Argument `2` of `` must be a list of positive integers."; StrangLinearAlgebra::notnzreallst = "Argument `2` of `1` must be a non-empty list of non-complex numbers."; StrangLinearAlgebra::badlnth = "Argument `2` of `1` must have `3` elements (the same number as argument `4`), not `5`."; StrangLinearAlgebra::notsym = "Argument `2` of `1` is expected to be a symbol, not `3`."; StrangLinearAlgebra::notmat = "Argument `2` of `1` must be a matrix."; StrangLinearAlgebra::notsqrmat = "Argument `2` of `1` must be a square matrix."; StrangLinearAlgebra::notnnint = "Argument `2` of `1` must be a non-negative integer not `3`."; StrangLinearAlgebra::notsqrmattry = "Argument `2` of `1` must be a square matrix; you might like to try plu."; StrangLinearAlgebra::rhs = "Argument 2 of `` must be a vector with `` entries: the number of rows of Argument 1."; StrangLinearAlgebra::rhs1 = "Argument 2 of `` must be a vector with one entry: the number of rows of Argument 1."; StrangLinearAlgebra::singmat = "Argument 1 of ``, is a singular matrix."; StrangLinearAlgebra::fldalgo = "The algorithm used in ``, has failed; you might like to try ``."; StrangLinearAlgebra::lindepcols = "The first `2` columns of argument `3` of `1` are linearly dependent."; StrangLinearAlgebra::zerocol = "The first column of argument `2` of `1` is a zero column."; StrangLinearAlgebra::notnl = "Argument `2` of `1` must be an integer, n, with `3` <= n <= `4` or a list {r,s} with `3` <= r <= s <= `4`; not `5`."; StrangLinearAlgebra::notzero = "Argument `2` of `1` must be 0 when the matrix in Argument 1 has no columns; not `3`."; StrangLinearAlgebra::nopiv = "`` found no pivots b."; StrangLinearAlgebra::zeropiv = "`1` has found a zero pivot at position (`2`, `2`)."; StrangLinearAlgebra::notopt = "Options expected (instead of `2`) beyond position `3` in `1`. An option must be a rule or a list of rules."; StrangLinearAlgebra::incons = "The input to `1` represents an inconsistent system of equations."; StrangLinearAlgebra::not2by2 = "Argument 1 of `1` must be a 2 by 2 matrix."; StrangLinearAlgebra::nocols = "The columns of the matrix, Argument 1 of `1`, are, or simplify to, zero. There seems to be no solution to the problem of the form required. The output returned is {Q,1,ip}, where Q is the `2` by 0 matrix, and ip is the inner product used: Q is in fact an orthonormal basis for the column space of the input matrix."; (*Messages and abbreviations*) msg[n_, v___] := Message[MessageName[StrangLinearAlgebra,n],fn,v]; thfld := Throw[$Failed]; fld = $Failed; (*Removing any unique symbols created*) removeUnique[a_String] := Remove@@ Array[StringJoin[a,ToString[#]]&, StringDrop[ToString[Unique[a]],StringLength[a]]//ToExpression ]; removeUnique[a___String]:= Scan[removeUnique,{a}]; (*Testing for and replacing by zero: tozeroQ[test,changed] gives a function, tz say, which for an expression x which is either a symbol or of the form c[[...]] returns True when test[x] is True and otherwise returns False; AND has the SIDE EFFECT: if x is not zero sets x to 0 and test[x] is true then changed is set to True; recording that an apparently non-zero expression has been simplified to and replaced by zero. tz has attribute HoldFirst so that it passes the expression x, not its value. If changed is to be passed between functions it must not be localized in their definitions by Block or Module etc*) Attributes[tozeroQ] = {HoldRest}; tozeroQ[test_,changed_] = Function[x, If[test[x], If[!TrueQ[x==0], If[FreeQ[x,Real],x = 0,x = 0.0]; changed = True]; True, False, False ], {HoldAll} ]; (*Extension of Transpose*) Unprotect[Transpose]; Transpose[v_?VectorQ] := v; Protect[Transpose]; (*Formatting printed ouput.*) printwide[e__] := Module[{pw}, pw = PageWidth/.Options[$Output]; SetOptions[$Output, PageWidth -> Infinity]; Print[e]; SetOptions[$Output, PageWidth -> pw]; ]; matform[opts___?OptionQ] := Which[ MatchQ[#,{{}..}], #1, MatrixQ[#1, (IntegerQ[#]||Head[#]==String)&], TableForm[#1, opts, TableAlignments -> {Center, Right}], MatrixQ[#1], TableForm[#1, opts, TableAlignments -> {Center, Center}], True, #1 ] &; (* basis *) basis[e___]:= With[{result = Block[{fn = basis},basisX[e]]}, result/;result =!= fld ]; basisX[A_,opts___] := Module[{ponopt,n,expln,zero,Pr,PrW,U,pivcol}, Which[ !MatrixQ[A], msg["notmat",1];thfld, (posnopt = Position[{opts}, a_/;!OptionQ[a],{1},1,Heads->False])=!= {}, msg["notopt", {opts}[[posnopt[[1,1]]]], 1];thfld ]; n = Dimensions[A][[2]]; (*Read options*) {expln, zero} = {Explain,ZeroTest}/.{opts}/.Options[basis]; (*If an explanation is required, turn on printing.*) If[expln, Pr = Print; PrW = printwide, Pr = PrW = Hold ]; Pr["We want to find a matrix whose columns form a basis for the"]; Pr["column space of the matrix"]; Pr[""]; PrW[" ", A//matform[]]; Pr[""]; Pr["Call this matrix A."]; Pr["The U part of the PLU factorization of A, with the"]; Pr["pivot columns headed \"p\", is"]; {U,pivcol} = plu[A, ZeroTest->zero][[{3,4}]]; Pr[""]; PrW[" ", U//matform[TableHeadings -> {None,MapAt["p\n"&,Table["",{n}], List/@pivcol]}] ]; Pr[""]; Pr["The pivot columns clearly form a basis for the column space"]; Pr["of U; so, since this was produced from A by row operations,"]; Pr["and row operations do not alter linear relations amongst the"]; Pr["columns, we see that the corresponding columns of A form a"]; Pr["basis for the column space of A"]; Pr[""]; PrW[" ", A//matform[TableHeadings -> {None,MapAt["p\n"&,Table["",{n}], List/@pivcol]}] ]; ans = (#[[pivcol]])&/@A; Pr[""]; Pr["So, the columns of the matrix below give a basis as required."] Pr[""]; PrW[" ", ans//matform[]]; Pr[""]; Pr["This matrix is returned as output."]; ans ]//Catch; (* Trap and warn of any other faulty inputs*) basisX[e___] := (msg["argm",Length[{e}],1];fld); (* cofactor *) cofactor[e___] := With[{result = Block[{fn = cofactor},cofactorX[e] ]}, result/;result =!=fld ]; cofactorX[A_,opts___]:= Module[{posnopt,expln,m,n,Pr,PrW,pat,r,c}, Which[ !MatrixQ[A]||!SameQ@@Dimensions[A], msg["notsqrmat", 1];thfld, (posnopt = Position[{opts}, a_/;!OptionQ[a],{1},1,Heads->False])=!= {}, msg["notopt", {opts}[[posnopt[[1,1]]]], 1];thfld ]; (*Read options.*) expln = Explain/.{opts}/.Options[cofactor]; {m,n} = Dimensions[A]; (*Turn on explanation if required.*) If[MatchQ[expln,True|PatternOnly], Pr = Print;PrW = printwide, Pr = PrW = Hold ]; If[expln === PatternOnly, pat = Map[If[MatchQ[#,0|0.|1|1.0],#,"+"]&,#,{2}]&, pat = Identity ]; (*Begin explanation*) Pr["Given the square matrix", If[expln === PatternOnly, " with pattern", "" ] ]; Pr[""]; PrW[" ", pat[A]//matform[] ]; Pr[""]; Pr["(call it A)"]; Pr["we want to find its cofactor matrix."]; Pr["This is the matrix C with the same dimensions as A and with"]; Pr[""]; Pr[" ", "C(i,j) = (-1)"^"i+j", "Det[M[i,j]]"]; Pr[""]; Pr["where M[i,j] is the submatrix of A found by removing "]; Pr["row i and column j"]; Pr["An important property is"]; Pr[""] Pr[" ", "C"^"T",".A"," = ","A.C"^"T"," = det[A].I."]; Pr[""]; r = Range[m]; c = Table[ (-1)^(i+j) Det[A[[Drop[r,{i}],Drop[r,{j}]]]], {i,m},{j,m} ]; Pr["C is equal to"]; Pr[""]; PrW[" ", pat[c]//matform[]]; Pr[""]; Pr["This is the output returned"]; c ]//Catch; cofactorX[e___]:= (msg["argm",Length[{e}],1];fld); (* cramer *) (*Failure should return the input unaltered*) cramer[e___]:= With[{result = Block[{fn = cramer},cramerX[e] ]}, result/;result =!= fld ]; cramerX[A_,b_,opts___] := Module[{posnopt,expln,zero,zeroQ,m,n,Pr,PrW,DA,tA,x}, Which[ !MatrixQ[A]||!SameQ@@({m,n}= Dimensions[A]), msg["notsqrmat",1];thfld, !VectorQ[b]||Length[b] =!= m, If[m == 1, msg["rhs1", m], msg["rhs",m] ];thfld, (posnopt = Position[{opts}, a_/;!OptionQ[a],{1},1,Heads->False])=!= {}, msg["notopt", {opts}[[posnopt[[1,1]]]], 2];thfld ]; {m,n}= Dimensions[A]; (*Read options.*) {expln, zero} = {Explain, ZeroTest}/.{opts}/.Options[cramer]; zeroQ = TrueQ[zero[#]]&; (*Turn on explanation if required.*) If[expln, Pr = Print;PrW = printwide, Pr = PrW = Hold ]; top[x_] := Prepend[Table[{""},{m-1}], {x}]; (*Begin explanation*) Pr["We want to use Cramer©s rule to solve the equation"]; Pr[""]; PrW[" ", MapThread[Join,{A, top[". x = "], List/@b}]//matform[] ]; Pr[""]; Pr[" (call this A.x = b, with x = {x1, x2, ...})"]; Pr[""]; Pr["The Rule, which only works if det[A] is not 0, is that"]; Pr[""]; Pr[" xi = ", "det[Bi]"/"det[A]"]; Pr[""]; Pr["Where Bi is the result of replacing column i in A by b."]; Pr[""]; If[zeroQ[DA = Det[A]], Pr["For the given A, det[A] = 0, so the rule cannot be applied."]; msg["singmat"];thfld ]; Pr["For the given A and b it gives"]; Pr[""]; Pr[" x = "]; tA = Transpose[A]; x = Table[Det[ReplacePart[tA,b,i]], {i,n}]/DA; PrW[" ", x//matform[] ]; Pr[""]; Pr["This list is returned as the output."]; x ]//Catch; (*Trap and warn of any other faulty inputs*) cramerX[e___] := (msg["argm",Length[{e}], 2];fld) (* determ *) (*Failure must return the input unaltered*) determ[e___]:= With[{result = Block[{fn = determ},determX[e]] }, result/;result =!= fld ]; determX[A_, opts___] := Module[{posnopt,expln,det,Pr,PrW,pat,U,sign,}, Which[ !MatrixQ[A]||!SameQ@@Dimensions[A], msg["notsqrmat",1];thfld, (posnopt = Position[{opts}, a_/;!OptionQ[a],{1},1,Heads->False])=!= {}, msg["notopt", {opts}[[posnopt[[1,1]]]], 1];thfld ]; (*Read options.*) {expln, zero} = {Explain, ZeroTest}/.{opts}/.Options[determ]; tzQ = tozeroQ[zero, changed];(* see notes on tozeroQ*) changed = False; (*Guard against receiving unasked value.*) (*Turn on explanation if required.*) If[MatchQ[expln,True|PatternOnly], Pr = Print;PrW = printwide, Pr = PrW = Hold ]; If[expln === PatternOnly, pat = Map[If[MatchQ[#,0|0.|1|1.0],#,"+"]&,#,{2}]&, pat = Identity ]; Pr["We wish to find the determinant of the matrix"]; If[expln === PatternOnly, Pr["with pattern"]]; Pr[""]; PrW[" ", pat[A]//matform[] ]; Pr[""]; Pr["The U part of its PLU factorization is"]; {U,sign} = plu[A, ZeroTest -> zero][[{3,5}]]; Pr[""]; Pr[" ", pat[U]//matform[] ]; Pr[""]; Pr["The determinant required is the sign of the permutation"]; Pr["used in the PLU factorization times the product of the"]; Pr["diagonal entries of U"]; Pr["That is"]; Pr[""]; changed = False; (*Initialize changed*) det =sign Times@@Transpose[U,{1,1}] ; (*-- Transpose[M,{1,1}] happens to give the diagonal.*) PrW[ " ",det]; tzQ[det]; If[changed, Pr["Which simplifies to"]; Pr[""]; PrW[ " ",det] ]; Pr[""]; Pr["This is the output returned."]; det ]//Catch; (* Trap and warn of any other faulty inputs*) determX[e___] := (msg["argm",Length[{e}],1];fld); (* eigen *) eigen[x___]:= With[{ result = Block[{fn = eigen}, eigenX[x]] }, result/;result =!= fld ]; eigenX[A_,opts___] := Module[{expln, zero,zeroQ,Pr,PrW,tr,det,pol,vals,vecs}, Which[ !MatrixQ[A]||!SameQ@@Dimensions[A], msg["notsqrmat",1];thfld, (posnopt = Position[{opts}, a_/;!OptionQ[a],{1},1,Heads->False])=!= {}, msg["notopt", {opts}[[posnopt[[1,1]]]], 1];thfld ]; (*Read options.*) {expln, zero} = {Explain, ZeroTest}/.{opts}/.Options[eigen]; zeroQ = TrueQ[zero[#]]&; (*Turn on explanation if required.*) If[expln, Pr = Print;PrW = printwide, Pr = PrW = Hold ]; (*Begin explanation*) Pr["Call the given matrix A."]; Pr[""]; PrW["A = ", A//matform[]]; Pr[""]; Pr["The trace and determinant are:"] Pr[""]; tr = Plus@@Transpose[A, {1,1}]; Pr[" ", tr ]; det = Det[A]; Pr[" ", det ]; Pr[""]; Pr["Det[A - x I] is:"]; Pr[""]; pol = Det[ A - "x" IdentityMatrix[Length[A]]]; Pr[" ", pol]; Pr[""]; Pr["The list of eigenvalues is:"]; Pr[""]; {vals,vecs} = Eigensystem[A]; Pr[" ", vals]; Pr[""]; Pr["The list of eigenvectors is:"] Pr[""]; Pr[" ", vecs]; Pr[""]; Pr["The output is the list"]; Pr[""]; Pr["{trace, determinant, Det[A - x I], eigenvalues, eigenvectors}"]; {tr, det, pol, vals, vecs} ]//Catch; (* Trap and warn of any other faulty inputs*) eigenX[x___] := (msg["argm",Length[{x}],1];fld); (* eigen2 *) eigen2[x___]:= With[{ result = Block[{fn = eigen2}, eigen2X[x]] }, result/;result =!= fld ]; eigen2X[A_,opts___] := Module[ { posnopt,expln,zero,zeroQ,Pr,a11,a12,a21,a22,d,t, x1,x2,v1,v2}, Which[ !MatrixQ[A]||!Dimensions[A] == {2,2}, msg["not2by2"]; thfld, (posnopt = Position[{opts}, a_/;!OptionQ[a],{1},1,Heads->False])=!= {}, msg["notopt", {opts}[[posnopt[[1,1]]]], 1];thfld ]; (*Read options.*) {expln, zero} = {Explain, ZeroTest}/.{opts}/.Options[eigen2]; zeroQ = TrueQ[zero[#]]&; (*Turn on explanation if required.*) If[expln, Pr = Print, Pr = Hold]; {{a11,a12},{a21,a22}} = A; d = a11 a22 - a12 a21; t = a11 + a22; {x1,x2} = Cancel[(t + Sqrt[t^2-4d]{1,-1})/2]; {v1,v2} = Which[ !zeroQ[a12], {{a12,x1-a11},{a12, x2-a11}}, !zeroQ[a21], {{x1 -a22, a21}, {x2 -a22, a21}}, True, {{1,0},{0,1}} ]//Cancel; Pr["Let the given matrix be called A:"]; Pr[""]; Pr["A = "]; Pr[" ", A//MatrixForm]; Pr[""]; Pr["The determinant of A - x I is:"]; Pr[""]; Pr[ " ","x"^2," - ", t"x", " + ", d, " = ", 0]; Pr[""]; Pr["The first eigenvalue and eigenvector are {x1, v1} = "] Pr[""]; Pr[" ", {x1,v1}]; Pr[""]; Pr["The second eigenvalue and eigenvector are {x1, v1} = "] Pr[""]; Pr[" ", {x2,v2}]; Pr[""]; Pr["The output is {{x1,x2},{v1,v2}}"]; {{x1,x2},{v1,v2}} ]//Catch; (* Trap and warn of any other faulty inputs*) eigen2X[e___] := (msg["argm",Length[{e}],1];fld); (* findpiv *) findpiv[e___]:= With[{result = Block[{fn = findpiv}, findpivX[e]]}, result/; result =!=fld ]; findpivX[A_, rr_, cc_, opts___]:= Module[{posnopt,m,n,r,s,c,d,expln,zero,tzQ,U,tab,Pr,PrW,pat,pos}, Which[ (*input checks and error messages*) !MatrixQ[A], msg["notmat",1];thfld, {m,n} = Dimensions[A]; !( VectorQ[rr,IntegerQ] && Length[rr] == 2 && 1<=rr[[1]]<=rr[[2]]<=m || IntegerQ[rr] && 1<=rr<=m ), msg["notnl", 2,1,m,rr];thfld, !(VectorQ[cc,IntegerQ] && Length[cc] == 2 && Min[n,1]<=cc[[1]]<=cc[[2]]<=n || IntegerQ[cc] && Min[n,1]<=cc<=n ), If[n>0,msg["notnl", 3,1,n,cc], msg["notzero",3,cc] ];thfld, !((posnopt = Position[{opts}, a_/;!OptionQ[a],{1},1,Heads->False])=== {}), msg["notopt", {opts}[[posnopt[[1,1]]]], 3];thfld ]; (*Normalise inputs if necessary; use local names.*) {r,s} = If[IntegerQ[rr],{rr,m},rr]; {c,d} = If[IntegerQ[cc],{cc,n},cc];; (*Read options.*) {expln, zero} = {Explain,ZeroTest}/.{opts}/.Options[findpiv]; tzQ = tozeroQ[zero, changed];(* see notes on tozeroQ*) changed = False; (*Guard against receiving of false value True*) (*Give a local name to A.*) U = A; (*Form the table of positions at which to seek a pivot*) tab = Table[{i,j}, {i,r,s},{j,Max[c,1],d}]; (*Turn on explanation if required.*) If[MatchQ[expln,True|PatternOnly], Pr = Print;PrW = printwide, Pr = PrW = Hold; ]; If[expln === PatternOnly, pat = Map[If[MatchQ[#,0|0.|1|1.0],#,"+"]&,#,{2}]&, pat = Identity ]; (*State inputand explain the method.*) If[expln === PatternOnly, Pr["The pattern of the given matrix is"], Pr["The given matrix is"] ]; Pr[""]; PrW[" ",pat[A]//matform[TableSpacing -> {0,3}]]; Pr[""]; If[{r,s} == {1,m} && {c,d} == {Min[1,n],n}, Pr["Search the matrix column by column for a non-zero entry."], Pr[StringForm["Search the submatrix A[``,``]",{r,s},{c,d}]]; Pr[""]; PrW[" ", MapAt[StringForm["<``>",#]&, pat[U], Join@@tab]//matform[] ]; Pr[""]; Pr["column by column for a non-zero entry."] ]; Pr[""]; (*Find if there is a pivot: in the course of the search replace with zero any non zero entries that test as zero.*) changed = False; (*Set indicator for having replaced by 0*) pos = Cases[ Transpose[tab],(*Transpose so that search is by rows*) {i_,j_}/;!tzQ[U[[i,j]]], {2}, 1 ]//Flatten; (*Describe what has been found*) If[pos =={}, Pr["No non-zero entry has been found."]; If[ changed, Pr[""]; Pr["The matrix, A, has simplified to"]; Pr[""]; PrW[" ", pat[U]//matform[] ]; Pr[""]; Pr["There no pivot."] ], (*else*) Pr[StringForm[ "The first non-zero entry``is at ``: ", If[ changed, " (after some simplification of A) ", " " ], pos ]]; Pr[""]; PrW[" ",MapAt[StringForm["<``>",#]&, pat[U], pos]//matform[]]; Pr[""]; Pr["So the required pivot position is ", pos,"."]; ]; (*Describe the output*) Pr[""]; If[changed, If[pos == {}, Pr["The output is the list {{}, simplified A }."], Pr["The output is the list { pivot position, simplified A }."] ], If[pos == {}, Pr["The output is the list { {}, A }."], Pr["The output is the list { pivot position, A }."] ] ]; {pos,U} ]//Catch; (*trap and warn of any other faulty inputs*) findpivX[e___] := (With[{l = Length[{e}]},Switch[l,1,msg["argmu",3],_,msg["argm",l,3]]]; fld); (* grams *) (*Failure should return the input unaltered*)grams[x___]:= With[{ result = Block[{fn = grams}, gramsX[x]]}, result/;result =!= fld ]; (*Although we are dealing with columns it is more convenient to work with rows in Mathematica.*) gramsX[A_, opts___]:= Module[ {posnopt,m,n,L,zero,expln,ip,zeroQ,zerovecQ,zerovec,d2,Pr,PrW,p, sp,rec,proj,subproj,ssubproj,ssubprojval,trA,na,cmp,nmlzdp,Q,R, pos,zeropos }, Which[ (*Input checks and error messages.*) !MatrixQ[A], msg["notmat",1];thfld, (posnopt = Position[{opts}, a_/;!OptionQ[a],{1},1,Heads->False])=!= {}, msg["notopt", {opts}[[posnopt[[1,1]]]], 1];thfld ]; {m,n} = Dimensions[A]; (*Read the options settings*) {zero, nmlzQ, expln, ip} = {ZeroTest, Normalized, Explain,InnerProduct}/.{opts}/.Options[grams]; zeroQ = TrueQ[zero[#]]&; zerovecQ[v_] := And@@(zeroQ/@v); zerovec = Table[0,{m}];(*Standard zero vector*) If[ip == Automatic, If[FreeQ[A,Complex], ip = (Transpose[#1].#2)&, ip = (Conjugate[Transpose[#1]].#2)& ] ]; d2[x_] := ip[x,x];(*Square of norm function.*) (*Begin explanation if requested*) If[expln, Pr = Print;PrW = printwide, Pr = PrW = Hold ]; Pr["The inner product being used is"]; Pr[""]; PrW[" x´y = ", matform[]//@Block[{Dot},ip["x","y"]]]; Pr[""]; Pr["The given matrix, A, is"]; Pr[""]; PrW[" ", A//matform[]]; Pr[""]; Pr["Call its columns a1 to ","a", n,"."]; Pr[""]; If[nmlzQ, Pr["We construct an orthonormal matrix Q and an upper triangular"]; Pr["matrix R such that A = QR."], Pr["We construct an orthogonal matrix Q and an upper triangular"]; Pr["matrix R such that A = QR."] ]; Pr[""]; (*Remove any previous instances of the working symbols - needed because we use Unique to generate them -Clearing is not sufficient.*) removeUnique["a","p"]; (*Initialize the sequences of working symbols and strings.*) p = sp = rec = {}; Pr["First find Q."]; Pr[""]; If[nmlzQ, Pr["Begin by constructing an orthogonal sequence p1, p2, ..."]; Pr["as follows:"], Pr["The columns p1, p2, ... of Q are constructed in sequence"]; Pr["as follows:"]; ]; Pr[" (proj[a,b] is the projection of a onto b, ", "b´a"/"b´b"," b )."]; (*Project a onto b*) proj[a_,b_] := ip[b,a]/ip[b,b] b; (*Subtract from a its projections onto the elements of the list q - this is done step by step*) subproj[a_,q_] := Fold[(#1 - proj[#1,#2])&, a, q]; (*Two explanatory forms of the preceding*) ssubproj[x_,y_] := Sequence[x, Sequence@@Map[ StringForm[" - proj[``,``]",x, #]&, y ]]; ssubprojval[x_,y_] := Sequence[x, Sequence@@Map[ StringForm[" - ``",proj[x,#]]&, y ]]; (*Construct p1,p2,...*) (*Switch to working with rows*) trA = Transpose[A]; Scan[(*Scan trA*) ( na = Unique["a"];sna = ToString[na]; na = #; cmp = subproj[na,p]; Pr[" "]; If[sp==={}, Pr[" ",ssubproj[sna,sp]," = ", cmp], PrW[" ",ssubproj[sna,sp]]; PrW[" "," = ", ssubproj[na,p]]; PrW[" "," = ", ssubprojval[na,p]]; PrW[" "," = ", cmp] ]; If[ zerovecQ[cmp], AppendTo[rec,0]; If[cmp != zerovec, cmp = zerovec PrW[" "," = ", cmp] ]; Pr[" ","Ignore this."], np = Unique["p"];snp = ToString[np]; np = cmp; AppendTo[p, np];AppendTo[sp, snp];AppendTo[rec,1]; PrW[" ",snp, " = ", np] ]; If[Length[sp] == m < n, Pr[""]; If[m =!= 1, Pr["We now have ",m," independent vectors in a ",m]; Pr["dimensional space, so no further additions can"]; Pr["be made to the sequence."], Pr["We now have one independent vector in an one"]; Pr["dimensional space, so no further additions can"]; Pr["be made to the sequence."] ]; Return[Null] ] )&, trA ]; If[p == {}, Pr[""]; Pr["We have an empty sequence. There seems to be no solution of the kind"]; Pr["that we want. The best we can do is to define Q to be the m by 0 matrix."]; Pr["It is orthonormal - it has no columns to be otherwise with - and it is an"]; Pr["orthonormal basis for the column space of A"]; Pr["The output is {Q, $Failed, ip}"]; msg["nocols",m]; Throw[{Table[{},{m}], fld, ip}] ]; (*normalize p1, p2, ... if requested*) If[nmlzQ, Pr[""]; Pr["Normalize the sequence to find the columns of Q."]; Pr[""]; nmlzdp = p/Sqrt[d2/@p]; If[ TrueQ[p == nmlzdp], Pr["The sequence is already orthonormal; this gives the columns of Q.."], p = nmlzdp ] ]; (*Define Q -- back to column vectors*) Pr[""]; Pr["Q is "]; Pr[""]; Q = Transpose[p]; PrW[" ",Q//matform[]]; (*Find R*) Pr[""]; Pr["Now find R:"]; Pr[""]; Which[ nmlzQ, Pr[" Q´A = Q´(QR) =, (Q´Q)R = IR = R"]; Pr[""]; Pr["So R is Q´A ="]; R = ip[Q,A], !nmlzQ, Pr[" Q´A = Q´(QR) = (Q´Q)R = DR"]; Pr["where D is the diagonal matrix whose diagonal entries are the"]; Pr["the squares of the lengths of the columns of Q (these are all"]; Pr["non-zero)."]; Pr[""]; Pr["So R = ","D"^"-1","(Q´A) ="]; Pr[""]; R = Inverse[ip[Q,Q]].ip[Q,A]; ]; (*Make the entries that must be zero equal to zero: the calculation may show them as not zero because of numerical approximation or lack of symbolic simplification.*) (*Construct the list of the positions of such entries R.*) pos = Flatten[Position[rec,1]]; zeropos = Cases[ Flatten[Array[{#1,#2}&,{Length[pos],n}],1], {i_,j_}/;jFalse])=!= {}, msg["notopt", {opts}[[posnopt[[1,1]]]], 1];thfld ]; (*Read options.*) {expln, zero} = {Explain, ZeroTest}/.{opts}/.Options[invert]; zeroQ = TrueQ[zero[#]]&; {m,n}=Dimensions[A]; (*If Explain is set to True, turn on the printing.*) If[expln, Pr = Print;PrW = printwide, Pr = PrW = Hold ]; Pr["We want to find the inverse of the matrix"]; Pr[""]; PrW[" ", A//matform[]]; Pr[""]; Pr["First form AI = [A I] where I is the m by m identity matrix."]; ai = Transpose[Join[Transpose[A],IdentityMatrix[m]]]; Pr[""]; PrW[" ", ai//matform[]]; Pr[""]; Pr["Find the reduced echelon form RAI of AI by the Gauss-Jordan"]; Pr[" method (using the function ref)."]; refai = ref[ai, ZeroTest -> zero][[1]]; (* -- pass zero to ref *) Pr[""]; Pr[" ", refai//matform[]]; Pr[""]; If[zeroQ[refai[[m,n]]], Pr["The matrix A cannot not have an inverse; for suppose it does have"] Pr["an inverse:"]; Pr["Let G be the matrix giving the row reduction used to get RAI."] Pr["Then"]; Pr[""]; Pr[" ", "G.[A I] = [G.A G.I]"]; Pr["On identifying G.A in RAI we see that"]; Pr["G.A = "]; Pr[""]; PrW[" ", Take[#,n]&/@refai//matform[] ]; Pr[""]; Pr["And"]; Pr[" ", " I = (G.A).","A"^"-1","G"^"-1"]; Pr[""]; Pr["Which is impossible since the right hand side has a zero last row"]; Pr[""]; msg["singmat",1];thfld ]; Pr["If G is the matrix giving the row reduction used then"]; Pr[""]; Pr[" ", "G.[A I] = [G.A G.I] = [I G]."]; Pr[""]; Pr["So, G.A = I and G is the inverse of A"]; Pr["On identifying G in RAI we get"]; G = Drop[#,n]&/@refai; Pr[""]; Pr["A"^"-1"," = G = "]; Pr[""]; PrW[" ", G//matform[]]; Pr[""]; Pr["This is the output returned."]; G ]//Catch; (* Trap and warn of any other faulty inputs*) invertX[e___] := (msg["argm",Length[{e}],1];fld); (* linefit *) linefit[x___]:= With[{ result = Block[{fn = linefit}, linefitX[x]]}, result/;result =!= fld ]; linefitX[ xs_, vals_, x_, opts___] := Module[{lx,lv,Pr,PrW,points,fmla,min,max,delta,minv,maxv,deltav}, Which[ !VectorQ[xs, (NumberQ[#]&&FreeQ[#,Complex])&]||(lx=Length[xs])<=0, msg["notnzreallst",1];thfld, !VectorQ[vals, (NumberQ[#]&&FreeQ[#,Complex])&]||(lv=Length[vals])<=0, msg["notnzreallst",2];thfld, lx=!=lv, msg["badlnth",2,lx,1,lv];thfld, Head[x]=!=Symbol, msg["notsym",3,x];thfld, (posnopt = Position[{opts}, a_/;!OptionQ[a],{1},1,Heads->False])=!= {}, msg["notopt", {opts}[[posnopt[[1,1]]]], 3];thfld ]; (*Turn on explanation if required.*) If[Explain/.{opts}, Pr = Print;PrW = printwide, Pr = PrW = Hold ]; Pr["The built-in function Fit is used to find the formula"]; Pr["for the best fit line. It gets"]; points = Transpose[{xs,vals}]; fmla = Fit[points, {1,x},x]; Pr[""]; Pr[" ", fmla]; Pr[""]; Pr["The function Plot is then used to graph this together with"]; Pr["the points given by the given data; these points are added by"]; Pr["using the Epilog option for Plot"]; Pr["You might like to look at the function Fit - it does a lot"]; Pr["more than we have used it for here. Also, you can use all the"]; Pr["various options for Plot to change picture - but beware that"]; Pr["if you use the option Epilog you will remove the data points."]; Pr[""]; Pr["The output returned is the formula found above."]; {min, max} = {Min[xs], Max[xs]}; delta = 0.1(max - min); {minv, maxv} = {Min[vals], Max[vals]}; deltav = 0.1(maxv - minv); plot = Plot[ fmla, {x, min-delta, max+delta}, Evaluate[FilterOptions[Plot,opts]], PlotStyle -> RGBColor[0,0,1], PlotRange -> {minv - deltav, maxv + deltav}, Epilog -> {PointSize[.02],RGBColor[1,0,0],Point/@points}, PlotLabel -> SequenceForm[ "y = ",N[fmla,2]/.x->"x"], AxesLabel -> {"x", "y"} ]; fmla ]//Catch; (*Trap and warn of any other faulty inputs*) linefitX[e___] := (msg["argm",Length[{e}],3];fld); (* lsq *) (* In case of faulty input or failure of computation return input unaltered and display messages. *) lsq[x___]:= With[{result = Block[{fn = lsq}, lsqX[x]]}, result/; result =!= $Failed ]; lsqX[A_,b_, opts___] := Module[{posnopt,m,n,expln,Sh,Pr,PrW,x,p,e}, Which[ !MatrixQ[A], msg["notmat",1];thfld, {m,n}=Dimensions[A]; !VectorQ[b]||Length[b] =!= m, If[m == 1, msg["rhs1", m], msg["rhs",m] ];thfld, (posnopt = Position[{opts}, a_/;!OptionQ[a],{1},1,Heads->False])=!= {}, msg["notopt", {opts}[[posnopt[[1,1]]]], 2];thfld ]; (*Read options.*) expln = Explain/.{opts}/.Options[lsq]; (*Turn on explanation if required.*) If[expln, Sh = Show; Pr = Print;PrW = printwide, Pr = PrW = Sh = Hold ]; top[x_] := Prepend[Table[{""},{m-1}], {x}]; (*Begin explanation: deal with the special case n = 0.*) If[n>0, Sh[fig]]; Pr["Often an equation, A.x = b, will not have an exact solution."]; Pr["We can then look for an x that makes the length of the"]; Pr["vector A.x - b as small as possible."]; If[n==0, Pr["In this case A has no columns, so the only solution for x is"]; Pr["the empty vector {}."]; Pr["The output is the triple {x, A.x, A.x - b}, that is"]; Pr[""]; Pr[" ", {{},0,-b}]; Throw[{{},0,-b}] ]; Pr["Geometrically, we need an x such that A.x is the projection"]; Pr["of b onto the linear space spanned by the columns of A, that"]; Pr["is, A x - b is orthogonal to each column in A."]; Pr["Algebraicly, this means that x is the solution of"]; Pr[""]; Pr[" ", "A"^"T",".(A.x - b) = 0" ]; Pr[""]; Pr["which is equivalent to"]; Pr[""]; Pr[" ", "A"^"T",".A.x = ", "A"^"T",".b", ", the \"normal equation\". "]; Pr[""]; Pr["Any such solution is called a \"least squares solution\""]; Pr["of the original equation because it minimizes the sum of"]; Pr["the squares of the entries in A x - b."]; Pr[""]; Pr["The given equation is"]; Pr[""]; PrW[" ", MapThread[Join,{A, top[" x ="], List/@b}]//matform[] ]; Pr[""]; Pr["Find a particular solution using solve"]; Pr[""]; Pr[" x = solve[Transpose[A].A, Transpose[A].b]"]; x = solve[Transpose[A].A, Transpose[A].b]; Pr[""]; PrW[" = ", x]; Pr[""]; Pr["Then"]; p = A.x; PrW[" p = A.x = ",p]; Pr[""]; Pr["So the error, the vector from b to p, is"]; Pr[""]; e = p -b; PrW[" e = p - b = ", e]; Pr[""]; PrW["Its length is ", N[Sqrt[e.e]]]; Pr[""]; Pr["The output returned is the triple {x, p, e}"]; {x,p,e} ]//Catch; (* Trap and warn of any other faulty inputs*) lsqX[e___] := (msg["argm",Length[{e}],2];fld); (*Code for the figure used in the explanation. The form fig := fig = ... saves the output of fig the first time fig is called and avoids the need for recomputation on subsequent calls*) fig := fig = With[ { z = {0,0}, b = {2,3}, a1 = {1,-1.8},a2 = {4,.5}, p = {2,-.7} }, Graphics[ { Pink, Polygon[{{0,0},a1,a1+a2,a2}], Blue,Arrow[z,b], Black,Text["b",b,{-2,-1}], Red,Arrow[z,a1],Arrow[z,a2], Black, Text["a1",a1,{0,3}], Text["a2",a2,{-2,-1}], Purple,Arrow[z,p], Black,Text["p = A.x",p, {-1.5,0}], Green,Arrow[b,p], Black, Text["e", .5(p+b),{3,0}], Text[ ColumnForm[ {SequenceForm["a",ColumnForm[{"T","",1}, Center,Center],".e = 0"], "\n\n\n", SequenceForm["a",ColumnForm[{"T","",2},Center, Center],".e = 0"], "\n\n\n", SequenceForm["A"^"T",".e = ","A"^"T",".(b - A x) = 0"] }, Left, Center ], (b+p)/2, {-1.5,-1} ] }, AspectRatio -> Automatic, PlotRange -> All ] ]; (* null *) null[x___]:= With[{ result = Block[{fn = null}, nullX[x]]}, result/;result =!= fld ]; nullX[A_, opts___] := Block[{posnopt,expln,zero,zeroQ,Pr,PrW,R,pivcol,k,l,vars,K,npR,Id}, (*input checks and error messages*) Which[ !MatrixQ[A], msg["notmat",1];thfld, !((posnopt = Position[{opts}, a_/;!OptionQ[a],{1},1,Heads->False])=== {}), msg["notopt", {opts}[[posnopt[[1,1]]]], 1, fn];thfld ]; {m,n} = Dimensions[A]; (*Read options.*) {expln, zero} = {Explain, ZeroTest}/.{opts}/.Options[null]; zeroQ = TrueQ[zero[#]]&; (*Turn on explanation if required.*) If[expln, Pr = Print;PrW = printwide, Pr = PrW = Hold ]; Pr["Given the matrix"]; Pr[""]; PrW[" ", A//matform[]]; Pr[""]; Pr[" (Call this A)"] Pr["We want to find an matrix N which has a basis for"] Pr["the null space of A as its columns."]; Pr["Find the reduced row echelon form R of A using the function"]; Pr["ref."]; Pr["Call its pivot variables p1,¼, and its free variables x1,¼ :"]; {R,pivcol} = ref[A, ZeroTest -> zero]; nonpiv = Complement[Range[n],pivcol]; p = Length[pivcol]; Pr[""]; Pr[" R = "]; Pr[""]; PrW[" ", R//matform[ TableHeadings -> {None, vars = ( k=1; MapAt[ ToString[SequenceForm["x",k++,"\n"]]&, l=1; MapAt[ ToString[SequenceForm["p",l++,"\n"]]&, Table["",{n}], List/@pivcol ], List/@nonpiv ] ) } ] ]; Pr[""]; PrW["The variables vector is v = ", vars]; Pr["We are seeking NON ZERO solutions to R.v = 0"]; Which[ p==n, Pr["In this case there are no such solutions."], 0 == p==n-1, Pr["Here there is, up to a scalar multiple, only one."]; Pr["It is x1 == 1."], p==n-1 , Pr["Here there is, up to a scalar multiple, only one."]; Pr["It is found by giving the value 1 to x1 and solving"]; Pr["for the pivot variables."], True, Pr["The first basis vector needed is found by giving the value 1"]; Pr["to x1, 0 to the other xi©s and then solving for the pivot"]; Pr["variables;and so on."] ]; npR = R[[Range[m],nonpiv]]; K = Array[0&,{n,n-p}]; Id = IdentityMatrix[n-p]; Do[K[[ pivcol[[i]] ]] = - npR[[i]], {i,p}]; Do[K[[ nonpiv[[i]] ]] = Id[[i]], {i,n-p}]; Pr[""]; Pr["The required matrix is"]; Pr[""]; PrW[" ", K//matform[]]; Pr[""]; Pr["This is the output returned."]; K ]//Catch; (*Trap and warn of any remaining kind of faulty inputs.*) nullX[e___] := (msg["argm", Length[{e}],1];fld); (* permdet *) permdet[x___]:= With[{ result = Block[{fn = permdet}, permdetX[x]]}, result/;result =!= fld ]; permdetX[P_, opts___] := Module[{posnopt,expln,Pr,PrW,mat,p,pi,pos}, Which[ !VectorQ[P,IntegerQ[#]&&Positive[#]&], msg["notposintlst",1];thfld, (posnopt = Position[{opts},a_/;!OptionQ[a],{1},1,Heads->False])=!= {}, msg["notopt", {opts}[[posnopt[[1,1]]]], 2];thfld ]; (*Read options.*) {expln} = {Explain}/.{opts}/.Options[permdet]; (*Turn on explanation if required.*) If[expln, Pr = Print;PrW = printwide, Pr = PrW = Hold ]; (*Begin explanation*) Pr["The given list is a list of positive integers."]; Pr[""]; Pr[" ", P]; Pr[""]; ln = Length[P]; Pr["of length ",ln,"."]; Pr["Call it P with entries p1, p2, ... ."]; Pr["Its permutation matrix is the ",ln," by ",ln, " matrix M in"]; Pr["which row i has a 1 in position pi and is zero everywhere else."]; Pr["It has the property that M.", Switch[ln, 0, {}, 1, {1}, 2, {1,2}, 3, {1,2,3}, _, StringForm["{1,2,...,``}", ln] ], " = P." ]; id = IdentityMatrix[Length[P]]; M = id[[ P[[#]] ]]&/@ Range[Length[P]]; Pr["It is"]; Pr[""]; PrW[" ", M//matform[]]; Pr[""]; Pr["The determinant of P is defined to be the determinant of M."]; det = Det[M]; Pr["Its value is ", det]; If[!MatchQ[det,-1|1], Pr["Since the values is not -1 or 1 P is not a permutation."] msg[notpermM, 1, det];thfld ]; Pr["This is the output returned."]; det ]//Catch; (* Trap and warn of any other faulty inputs*) permdetX[e___] := (msg["argm",Length[{e}],1];fld); (* plu *) (* in case of faulty input or failure of computation return input unaltered and display messages *) plu[x___]:= With[{ result = Block[{fn = plu}, pluX[x]]}, result/;result =!= fld ]; pluX[A_,opts___]:= Module[ {posnopt,m,n,P,L,U,pivcol,sign,expln,zero,tzQ,Pr,PrW,pat, tf3,pos,r,c,rowredQ}, (*input checks and error messages*) Which[ !MatrixQ[A], msg["notmat",1];thfld, !((posnopt = Position[{opts}, a_/;!OptionQ[a],{1},1,Heads->False])=== {}), msg["notopt", {opts}[[posnopt[[1,1]]]], 1, fn];thfld ]; {m,n} = Dimensions[A]; (* initialize quantities to be computed*) P = L = IdentityMatrix[m]; U = A; pivcol = {}; sign = 1; (*Read options*) {expln, zero} = {Explain, ZeroTest}/.{opts}/.Options[plu]; tzQ = tozeroQ[zero, changed];(* see notes on tozeroQ*) changed = False; (*Guard against receiving the value True*) (*Turn on explanation if required.*) If[MatchQ[expln,True|PatternOnly], Pr = Print;PrW = printwide; tf3[x___] := TableForm[{x}, TableDirections->{Row,Column,Row}, TableSpacing ->{8,0,3}, TableHeadings -> {SequenceForm[ SpaceForm[Floor[(3m-1)/2]],#]&/@{"P\n","L","U"},None,None}, TableAlignments -> If[TrueQ[MatrixQ[MM,IntegerQ]], {Center,Center,Right}, {Center,Center,Center} ] ], Pr = PrW = tf3 = Hold ]; If[expln === PatternOnly, pat = Map[If[MatchQ[#,0|0.|1|1.0],#,"+"]&,#,{2}]&, pat = Identity ]; (*State input,and kind of commentary to be given -- explanations to be printed all begin with Pr and are fairly self-explanatory.*) If[expln === PatternOnly, Pr["The pattern of the given matrix is"], Pr["The given matrix is"] ]; Pr[""]; PrW[" ",pat[U]//matform[TableSpacing -> {0,3}]]; Pr[""]; Pr["Call it A."]; Pr["We construct matrices P, L, U such that P.A = L.U and"]; Pr[" P is a permutation matrix;"]; Pr[" L is square lower triangular with ones on the diagonal;"]; Pr[" U is upper triangular."]; Pr[""]; Pr["The construction starts with the triple below: the \ncommentary will be concerned with U; consequential changes \nin P and L are shown without comment"]; Pr[""]; PrW[tf3[P,L,pat[U]]]; Pr[""]; (*For k from 1 to m look for a pivot in row k*) Do[ changed = False; (*changed will be set to True if tzQ replaces any nonzero entry with zero*) {pos,U} = findpiv[U, k, Min[Max[Min[1,n],pivcol+1],n], ZeroTest ->zero]; (*-- Pass zero into findpiv, and pass any changes to the U inside findpiv out to the U on the left side.*) If[ pos =!= {}, (*if there is a pivot*) (*Name the row and column indices.*) {r,c} = pos; (*Append the column index to pivcol*) AppendTo[pivcol, c]; Pr[""]; Pr["Pivot ",k," in U is at ",{r,c},"."]; Pr[""]; PrW[tf3[ P, pat[L], MapAt[StringForm["<``>",#]&,pat[U],{r,c}] ]]; Pr[""]; (*If necessary swap rows to put this pivot in row k.*) If[r > k, (*Swap rows r and k in U*) {U[[k]],U[[r]]} = {U[[r]],U[[k]]}; (*make the consequential adjustments to P and L*) (*Swap rows r and k in P*) {P[[k]],P[[r]]} = {P[[r]],P[[k]]}; (*Swap rows r and k in L - but only up to position k-1.*) Do[ {L[[k,i]],L[[r,i]]} = {L[[r,i]],L[[k,i]]}, {i,k-1} ]; (*change sign to record the swap*) sign = -sign; Pr["Swap rows ",k," and ",r," to change it to ",{k,c},"."]; Pr[""]; PrW[tf3[ P, pat[L], MapAt[StringForm["<``>",#]&,pat[U],{k,c}] ]]; Pr[""]; Pr["Pivot ",k," is now at ",{k,c},"."]; ]; (* Check the entries below the pivot. Assign value zero to any that simplify to zero. Row reduce the others to zero *) rowredQ = False; (*-- rowredQ will be set to True if there are any row reductions *) Do[ If[!tzQ[U[[i,c]]], rowredQ = True; prech = changed;(*note the values of changed*) With[{mult = U[[i,c]]/U[[k,c]]}, (* Subtract mult = U[[i,c]]/U[[k,c]] times row k from row i. *) U[[i]] = U[[i]] - mult U[[k]]; tzQ[U[[i,c]]]; (*to deal with rounding problems*) (*Set L[[i,k]] to mult: the effect of adding mult times column i to column k. *) L[[i,k]] = mult ]; changed = prech(*restore the values of changed*) ], {i,k+1,m} ]; (*Explain, if k < m.*) If[ k",#]&,pat[U],{k,c}] ]]; Pr[""], (*else*) If[changed, Pr["After some simplification, the entries below the", "pivot "]; Pr[{k,c}," in U are all zero:"]; Pr[""]; PrW[tf3[ P, pat[L], MapAt[StringForm["<``>",#]&,pat[U],{k,c}] (* MapAt[StringForm["[``]",#]&,pat[U],{#,c}&/@Range[k+1,m]] *) ]]; Pr[""], Pr["The entries below the pivot ",{k,c}, " in U are already zero." ]; ] ] ], (*Else*) (*Exit from If and Do - there is no point in looking at any remaining rows.*) Throw[""] ], {k,1,m } ]//Catch; (*Explain that the calculation is complete.*) Pr[""]; If[changed, Pr["After some simplification we get."]; Pr[""]; PrW[tf3[P,pat[L],pat[U]]]; Pr[""] ]; Pr["The matrix U is upper triangular."]; Pr["The calculation is complete."]; Pr[""]; Pr["The output is the list {P,L,U,pivcol,sign} where:"]; Pr[" - pivcol is the list pivot column positions;"]; Pr[" - sign is -1^(number of swaps used)."]; {P,L,U,pivcol,sign} (*The output.*) ]//Catch; (*Trap and warn of any other faulty inputs.*) pluX[e___] := (msg["argm", Length[{e}],1];fld); (* projmat *) (*in case of faulty input or failure return input unaltered and with error messages *) projmat[x___]:= With[{result = Block[ {fn = projmat}, projmatX[x]]}, result/; result =!= fld ]; projmatX[A_,opts___] := Module[{posnopt,m,n,expln,zero,zeroQ,Pr,PrW,G,ans,B}, (*input checks and error messages*) Which[ !MatrixQ[A], msg["notmat",1]; thfld, (posnopt = Position[{opts}, a_/;!OptionQ[a],{1},1,Heads->False])=!= {}, msg["notopt", {opts}[[posnopt[[1,1]]]], 2];thfld ]; {m,n} = Dimensions[A]; (*Read options.*) {expln, zero} = {Explain, ZeroTest}/.{opts}/.Options[projmat]; zeroQ[x_] := TrueQ[zero[x]]; (*Turn on the explanation if required.*) If[expln, Pr = Print;PrW = printwide, Pr = PrW = Hold ]; (*Begin explanation*) Pr["Given the matrix "]; Pr[""]; PrW[" ", A//matform[]]; Pr[""]; Pr[" (call this A)."]; Pr[""]; Pr["We construct a matrix P such that the projection of any"]; Pr["vector b onto the column space of A is P.b"]; Pr["(P is necessarily square)."]; If[n==0, Pr["In this case A has no columns, so the column space of A"]; Pr["consists of just the zero vector, and P must be the"]; Pr[m, " by ", m, " zero matrix."]; Pr["This is the output returned."]; Throw[ Array[0&,{m,m}] ] ]; Pr[""]; Pr["Since P.b is in the column space of A"]; Pr[""]; Pr[" P.b = A.x for some x."]; Pr[""]; Pr["Since A.x is the projection of b onto the column space of A"]; Pr[""]; Pr[" A x - b is orthogonal to each column in A."]; Pr[""]; Pr["Algebraicly, this means that x is the solution of"]; Pr[""]; Pr[" ", "A"^"T",".(A.x - b) = 0" ]; Pr[""]; Pr["which is equivalent to"]; Pr[""]; Pr[" ", "A"^"T",".A.x = ", "A"^"T",".b", ", the \"normal equation\". " ]; Pr[""]; Pr["If ", "A"^"T",".A is invertible then"]; Pr[""]; Pr[" x = ","(A"^"T",".A)"^("-1"),".A"^"T",".b"]; Pr[""]; Pr["and"]; Pr[""]; Pr[" P.b = A.x"]; Pr[" = A.(A"^"T",".A)"^("-1"),".A"^"T",".b"]; Pr[""]; Pr["So,"]; Pr[" P = A.(A"^"T",".A)"^("-1"),".A"^"T"]; Pr[""]; Pr["is the unique solution for P."]; Pr[""]; Pr[ "If ", "A"^"T", ".A is not invertible then we can begin by replacing it"]; Pr["by a matrix B which has a basis for the column space of A as"]; Pr["its columns."]; Pr["Then"]; Pr[" P = B.(B"^"T",".B)"^("-1"),".B"^"T"]; Pr[""]; Pr[" ****************** "]; Pr[""]; G = Transpose[A].A; If[!zeroQ[Det[G]], Pr["For the given A the determinant of (A"^"T",".A) is not zero."]; Pr["So A"^"T",".A is invertible and"]; Pr[""]; Pr[" P = A.(A"^"T",".A)"^("-1"),".A"^"T"]; Pr[""]; ans = A.Inverse[G].Transpose[A]; PrW[" = ", ans//matform[]], (*else*) Pr["For the given A the determinant of (A"^"T",".A is zero;"]; Pr["so A"^"T",".A is not invertible, and we find B = "]; Pr[""]; B = basis[A, ZeroTest -> zero]; Pr[" ", B//matform[] ]; Pr[""]; Pr["And then"]; Pr[""]; Pr[" P = B.(B"^"T",".B)"^("-1"),".B"^"T"]; Pr[""]; ans = B.Inverse[Transpose[B].B,ZeroTest -> zero].Transpose[B]; PrW[" = ", ans//matform[]]; ]; Pr[""]; Pr["This is the output returned"]; ans ]//Catch; (* trap and warn of any other bad inputs*) projmatX[e___] := (msg["argrx",Length[{e}],1];fld); (* randperm *) randperm[x___]:= With[{ result = Block[{fn = randperm}, randpermX[x]]}, result/;result =!= fld ]; (* I use Block instead of Module so that the local variable will print without an appended module number*) randpermX[n_,opts___] := Block[{posnopt,expln,Pr,t,s, ans}, Which[ !IntegerQ[n]||n<0, msg["notnnint", 1, n];thfld, (posnopt = Position[{opts},a_/;!OptionQ[a],{1},1,Heads->False])=!= {}, msg["notopt", {opts}[[posnopt[[1,1]]]], 2];thfld ]; {expln} = {Explain}/.{opts}/.Options[randperm]; (*Turn on explanation if required.*) If[expln, Pr = Print]; (* Pr is not set to Hold if expln is not True:this allows evaluation of te terms in expressions ofr the form Pr[...], as is needed below*) Pr["We construct a random permutation of the integers 1 to ",n,"."]; Pr["First make a table of ",n," pairs, of the form"]; Pr[""]; Pr[" {random real, i}, i = 1 to ",n, ":"]; Pr[""]; Pr[" t = Table[{Random[],i}, {i,",n,"}] ="]; Pr[""]; Pr[" ",(t = Table[{Random[],i}, {i,n}])//matform[]]; Pr[""]; Pr["Then sort this - the sorting is by the first entry of"]; Pr["the pairs, so the random ordering of the first entries"]; Pr["is converted to a random ordering of the second entries"]; Pr[""]; Pr[" s = Sort[t] ="]; Pr[""]; Pr[" ",(s = Sort[t])//matform[]]; Pr[""]; Pr["Take the second entries"]; Pr[""]; Pr[" ans = Last/@s ="]; Pr[""]; Pr[" ", ans = Last/@s]; Pr[""]; Pr["This is the output returned"]; ans ]//Catch; (* Trap and warn of any other faulty inputs*) randpermX[e___] := (msg["argm",Length[{e}],1];fld); (* rats *) rats[x_] := Rationalize[x,10^-6]; (* ref *) (*Arrange for failure to result in return of the input unaltered and with error messages messages to be issued*) ref[x___]:= With[{result = Block[{fn = ref},refX[x]] }, result/; result =!= fld ]; refX[A_,opts___] := Module[ {posnopt,expln,zero,zeroQ,tzQ,Pr,PrW,pat,mf, showpiv,R,pivcol,c,rowredQ }, (*Checks inputs and issue any error messages needed*) Which[ !MatrixQ[A], msg["notmat",1];thfld, !((posnopt = Position[{opts}, a_/;!OptionQ[a],{1},1,Heads->False])=== {}), Message[Strang::"nonopt", {opts}[[posnopt[[1,1]]]], 1, fn];thfld ]; (*Read options*) {expln, zero} = {Explain, ZeroTest}/.{opts}/.Options[ref]; zeroQ = TrueQ[zero[#]]&; tzQ = tozeroQ[zero, changed];(* see notes on tozeroQ*) changed = False; (*Guard against receiving unasked value.*) (*Turn on explanation if required.*) If[MatchQ[expln,True|PatternOnly], Pr = Print;PrW = printwide, Pr = PrW = Hold ]; If[expln === PatternOnly, pat = Map[If[MatchQ[#,0|0.|1|1.0],#,"x"]&,#,{2}]&, pat = Identity ]; mf = matform[TableSpacing -> {0,3}]; (*Describe the input matrix*) If[expln === PatternOnly, Pr["We find the Reduced Echelon form of a matrix with pattern"], Pr["We find the Reduced Echelon form of the matrix "] ]; Pr[""]; PrW[" ", pat[A]//mf]; Pr[""]; (*Begin computation.*) (*Take U part and pivcol from the plufactorization of A *) Pr["Take the U part of its plu factorization:"]; {R,pivcol} = plu[A, ZeroTest->zero][[{3,4}]]; Pr[""]; PrW[" ", pat[R]//mf]; Pr[""]; (*Scale the rows to make all pivots equal 1*) showpiv[pos_] := ( Pr[""]; PrW[" ",MapAt[StringForm["<``>",#]&,pat[R], {k,c}]//mf]; Pr[""] ); Do[(*for k = 1 to length of pivcol*) c = pivcol[[k]]; Pr["Find the pivot in row ",k,"."]; showpiv[{k,c}]; If[ !TrueQ[R[[k,c]] == 1], Pr["Make this pivot equal 1, by scaling the row."]; R[[k]] = R[[k]]/R[[k,c]]; showpiv[{k,c}] ]; (*Eliminate non-zeros above the pivots*) If[k>1, changed = False; (*changed will be set to True if tzQ replaces any nonzero entry with zero*) rowredQ = False; (*-- rowredQ will be set to True if there are any row reductions*) Do[ If[ !tzQ[R[[i,c]]], rowredQ = True; R[[i]] = R[[i]] - R[[i,c]] R[[k]] ], {i,k-1} ]; Which[ rowredQ&&changed, Pr["Make the entries above the pivot zero by simplification and"] Pr["row operations."]; showpiv[{k,c}], rowredQ, Pr["Make the entries above the pivot zero by row operations."]; showpiv[{k,c}], changed, Pr["After some simplification the entries above the pivot are all"] Pr["zero."]; showpiv[{k,c}], True, Pr["The entries above the pivot are already zero."] ] ], {k,Length[pivcol]} ]; Pr["This is the required reduced echelon form, R."]; (*Computation completed.*) (*Describe the output.*) Pr["The output is the pair {R,pivcol}; where pivcol is the list of"]; Pr["the positions of the pivot columns of U."]; {R,pivcol} ]//Catch; (*trap and warn of any remaining kind of bad inputs*) refX[e___] := (msg["argt",Length[{e}],1,2];fld); (* signperm *) signperm[x___]:= With[{ result = Block[{fn = signperm}, signpermX[x]]}, result/;result =!= fld ]; signpermX[P_,opts___] := Module[{posnopt, expln,Pr,s,p,pi,pos}, Which[ !VectorQ[P,IntegerQ[#]&&Positive[#]&], msg["notposintlst",1];thfld, (posnopt = Position[{opts},a_/;!OptionQ[a],{1},1,Heads->False])=!= {}, msg["notopt", {opts}[[posnopt[[1,1]]]], 2];thfld ]; (*Read options.*) {expln} = {Explain}/.{opts}/.Options[signperm]; (*If Explain is set to True, turn on the printing.*) If[expln, Pr = Print, Pr = Hold]; Pr["The given list is"]; Pr[""]; Pr[" ", P]; Pr[""]; Pr["If it is not a permutation our calculation will find it out."]; Pr["Call it P with entries p1, p2, ..."]; Pr["We scan P from left to right and if pi is not i we look for"]; Pr["a j such that pj = i."]; Pr["If we find one then we swap the entries in positions i and j,"] Pr["and, if we are not at the end of P, we continue the scan. "]; Pr["If we do not find one then P is not a permutation (why?)."]; Pr["The answer is -1^ number of swaps used"]; Pr[""]; Pr["The changes in P starting from the input value are shown below."]; Pr[""]; Pr[" ",P]; s=1; (*initialize s*) p = P; (*give P a local name*) Do[ If[ (pi = p[[i]]) =!= i, pos = Position[p, i, {1},1]; If[pos === {}, Pr[""]; Pr["There is no pi with value ",i, " so P is not a permutation."]; msg["notperm",1,i,Length[p]];thfld ]; k = pos[[1,1]]; {p[[i]], p[[k]]} = {i,pi}; Pr[" ",p]; s=-s ], {i,Length[p]} ]; Pr[""]; Pr["The answer is ", s]; Pr["This is the output returned"]; s ]//Catch; signpermX[e___] := (msg["argm",Length[{e}],1];fld); (* slu *) (* in case of faulty input or failure of computation return input unaltered and display messages *) slu[x___]:= With[{ result = Block[{fn = slu}, sluX[x]]}, result/;result =!= fld ]; sluX[A_,opts___]:= Module[{posnopt,m,n,L,U,sign,expln,zero,tzQ,Pr,PrW,tf2,pat,prech}, (*input checks and error messages*) Which[ !MatrixQ[A], msg["notmat",1];thfld, !Equal@@Dimensions[A], msg["notsqrmattry",1];thfld, !((posnopt = Position[{opts}, a_/;!OptionQ[a],{1},1,Heads->False])=== {}), msg["notopt", {opts}[[posnopt[[1,1]]]], 1, fn];thfld ]; {m,n} = Dimensions[A]; (* initialize quantities to be computed*) L = IdentityMatrix[m]; U = A; (*Read options*) {expln, zero} = {Explain, ZeroTest}/.{opts}/.Options[slu]; tzQ = tozeroQ[zero, changed];(* see notes on tozeroQ*) changed = False; (*Guard against receiving un-asked-for value.*) (*If Explain is set to True, turn on printing and define the printing format of the tables used.*) If[MatchQ[expln,True|PatternOnly], Pr = Print;PrW = printwide; tf2[x___] := TableForm[{x}, TableDirections->{Row,Column,Row}, TableSpacing ->{8,0,3}, TableHeadings -> {SequenceForm[ SpaceForm[Floor[(3m-1)/2]],#]&/@{"L\n","U"},None,None}, TableAlignments -> If[MatrixQ[ Last[{x}],MatchQ[#,_Integer|_StringForm|"[0]"]& ], {Center,Center,Right}, {Center,Center,Center} ] ], Pr = PrW = tf2 = Hold ]; If[expln === PatternOnly, pat = Map[If[MatchQ[#,0|0.|1|1.0],#,"+"]&,#,{2}]&, pat = Identity ]; (*state input,and kind of commentary to be given -- explanations to be printed all begin with Pr and are fairly self-explanatory *) If[expln === PatternOnly, Pr["The pattern of the given matrix is"], Pr["The given matrix is"] ]; Pr[""]; Pr[" ",pat[U]//matform[TableSpacing -> {0,3}]]; Pr[""]; Pr["Call it A."]; Pr["We try to construct square matrices L, U such that"]; Pr["A = L.U and"]; Pr[" L is lower triangular with ones on the diagonal;"]; Pr[" U is upper triangular."]; Pr[""]; Pr["The construction starts with the pair below;the \ncommentary will be concerned with U; consequential changes \nin L are shown without comment" ]; Pr[""]; PrW[tf2[L,pat[U]]]; Pr[""]; Do[(*For k from 1 to m *) changed = False; (*changed will be set to True if tzQ,possibly inside findpiv, replaces any nonzero entry with zero*) If[ !tzQ[U[[k,k]]], (*-- if there is a pivot at {k,k}*) Pr[""]; Pr["Pivot ",k," is at ",{k,k}," in U."]; Pr[""]; PrW[tf2[pat[L], MapAt[StringForm["<``>",#]&,pat[U],{k,k}]]]; Pr[""]; (* Check the entries below the pivot. Assign value zero to any that simplify to zero. Row reduce the others to zero *) rowredQ = False; (* -- rowredQ records if any reductions were needed.*) Do[ If[!tzQ[U[[i,k]]], rowredQ = True; prech = changed;(*note the values of changed*) With[{mult = U[[i,k]]/U[[k,k]]}, (*Subtract mult times row k from row i.*) U[[i]] = U[[i]] - mult U[[k]]; tzQ[U[[i,k]]]; (*to deal with rounding problems*) (*Set L[[i,k]] to mult: the effect of adding mult times column i to column k. *) L[[i,k]] = mult ]; changed = prech(*restore the values of changed*) ], {i,k+1,m} ]; (*Explain, if k < m*) If[ k",#]&,pat[U],{k,k}]]]; Pr[""] ], (*Else*) (*Explain that the method has failed*) If[changed, Pr["But, after simplification,"], Pr["But"], ]; Pr[" U"[k,k]," = 0;"]; Pr["so there is no pivot at ",{k,k}," and this algorithm fails." ]; Pr["You might like to try splu or plu."]; (*Warn of failure*) msg["fldalgo", "splu or plu"];thfld ], {k,1,m } ] (*Explain that the calculation is complete.*) Pr[""]; Pr["The matrix U is upper triangular."]; Pr["The calculation is complete."]; Pr[""]; (*Describe the form of the output.*) Pr["The output is the list {L,U}."]; {L,U} (*-- The output.*) ]//Catch; (*Trap and warn of any other faulty inputs.*) sluX[e___] := (msg["argm", Length[{e}],1];fld); (* slv *) (*Failure must return the input unaltered*) slv[x___]:= With[{result = Block[{fn = slv},slvX[x]]}, result/;result =!= fld ]; slvX[A_,b_,opts___] := Module[{posnopt,m,n,expln,zero,Pr,PrW,pat,temp,L,U,pivcol,r,pb,c,x,j}, (*input checks and error messages*) Which[ !MatrixQ[A], msg["notmat",1];thfld, {m,n} = Dimensions[A]; !VectorQ[b]||Length[b] =!= m, If[m == 1, msg["rhs1", m], msg["rhs",m] ];thfld, !((posnopt = Position[{opts}, a_/;!OptionQ[a],{1},1,Heads->False])=== {}), msg["notopt", {opts}[[posnopt[[1,1]]]], 2];thfld ]; (*Read options.*) {expln, zero} = {Explain, ZeroTest}/.{opts}/.Options[slv]; (*Turn on explanation if required.*) If[MatchQ[expln,True|PatternOnly], Pr = Print;PrW = printwide, Pr = PrW = Hold ]; If[expln === PatternOnly, pat = Map[If[MatchQ[#,0|0.|1|1.0],#,"+"]&,#,{2}]&, pat = Identity ]; top[x_] := Prepend[Table[{""},{m-1}], {x}]; Pr["We try to solve A.x = b."]; If[expln === PatternOnly, Pr["The pattern of which is"]; ]; Pr[""]; PrW[" ", MapThread[Join,{pat[A], top[". x = "],pat[List/@b]}]// matform[] ]; Pr[""]; Pr["by using the slu factorization of A, A = L.U to convert this"]; Pr["first to"]; Pr[""]; Pr[" ", "L.U.x = P.b"]; Pr[""]; If[Head[temp = slu[A, ZeroTest -> zero]] === slu, Pr["The algorithm used has failed, you might like to try splv."]; msg["fldalgo", "splv or solve"];thfld, {L,U} = temp ]; Pr["and then to the pair"]; Pr[""]; Pr[" ", "L.c = b "]; Pr[" ", "U.x = c "]; Pr[""]; Pr["First solve L.c = b, that is"]; Pr[""]; PrW[" ", MapThread[Join,{pat[L], top[". c = "],pat[List/@b]}]// matform[] ]; Pr[""]; Pr["L is lower triangular with diagonal entries all equal to 1;"]; Pr["so we can solve by forward substitution:"]; Pr[""]; (* Solve L.c = pb by forward elimination. The diagonal elements of L are all equal to 1. *) c = Table[0,{m}];(*Initialize c.*) Do[c[[k]] = b[[k]] - L[[k]].c, {k,m}]; Pr[" ","c = ", c]; Pr[""]; Pr["Next, solve U.x = c, that is"]; Pr[""]; PrW[" ", MapThread[Join,{pat[U], top[". x = "],pat[List/@c]}]// matform[] ]; Pr[""]; Pr["U is upper triangular, with diagonal entries all non-zero;"]; Pr["and can be solved by back substitution:"]; (* Solve U.x = x by back substitution. The diagonal elements of U are all non-zero. *) x = Table[0,{n}]; Do[ x[[k]] = (c[[k]] - U[[k]].x )/U[[k,k]], {k,m,1,-1} ]; Pr[""]; PrW[" x = ", x]; Pr[""]; Pr["This is the unique solution to A.x = b."]; Pr["It is the output returned."]; x ]//Catch; (* Trap and warn of any other faulty inputs*) slvX[e___] := (msg["argm",Length[{e}],2];fld); (* splu *) (* in case of faulty input or failure of computation return input unaltered and display messages *) splu[x___]:= With[{ result = Block[{fn = splu}, spluX[x]]}, result/;result =!= fld ]; spluX[A_,opts___]:= Module[{posnopt,m,n,P,L,U,sign,expln,zero,tzQ,Pr,PrW,tf3,pat,pos,r,c,prech}, (*input checks and error messages*) Which[ !MatrixQ[A], msg["notmat",1];thfld, !Equal@@Dimensions[A], msg["notsqrmattry",1];thfld, !((posnopt = Position[{opts}, a_/;!OptionQ[a],{1},1,Heads->False])=== {}), msg["notopt", {opts}[[posnopt[[1,1]]]], 1, fn];thfld ]; {m,n} = Dimensions[A]; (* initialize quantities to be computed*) P = L = IdentityMatrix[m]; U = A; sign = 1; (*Read options*) {expln, zero} = {Explain, ZeroTest}/.{opts}/.Options[splu]; tzQ = tozeroQ[zero, changed];(* see notes on tozeroQ*) changed = False; (*Guard against receiving un-asked-for value.*) (*If Explain is set to True, turn on printing and define the printing format of the tables used.*) If[MatchQ[expln,True|PatternOnly], Pr = Print;PrW = printwide; tf3[x___] := TableForm[{x}, TableDirections->{Row,Column,Row}, TableSpacing ->{8,0,3}, TableHeadings -> {SequenceForm[ SpaceForm[Floor[(3m-1)/2]],#]&/@{"P\n","L","U"},None,None}, TableAlignments -> If[TrueQ[MatrixQ[ {x}[[3]],MatchQ[#,_Integer|_StringForm|"[0]"]&]], {Center,Center,Right}, {Center,Center,Center} ] ], Pr = PrW = tf3 = Hold ]; If[expln === PatternOnly, pat = Map[If[MatchQ[#,0|0.|1|1.0],#,"+"]&,#,{2}]&, pat = Identity ]; (*State input,and kind of commentary to be given -- explanations to be printed all begin with Pr and are fairly self-explanatory.*) If[expln === PatternOnly, Pr["The pattern of the given matrix is"], Pr["The given matrix is"] ]; Pr[""]; PrW[" ",pat[U]//matform[TableSpacing -> {0,3}]]; Pr[""]; Pr["Call it A."]; Pr["We try to construct square matrices P, L, U such that"]; Pr["P.A = L.U and"]; Pr[" P is a permutation matrix;"]; Pr[" L is lower triangular with ones on the diagonal;"]; Pr[" U is upper tringular."]; Pr[""]; Pr["The construction starts with the triple below: the \ncommentry will be concerned with U; consequential changes \nin P and L are shown without comment"]; Pr[""]; PrW[tf3[P,L,pat[U]]]; Pr[""]; (*For k from 1 to m look for a pivot in column k*) Do[ changed = False; (*changed will be set to True if tzQ,possibly inside findpiv, replaces any nonzero entry with zero*) {pos,U} = findpiv[U, k, Min[1,n], ZeroTest ->zero]; (*-- Pass zero into findpiv, and pass any changes to the U in findpiv out to the U on the left side.*) If[ !Less@@pos, (*if there is a suitable pivot in column k*) (*Name the row and column indices.*) {r,c} = pos; (*in fact c is necessarily equal to k; but keeping the notation c may make comparison with the code for plu easier*) Pr[""]; Pr["Pivot ",k," in U is at ",{r,c},"."]; Pr[""]; PrW[tf3[ P, pat[L], MapAt[StringForm["<``>",#]&,pat[U],{r,c}] ]]; Pr[""]; (*If necessary swap rows to put this pivot on the diagonal.*) If[r > k, {U[[k]],U[[r]]} = {U[[r]],U[[k]]}; (*make the consequential adjustments to P and L*) (*Swap rows r and k in P*) {P[[k]],P[[r]]} = {P[[r]],P[[k]]}; (*Swap rows r and k in L - but only up to position k-1.*) Do[ {L[[k,i]],L[[r,i]]} = {L[[r,i]],L[[k,i]]}, {i,k-1} ]; (*change sign to record the swap*) sign = -sign; Pr["Swap rows ",k," and ",r," to change it to ",{k,c},"."]; Pr[""]; PrW[tf3[ P,pat[L], MapAt[StringForm["<``>",#]&,pat[U],{k,c}] ]]; Pr[""]; Pr["Pivot ",k," is now at ",{k,c},"."]; ]; (* Check the entries below the pivot. Assign value zero to any that simplify to zero. Row reduce the others to zero *) rowredQ = False; (* -- rowredQ records if any reductions were needed.*) Do[ If[!tzQ[U[[i,c]]], rowredQ = True; prech = changed;(*note the values of changed*) With[{mult = U[[i,c]]/U[[k,c]]}, (* Subtract mult = U[[i,c]]/U[[k,c]] times row k from row i. *) U[[i]] = U[[i]] - mult U[[k]]; tzQ[U[[i,c]]]; (*to deal with rounding problems*) (*Set L[[i,k]] to mult: the effect of adding mult times column i to column k. *) L[[i,k]] = mult ]; changed = prech(*restore the value of changed*) ], {i,k+1,m} ]; (*Explain, if k < m*) If[ k",#]&,pat[U],{k,c}] (* MapAt[StringForm["[``]",#]&,pat[U],{#,c}&/@Range[k+1,m]] *) ]]; Pr[""] ], (*Else*) (*Explain that the method has failed*) If[ changed, Pr["But, after simplification, we have:"]; Pr[""]; PrW[tf3[P,pat[L],pat[U]]]; Pr[""]; Pr["so there is no pivot available in positions ",k, " to ",m ]; Pr["in column ",k," of U, and this algorithm fails."], Pr["But there is no pivot available in positions ",k," to " ,m ]; Pr["in column ",k," of U, and this algorithm fails."] ]; Pr["You might like to try plu."]; (*Warn of failure*) msg["fldalgo", "plu"];thfld ], {k,1,m } ]; (*Explain that the calculation is complete.*) Pr[""]; Pr["The calculation is complete."]; Pr[""]; (*Describe the form of the output.*) Pr["The output is the list {P,L,U,sign} where "]; Pr["sign is -1^(number of swaps used)."]; {P,L,U,sign} (*-- The output.*) ]//Catch; (*Trap and warn of any other faulty inputs.*) spluX[e___] := (msg["argm", Length[{e}],1];fld); (* splv *) (*Failure must return the input unaltered*) splv[x___]:= With[{result = Block[{fn = splv},splvX[x]]}, result/;result =!= fld ]; splvX[A_,b_,opts___] := Module[{posnopt,m,n,explnQ,zero,Pr,PrW,pat,temp,P,L,U,r,pb,c,x}, (*input checks and error messages*) Which[ !MatrixQ[A], msg["notmat",1];thfld, {m,n} = Dimensions[A]; !VectorQ[b]||Length[b] =!= m, If[m == 1, msg["rhs1", m], msg["rhs",m] ];thfld, !((posnopt = Position[{opts}, a_/;!OptionQ[a],{1},1,Heads->False])=== {} ), msg["notopt", {opts}[[posnopt[[1,1]]]], 2];thfld ]; (*Read options.*) {expln, zero} = {Explain, ZeroTest}/.{opts}/.Options[splv]; (*Turn on explanation if required.*) If[MatchQ[expln,True|PatternOnly], Pr = Print;PrW = printwide, Pr = PrW = Hold ]; If[expln === PatternOnly, pat = Map[If[MatchQ[#,0|0.|1|1.0],#,"+"]&,#,{2}]&, pat = Identity ]; top[x_] := Prepend[Table[{""},{m-1}], {x}]; Pr["We try to solve A.x = b."]; If[expln === PatternOnly, Pr["The pattern of which is"]; ]; Pr[""]; PrW[" ", MapThread[Join,{pat[A], top[". x = "],pat[List/@b]}]// matform[] ]; Pr[""]; Pr["by using the splu factorization of A, P.A = L.U to convert it"]; Pr["first to"]; Pr[""]; Pr[" ", "L.U.x = P.b"]; Pr[""]; If[Head[temp = splu[A, ZeroTest -> zero]] === splu, Pr["The algorithm used has failed, you might like to try solve."]; msg["fldalgo", "solve"];thfld, {P,L,U} = Take[temp,3] ]; pb = P.b; Pr["and then to the pair"]; Pr[""]; Pr[" ", "L.c = P.b "]; Pr[" ", "U.x = c "]; Pr[""]; Pr["First solve L.c = P.b, that is"]; Pr[""]; PrW[" ", MapThread[Join,{pat[L], top[". c = "],pat[List/@pb]}]// matform[] ]; Pr[""]; Pr["L is lower triangular with diagonal entries all equal to 1;"]; Pr["so we can solve by forward substitution:"]; Pr[""]; (* Solve L.c = pb by forward elimination. The diagonal elements of L are all equal to 1. *) c = Table[0,{m}];(*Initialize c.*) Do[c[[k]] = pb[[k]] - L[[k]].c, {k,m}]; Pr[" ","c = ", c]; Pr[""]; Pr["Next, solve U.x = c, that is"]; Pr[""]; PrW[" ", MapThread[Join,{pat[U], top[". x = "],pat[List/@c]}]// matform[] ]; Pr[""]; Pr["U is upper triangular, with diagonal entries all non-zero"]; Pr["and can be solved by back substitution."]; (* Solve U.x = x by back substitution. The diagonal elements of U are all non-zero. *) x = Table[0,{n}];(*Initialize x.*) Do[ x[[k]] = (c[[k]] - U[[k]].x )/U[[k,k]], {k,m,1,-1} ]; Pr[""]; PrW[" x = ", x]; Pr[""]; Pr["This is the unique solution to A.x = b."]; Pr["It is the output returned."]; x ]//Catch; (* Trap and warn of any other faulty inputs*) splvX[e___] := (msg["argm",Length[{e}],2];fld); (* solve *) (*Arrange for failure to result in return of the input unaltered and with error messages messages to be issued*) solve[x___]:= With[{result = Block[{fn = solve},solveX[x]]}, result/;result =!= fld ]; solveX[A_,b_,opts___] := Module[{posnopt,m,n,explnQ,zero,tzQ,Pr,PrW,pat,P,L,U,pivcol,r,pb,c,x,j}, (*input checks and error messages*) Which[ !MatrixQ[A], msg["notmat",1];thfld, {m,n} = Dimensions[A]; !VectorQ[b]||Length[b] =!= m, If[m == 1, msg["rhs1", m], msg["rhs",m] ];thfld, !((posnopt = Position[{opts}, a_/;!OptionQ[a],{1},1,Heads->False])=== {}), msg["notopt", {opts}[[posnopt[[1,1]]]], 2];thfld ]; (*Read options.*) {expln, zero} = {Explain, ZeroTest}/.{opts}/.Options[solve]; tzQ = tozeroQ[zero,changed];(*See notes on tozeroQ.*) changed = False; (*Guard against receiving un-asked-for value.*) (*Turn on explanation if required.*) If[MatchQ[expln,True|PatternOnly], Pr = Print;PrW = printwide, Pr = PrW = Hold ]; If[expln === PatternOnly, pat = Map[If[MatchQ[#,0|0.|1|1.0],#,"+"]&,#,{2}]&, pat = Identity ]; top[x_] := Prepend[Table[{""},{m-1}], {x}]; (*Begin explanation*) Pr["We try to find a particular solution of the system A.x = b."]; If[expln === PatternOnly, Pr["The pattern of which is"]; ]; Pr[""]; PrW[" ", MapThread[Join,{pat[A], top[". x = "],pat[List/@b]}]// matform[] ]; Pr[""]; Pr["Use the plu factorization of A, P.A = L.U we convert this"]; Pr["first to"]; Pr[""]; Pr[" ", "L.U.x = P.b"]; {P,L,U,pivcol} = Take[plu[A, ZeroTest -> zero],4]; pb = P.b; Pr[""]; Pr["and then to the pair"]; Pr[""]; Pr[" ", "L.c = P.b "]; Pr[" ", "U.x = c "]; Pr[""]; Pr["First solve L.c = P.b"]; Pr[""]; PrW[" ", MapThread[Join,{pat[L], top[". c = "],pat[List/@pb]}]// matform[] ]; Pr[""]; Pr["L is lower triangular with diagonal entries all equal to 1;"]; Pr["so we can solve by forward substitution:"]; Pr[""]; (* Forward eliminate to solve L c = pb. Note that L is lower triangular with 1©s on the diagonal. *) c = Table[0,{m}]; (*Initialize c.*) Do[c[[k]] = pb[[k]] - L[[k]].c, {k,m}]; Pr[" ","c = ", c]; Pr[""]; r = Length[pivcol]; (* Check if equations are consistent. The last m-r entries of the new right side c should be zero. *) Pr["Next try to solve U.x = c, but check first for consistency."]; Pr[""]; PrW[" ", MapThread[Join,{pat[U], top[". x = "],pat[List/@c]}]// matform[] ]; Pr[""]; Which[ m-r == 1, Pr[" ","the last entry in c should be zero."], m-r > 1, Pr[" ","the last ",m-r, " entries in c should be zero."] ]; (*Use txQ to Find if there is an i > r for which c[[i]] does not simplify to zero. Stop if one is found. *) changed = False; (* changed will be set to zero if tzQ replaces an non-zero element by zero*) If[ Position[ Range[r+1, m], i_/;!tzQ[c[[i]]], {1}, 1, Heads->False ] =!= {}, Pr[" ","The equation is inconsistent:"]; msg["incons"];thfld, (*else*) If[ changed, Pr[" ", "Simplification shows that the equation is in fact consistent"]; Pr[""]; PrW[" ", pat[U]//matform[]," .x = ", pat[List/@c]//matform[]]; Pr[""], Pr["-- clearly consistent."] ]; ]; (* If the equations are consistent, use back substitution to find the particular solution of U x = c with the entries of x that are not in the positions of the pivot columns being equal to zero. *) Pr[""]; Pr["U is upper triangular, and can be solved by back substitution."]; Pr["In fact we find the unique solution which has all the entries"]; Pr["corresponding to non-pivot columns of U equal to zero;it is"]; x = Table[0,{n}]; Do[ j = pivcol[[k]]; x[[j]] = (c[[k]] - U[[k]].x )/U[[k,j]], {k,r,1,-1} ]; Pr[""]; PrW[" x = ", x]; Pr[""]; Pr["This is a particular solution to A.x = b."]; Pr["The output returned is the value of x"]; x ]//Catch; (* Trap and warn of any other faulty inputs*) solveX[e___] := (msg["argm",Length[{e}],2];fld); (*Finish Private*) End[]; (*Finish*) Protect["`*"]; EndPackage[];