subroutine mat_help_text()
G_HELP_TEXT=[ CHARACTER(LEN=128) :: &
'================================================================================',&
'LALA USERS'' GUIDE',&
'',&
' LALA (May, 1981 - Apr, 2021)',&
'',&
' The Linear Algebra Fortran Facility (LALA) is a collection of Fortran',&
' procedures that serves as a convenient tool for Fortran programs to',&
' interact with their data (interactively or in batch mode) with a',&
' tool that acts as a basic "laboratory" for computations involving',&
' matrices.',&
'',&
' It provides easy access to matrix software developed by the LINPACK',&
' and EISPACK projects.',&
'',&
' It is based on the Los Alamos procedure MATLAB, and owes much to',&
' Cleve Moler, Department of Computer Science, University of New Mexico.',&
'',&
' CONTENTS',&
'',&
' - Elementary operations',&
' - LALA functions',&
' - Rows, columns and submatrices',&
' - "for", "while" and "if"',&
' - Characters, text, files and macros',&
' - The numerical algorithms',&
' - "flop" and "chop"',&
' - Census example',&
' - Partial differential equation example',&
' - Eigenvalue sensitivity example',&
' - Communicating with other programs',&
' - Appendix (The HELP file)',&
'',&
' The capabilities range from standard tasks such as solving simultaneous',&
' linear equations and inverting matrices, through symmetric and',&
' nonsymmetric eigenvalue problems, to fairly sophisticated matrix',&
' tools such as the singular value decomposition.',&
'',&
' LALA should be useful in applied linear algebra, as well as more',&
' general numerical analysis, matrix theory, statistics and applications',&
' of matrices to other disciplines.',&
'',&
' LALA can serve as a "desk calculator" for the quick solution of small',&
' problems involving matrices.',&
'',&
' The program is written in Fortran and is designed to be readily',&
' installed under any operating system which permits interactive',&
' execution of Fortran programs. The resources required are fairly',&
' modest.',&
'',&
' The size of the matrices that can be handled in LALA depends upon the',&
' amount of storage available on the supporting platform and the optional',&
' word count that can be supplied on an initial call to LALA(3f).',&
'',&
' In some ways, LALA resembles SPEAKEASY [4] and, to a lesser extent,',&
' APL. All are interactive terminal languages that ordinarily accept',&
' single-line commands or statements, process them immediately, and print',&
' the results. All have arrays or matrices as principal data types. But',&
' for LALA, the matrix is the only data type (although scalars, vectors',&
' and text are special cases), the underlying system is portable and',&
' requires fewer resources, and the supporting subroutines are more',&
' powerful and in some cases, have better numerical properties.',&
'',&
' Together, LINPACK and EISPACK provide for powerful matrix',&
' computation. EISPACK is a package of over 70 Fortran subroutines for',&
' various matrix eigenvalue computations that are based for the most',&
' part on Algol procedures published by Wilkinson, Reinsch and their',&
' colleagues [5]. LINPACK is a package of 40 Fortran subroutines (in',&
' each of four data types) for solving and analyzing simultaneous linear',&
' equations and related matrix problems. Since LALA is not primarily',&
' concerned with either execution time efficiency or storage savings,',&
' it ignores most of the special matrix properties that LINPACK and',&
' EISPACK subroutines use to advantage. Consequently, only 8 subroutines',&
' from LINPACK and 5 from EISPACK are actually involved.',&
'',&
' In more advanced applications, LALA can be used in conjunction with',&
' other programs in several ways. It is possible to define new LALA',&
' functions and add them to the system. it is possible to use the local',&
' file system to pass matrices between LALA and other programs. LALA',&
' command and statement input can be obtained from a local file instead',&
' of from the terminal. The most power and flexibility is obtained by',&
' using LALA as a subroutine which is called by other programs.',&
'',&
' This document first gives an overview of LALA from the user''s',&
' point of view. Several extended examples involving data fitting,',&
' partial differential equations, eigenvalue sensitivity and other',&
' topics are included. The system was designed and programmed using',&
' techniques described by Wirth [6], implemented in nonrecursive,',&
' portable Fortran. There is a brief discussion of some of the matrix',&
' algorithms and of their numerical properties. A final section describes',&
' how LALA can be used with other programs. The appendix includes the',&
' HELP documentation available on-line.',&
'',&
'================================================================================',&
'ELEMENTARY OPERATIONS',&
'',&
' LALA works with essentially only one kind of object, a rectangular',&
' matrix with complex elements. If the imaginary parts of the elements',&
' are all zero, they are not printed, but they still occupy storage. In',&
' some situations, special meaning is attached to 1 by 1 matrices,',&
' that is scalars, and to 1 by n and m by 1 matrices, that is row and',&
' column vectors.',&
'',&
' Matrices can be introduced into LALA in four different',&
' ways:',&
'',&
' -- Explicit list of elements,',&
' -- Use of "for" and "while" statements,',&
' -- Read from an external file,',&
' -- Execute an external Fortran program.',&
'',&
' The explicit list is surrounded by angle brackets, ''<'' and ''>'' or',&
' braces ''['' and '']'', and uses the semicolon '';'' to indicate the ends',&
' of the rows. For example, the input line',&
'',&
' A = <1 2 3; 4 5 6; 7 8 9>',&
'',&
' will result in the output',&
'',&
' A =',&
'',&
' 1. 2. 3.',&
' 4. 5. 6.',&
' 7. 8. 9.',&
'',&
' The matrix A will be saved for later use. The individual elements',&
' are separated by commas or blanks and can be any LALA expressions,',&
' for example',&
'',&
' x = < -1.3, 4/5, 4*atan(1) >',&
'',&
' results in',&
'',&
' x =',&
'',&
' -1.3000 0.8000 3.1416',&
'',&
' The elementary functions available include sqrt, log, exp, sin, cos,',&
' atan, abs, round, real, imag, and conjg.',&
'',&
' Large matrices can be spread across several input lines, with the',&
' carriage returns replacing the semicolons. The above matrix could',&
' also have been produced by',&
'',&
' A = < 1 2 3',&
' 4 5 6',&
' 7 8 9 >',&
'',&
' Matrices can be input from the local file system. Say a file named',&
' ''xyz'' contains five lines of text,',&
'',&
' A = <',&
' 1 2 3',&
' 4 5 6',&
' 7 8 9',&
' >;',&
'',&
' then the LALA statement exec(''xyz'') reads the matrix and assigns it',&
' to A .',&
'',&
' The "for" statement allows the generation of matrices whose elements',&
' are given by simple formulas. Our example matrix A could also have',&
' been produced by',&
'',&
' for i = 1:3, for j = 1:3, A(i,j) = 3*(i-1)+j;',&
'',&
' The semicolon at the end of the line suppresses the printing, which',&
' in this case would have been nine versions of A with changing elements.',&
'',&
' Several statements may be given on a line, separated by semicolons',&
' or commas.',&
'',&
' Two consecutive periods anywhere on a line indicate continuation. The',&
' periods and any following characters are deleted, then another line',&
' is input and concatenated onto the previous line.',&
'',&
' Two consecutive slashes anywhere on a line cause the remainder of',&
' the line to be ignored. This is useful for inserting comments.',&
'',&
' Names of variables are formed by a letter, followed by any of',&
' letters, digits and underscores, up to 63 characters in length.',&
'',&
' The special character prime ('') is used to denote the transpose of',&
' a matrix, so',&
'',&
' X = X''',&
'',&
' changes the row vector above into the column vector',&
'',&
' X =',&
'',&
' -1.3000',&
' 0.8000',&
' 3.1416',&
'',&
' Individual matrix elements may be referenced by enclosing their',&
' subscripts in parentheses. When any element is changed, the entire',&
' matrix is reprinted. For example, using the above matrix,',&
'',&
' B(3,3) = B(1,3) + B(3,1)',&
'',&
' results in',&
'',&
' B =',&
'',&
' 1. 2. 3.',&
' 4. 5. 6.',&
' 7. 8. 10.',&
'',&
' Addition, subtraction and multiplication of matrices are denoted by',&
' +, -, and * . The operations are performed whenever the matrices',&
' have the proper dimensions. For example, with the above A and x,',&
' the expressions A + X and X*A are incorrect because A is 3 by 3 and',&
' X is now 3 by 1. However,',&
'',&
' B = A*B',&
'',&
' is correct and results in the output',&
'',&
' B =',&
'',&
' 9.7248',&
' 17.6496',&
' 28.7159',&
'',&
' Note that both upper and lower case letters are allowed for input',&
' (on those systems which have both).',&
'',&
' There are two "matrix division" symbols in LALA, \ and / . If A and',&
' B are matrices, then A\B and B/A correspond formally to left and right',&
' multiplication of B by the inverse of A, that is inv(A)*B and B*inv(A),',&
' but the result is obtained directly without the computation of the',&
' inverse. In the scalar case, 3\1 and 1/3 have the same value, namely',&
' one-third. In general, A\B denotes the solution X to the equation A*X =',&
' B and B/A denotes the solution to X*A = B.',&
'',&
' Left division, A\B, is defined whenever B has as many rows as A. If A',&
' is square, it is factored using Gaussian elimination. The factors are',&
' used to solve the equations A*X(:,j) = B(:,j) where B(:,j) denotes the',&
' j-th column of B. The result is a matrix X with the same dimensions',&
' as B. If A is nearly singular (according to the LINPACK condition',&
' estimator, RCOND(3f)), a warning message is printed. If A is not',&
' square, it is factored using Householder orthogonalization with column',&
' pivoting. The factors are used to solve the under- or overdetermined',&
' equations in a least squares sense. The result is an M by N matrix X',&
' where M is the number of columns of A and N is the number of columns',&
' of B . Each column of X has at most K nonzero components, where K is',&
' the effective rank of A .',&
'',&
' Right division, B/A, can be defined in terms of left division by B/A =',&
' (A''\B'')''.',&
'',&
' For example, since our vector b was computed as A*X, the statement',&
'',&
' Y = A\B',&
'',&
' results in',&
'',&
' Y =',&
'',&
' -1.3000',&
' 0.8000',&
' 3.1416',&
'',&
' Of course, Y is not exactly equal to X because of the roundoff errors',&
' involved in both A*X and A\B , but we are not printing enough digits',&
' to see the difference. The result of the statement',&
'',&
' E = X - Y',&
'',&
' depends upon the particular computer being used. In one case it',&
' produces',&
'',&
' E =',&
'',&
' 1.0e-15 *',&
'',&
' .3053',&
' -.2498',&
' .0000',&
'',&
' The quantity 1.0e-15 is a scale factor which multiplies all the',&
' components which follow. Thus our vectors X and Y actually',&
' agree to about 15 decimal places on this computer.',&
'',&
' It is also possible to obtain element-by-element',&
' multiplicative operations. If A and B have the same dimensions,',&
' then A .* B denotes the matrix whose elements are simply the',&
' products of the individual elements of A and B . The expressions',&
' A ./ B and A .\ B give the quotients of the individual elements.',&
'',&
' There are several possible output formats. The statement',&
'',&
' long, X',&
'',&
' results in',&
'',&
' X =',&
'',&
' -1.300000000000000',&
' .800000000000000',&
' 3.141592653589793',&
'',&
' The statement',&
'',&
' short',&
'',&
' restores the original format.',&
'',&
' The expression A**p means A to the p-th power. It is',&
' defined if A is a square matrix and p is a scalar. If p is an',&
' integer greater than one, the power is computed by repeated',&
' multiplication. For other values of p the calculation involves',&
' the eigenvalues and eigenvectors of A.',&
'',&
' Previously defined matrices and matrix expressions can be',&
' used inside brackets to generate larger matrices, for example',&
'',&
' C = <A, B; <4 2 0>*X, X''>',&
'',&
' results in',&
'',&
' C =',&
'',&
' 1.0000 2.0000 3.0000 9.7248',&
' 4.0000 5.0000 6.0000 17.6496',&
' 7.0000 8.0000 10.0000 28.7159',&
' -3.6000 -1.3000 0.8000 3.1416',&
'',&
' There are four predefined variables, "eps", "flop", "rand" and',&
' "eye". The variable "eps" is used as a tolerance is determining such',&
' things as near singularity and rank. Its initial value is the distance',&
' from 1.0 to the next largest floating point number on the particular',&
' computer being used. The user may reset this to any other value,',&
' including zero. "eps" is changed by "chop", which is described later',&
' in this manual.',&
'',&
' The value of "rand" is a random variable, with a choice of a uniform',&
' or a normal distribution.',&
'',&
' The name "eye" is used in place of "i" to denote identity matrices',&
' because "i" is often used as a subscript or as sqrt(-1). The dimensions',&
' of "eye" are determined by context. For example,',&
'',&
' B = A + 3*eye',&
'',&
' adds 3 to the diagonal elements of A and',&
'',&
' X = eye/A',&
'',&
' is one of several ways in LALA to invert a matrix.',&
'',&
' "flop" provides a measure of the number of floating point operations,',&
' or "flops", required for each calculation by reporting the CPU time',&
' consumed.',&
'',&
' A statement may consist of an expression alone, in which case a',&
' variable named "ans" is created and the result stored in "ans" for',&
' possible future use. Thus',&
'',&
' A\A - eye',&
'',&
' is the same as',&
'',&
' ans = A\A - eye',&
'',&
' (Roundoff error usually causes this result to be a matrix of "small"',&
' numbers, rather than all zeros.)',&
'',&
' All computations are done using double precision real arithmetic. There',&
' is no mixed-precision arithmetic. The Fortran COMPLEX data type',&
' is not used because many systems create unnecessary underflows and',&
' overflows with complex operations.',&
'',&
'================================================================================',&
'FUNCTIONS',&
'',&
' Much of LALA''s computational power comes from the various',&
' matrix functions available. The current list includes:',&
'',&
' inv(A) - Inverse.',&
' det(A) - Determinant.',&
' cond(A) - Condition number.',&
' rcond(A) - A measure of nearness to singularity.',&
' eig(A) - Eigenvalues and eigenvectors.',&
' schur(A) - Schur triangular form.',&
' hess(A) - Hessenberg or tridiagonal form.',&
' poly(A) - Characteristic polynomial.',&
' svd(A) - Singular value decomposition.',&
' pinv(A,eps) - Pseudo-inverse with optional tolerance.',&
' rank(A,eps) - Matrix rank with optional tolerance.',&
' lu(A) - Factors from Gaussian elimination.',&
' chol(A) - Factor from Cholesky factorization.',&
' qr(A) - Factors from Householder orthogonalization.',&
' rref(A) - Reduced row echelon form.',&
' orth(A) - Orthogonal vectors spanning range of A.',&
' exp(A) - e to the A.',&
' log(A) - Natural logarithm.',&
' sqrt(A) - Square root.',&
' sin(A) - Trigonometric sine.',&
' cos(A) - Cosine.',&
' atan(A) - Arctangent.',&
' round(A) - Round the elements to nearest integers.',&
' abs(A) - Absolute value of the elements.',&
' real(A) - Real parts of the elements.',&
' imag(A) - Imaginary parts of the elements.',&
' conjg(A) - Complex conjugate.',&
' sum(A) - Sum of the elements.',&
' prod(A) - Product of the elements.',&
' diag(A) - Extract or create diagonal matrices.',&
' tril(A) - Lower triangular part of A.',&
' triu(A) - Upper triangular part of A.',&
' norm(A,p) - Norm with p = 1, 2 or ''Infinity''.',&
' eye(m,n) - Portion of identity matrix.',&
' rand(m,n) - Matrix with random elements.',&
' ones(m,n) - Matrix of all ones.',&
' magic(n) - Interesting test matrices.',&
' invh(n) - Inverse Hilbert matrices.',&
' roots(C) - Roots of polynomial with coefficients C.',&
' display(A,p) - Print base p representation of A.',&
' kron(A,B) - Kronecker tensor product of A and B.',&
' plot(X,Y) - Plot Y as a function of X .',&
' rat(A) - Find "simple" rational approximation to A.',&
' user(A) - Function defined by external program.',&
'',&
' Some of these functions have different interpretations when the',&
' argument is a matrix or a vector and some of them have additional',&
' optional arguments. Details are given in the HELP document in the',&
' appendix.',&
'',&
' Several of these functions can be used in a generalized assignment',&
' statement with two or three variables on the left hand side. For',&
' example',&
'',&
' <X,D> = eig(A)',&
'',&
' stores the eigenvectors of A in the matrix X and a diagonal matrix',&
' containing the eigenvalues in the matrix D. The statement',&
'',&
' eig(A)',&
'',&
' simply computes the eigenvalues and stores them in "ans".',&
'',&
' Future versions of LALA will probably include additional functions,',&
' since they can easily be added to the system.',&
'',&
'================================================================================',&
'ROWS COLUMNS AND SUBMATRICES',&
'',&
' Individual elements of a matrix can be accessed by giving their',&
' subscripts in parentheses, eg. A(1,2), x(i), TAB(ind(k)+1).',&
' An expression used as a subscript is rounded to the nearest integer.',&
'',&
' Individual rows and columns can be accessed using a colon '':'' (or a',&
' ''|'') for the free subscript. For example, A(1,:) is the first row of',&
' A and A(:,j) is the j-th column. Thus',&
'',&
' A(i,:) = A(i,:) + c*A(k,:)',&
'',&
' adds c times the k-th row of A to the i-th row.',&
'',&
' The colon is used in several other ways in LALA, but all of the uses',&
' are based on the following definition.',&
'',&
' j:k is the same as <j, j+1, ..., k>',&
' j:k is empty if j > k .',&
' j:i:k is the same as <j, j+i, j+2i, ..., k>',&
' j:i:k is empty if i > 0 and j > k or if i < 0 and j < k .',&
'',&
' The colon is usually used with integers, but it is possible to',&
' use arbitrary real scalars as well. Thus',&
'',&
' 1:4 is the same as <1, 2, 3, 4>',&
' 0: 0.1: 0.5 is the same as <0.0, 0.1, 0.2, 0.3, 0.4, 0.5>',&
'',&
' In general, a subscript can be a vector. If X and V are vectors,',&
' then X(V) is <X(V(1)), X(V(2)), ..., X(V(n))> . This can also be',&
' used with matrices. If V has m components and W has n components,',&
' then A(V,W) is the m by n matrix formed from the elements of A whose',&
' subscripts are the elements of V and W. Combinations of the colon',&
' notation and the indirect subscripting allow manipulation of various',&
' submatrices. For example,',&
'',&
' A(<1,5>,:) = A(<5,1>,:) interchanges rows 1 and 5 of A.',&
' A(2:k,1:n) is the submatrix formed from rows 2 through k',&
' and columns 1 through n of A .',&
' A(:,<3 1 2>) is a permutation of the first three columns.',&
'',&
' The notation A(:) has a special meaning. On the right hand side of an',&
' assignment statement, it denotes all the elements of A, regarded as',&
' a single column. When an expression is assigned to A(:), the current',&
' dimensions of A, rather than of the expression, are used.',&
'',&
'================================================================================',&
'FOR WHILE AND IF',&
'',&
' The "for" clause allows statements to be repeated a specific',&
' number of times. The general form is',&
'',&
' for variable = expr, statement, ..., statement, end',&
'',&
' The "end" and the comma before it may be omitted. In general, the',&
' expression may be a matrix, in which case the columns are stored one',&
' at a time in the variable and the following statements, up to the',&
' "end" or the end of the line, are executed. The expression is often',&
' of the form j:k, and its "columns" are simply the scalars from j to',&
' k. Some examples (assume n has already been assigned a value):',&
'',&
' for i = 1:n, for j = 1:n, A(i,j) = 1/(i+j-1);',&
'',&
' generates the Hilbert matrix.',&
'',&
' for j = 2:n-1, for i = j:n-1, ...',&
' A(i,j) = 0; end; A(j,j) = j; end; A',&
'',&
' changes all but the "outer edge" of the lower triangle and then',&
' prints the final matrix.',&
'',&
' for h = 1.0: -0.1: -1.0, (<h, cos(pi*h)>)',&
'',&
' prints a table of cosines.',&
'',&
' <X,D> = eig(A); for v = X, v, A*v',&
'',&
' displays eigenvectors, one at a time.',&
'',&
' The "while" clause allows statements to be repeated an',&
' indefinite number of times. The general form is',&
'',&
' while expr relop expr, statement,..., statement, end',&
'',&
' where relop is =, <, >, <=, >=, or <> (not equal). The statements are',&
' repeatedly executed as long as the indicated comparison between the',&
' real parts of the first components of the two expressions is true. Here',&
' are two examples. (Exercise for the reader: What do these segments do?)',&
'',&
' eps = 1;',&
' while 1 + eps > 1, eps = eps/2;',&
' eps = 2*eps',&
'',&
' E = 0*A; F = E + eye; n = 1;',&
' while norm(E+F-E,1) > 0, E = E + F; F = A*F/n; n = n + 1;',&
' E',&
'',&
' The IF clause allows conditional execution of statements. The general',&
' form is',&
'',&
' if expr relop expr, statement, ..., statement,',&
' else statement, ..., statement',&
'',&
' The first group of statements are executed if the relation is true and',&
' the second group are executed if the relation is false. The "else"',&
' and the statements following it may be omitted. For example,',&
'',&
' if abs(i-j) = 2, A(i,j) = 0;',&
'',&
'================================================================================',&
'CHARACTERS AND TEXTFILES AND MACROS',&
'',&
' LALA has several commands which control the output format and the',&
' overall execution of the system.',&
'',&
' The "help" command allows on-line access to short portions of text',&
' describing various operations, functions and special characters. The',&
' entire "help" document is reproduced in an appendix.',&
'',&
' Results are usually printed in a scaled fixed point format that shows',&
' 4 or 5 significant figures. The commands "short", "long", "short e",',&
' "long e" and "long z" alter the output format, but do not alter the',&
' precision of the computations or the internal storage.',&
'',&
' The "who" command provides information about the functions and',&
' variables that are currently defined.',&
'',&
' The "clear" command erases all variables, except "eps", "flop",',&
' "rand" and "eye". The statement A = <> indicates that a "0 by 0"',&
' matrix is to be stored in A. This causes A to be erased so that its',&
' storage can be used for other variables.',&
'',&
' The "quit" and "exit" commands cause return to the underlying operating',&
' system through the Fortran RETURN statement.',&
'',&
' LALA has a limited facility for handling text. Any string of characters',&
' delineated by quotes (with two quotes used to allow one quote within',&
' the string) is saved as a vector of integer values that are the ADE',&
' (Ascii Decimal Equivalent) value of the character, with special',&
' equivalencing of the characters {}[]" into ()<>'' in expressions. It',&
' is important to know you use those characters as part of an expression',&
' or command without treating them as equivalent outside of strings.',&
'',&
' (The complete list is in the appendix under "CHAR".) For example',&
'',&
' ''2*A + 3'' is the same as < 50 42 65 32 43 32 51 >.',&
'',&
' It is possible, though seldom very meaningful, to use such',&
' strings in matrix operations. More frequently, the text is used',&
' as a special argument to various functions.',&
'',&
' norm(A,''inf'') computes the infinity norm of A .',&
' display(T) prints the text stored in T .',&
' exec(''file'') obtains LALA input from an external file.',&
' save(''file'') stores all the current variables in a file.',&
' load(''file'') retrieves all the variables from a file.',&
' print(''file'',X) prints X on a file.',&
' diary(''file'') makes a copy of the complete LALA session.',&
'',&
' The text can also be used in a limited string substitution',&
' macro facility. If a variable, say T, contains the source text',&
' for a LALA statement or expression, then the construction',&
'',&
' > T <',&
'',&
' causes T to be executed or evaluated. For example',&
'',&
' T = ''2*A + 3'';',&
' S = ''B = >T< + 5''',&
' A = 4;',&
' > S <',&
'',&
' produces',&
'',&
' B =',&
'',&
' 16.',&
'',&
' Some other examples are given under MACROS in the appendix. This',&
' facility is useful for fairly short statements and expressions.',&
' More complicated LALA "programs" should use the "exec" facility.',&
'',&
'================================================================================',&
'NUMERICAL ALGORITHMS',&
'',&
' The algorithms underlying the basic LALA functions are described in',&
' the LINPACK and EISPACK guides [1-3]. The following list gives the',&
' subroutines used by these functions.',&
'',&
' inv(A) - CGECO,CGEDI',&
' det(A) - CGECO,CGEDI',&
' lu(A) - CGEFA',&
' rcond(A) - CGECO',&
' chol(A) - CPOFA',&
' svd(A) - CSVDC',&
' cond(A) - CSVDC',&
' norm(A,2) - CSVDC',&
' pinv(A,eps) - CSVDC',&
' rank(A,eps) - CSVDC',&
' qr(A) - CQRDC,CQRSL',&
' orth(A) - CQRDC,CSQSL',&
' A\B and B/A - CGECO,CGESL if A is square.',&
' - CQRDC,CQRSL if A is not square.',&
' eig(A) - HTRIDI,IMTQL2,HTRIBK if A is Hermitian.',&
' - CORTH,COMQR2 if A is not Hermitian.',&
' schur(A) - same as EIG.',&
' hess(A) - same as EIG.',&
'',&
' Minor modifications were made to all these subroutines. The LINPACK',&
' routines were changed to replace the Fortran complex arithmetic',&
' with explicit references to real and imaginary parts. Since most',&
' of the floating point arithmetic is concentrated in a few low-level',&
' subroutines which perform vector operations (the Basic Linear Algebra',&
' Subprograms), this was not an extensive change. It also facilitated',&
' implementation of the "flop" and "chop" features which count and',&
' optionally truncate each floating point operation.',&
'',&
' The EISPACK subroutine COMQR2 was modified to allow access to the',&
' Schur triangular form, ordinarily just an intermediate result. IMTQL2',&
' was modified to make computation of the eigenvectors optional. Both',&
' subroutines were modified to eliminate the machine-dependent accuracy',&
' parameter and all the EISPACK subroutines were changed to include',&
' "flop" and "chop".',&
'',&
' The algorithms employed for the "poly" and "roots" functions',&
' illustrate an interesting aspect of the modern approach to eigenvalue',&
' computation. "poly(A)" generates the characteristic polynomial of',&
' A and "roots(poly(A))" finds the roots of that polynomial, which',&
' are, of course, the eigenvalues of A . But both "poly" and "roots"',&
' use EISPACK eigenvalues subroutines, which are based on similarity',&
' transformations. So the classical approach which characterizes',&
' eigenvalues as roots of the characteristic polynomial is actually',&
' reversed.',&
'',&
' If A is an n by n matrix, "poly(A)" produces the coefficients C(1)',&
' through C(n+1), with C(1) = 1, in',&
'',&
' det(z*eye-A) = C(1)*z**n + ... + C(n)*z + C(n+1) .',&
'',&
' The algorithm can be expressed compactly using LALA:',&
'',&
' Z = eig(A);',&
' C = 0*ones(n+1,1); C(1) = 1;',&
' for j = 1:n, C(2:j+1) = C(2:j+1) - Z(j)*C(1:j);',&
' C',&
'',&
' This recursion is easily derived by expanding the product',&
'',&
' (z - z(1))*(z - z(2))* ... * (z-z(n)) .',&
'',&
' It is possible to prove that "poly(A)" produces the coefficients in',&
' the characteristic polynomial of a matrix within roundoff error of',&
' A. This is true even if the eigenvalues of A are badly conditioned. The',&
' traditional algorithms for obtaining the characteristic polynomial',&
' which do not use the eigenvalues do not have such satisfactory',&
' numerical properties.',&
'',&
' If C is a vector with n+1 components, "roots(C)" finds the roots of',&
' the polynomial of degree n ,',&
'',&
' p(z) = C(1)*z**n + ... + C(n)*z + C(n+1) .',&
'',&
' The algorithm simply involves computing the eigenvalues of the',&
' companion matrix:',&
'',&
' A = 0*ones(n,n)',&
' for j = 1:n, A(1,j) = -C(j+1)/C(1);',&
' for i = 2:n, A(i,i-1) = 1;',&
' eig(A)',&
'',&
' It is possible to prove that the results produced are the exact',&
' eigenvalues of a matrix within roundoff error of the companion matrix',&
' A, but this does not mean that they are the exact roots of a polynomial',&
' with coefficients within roundoff error of those in C . There are',&
' more accurate, more efficient methods for finding polynomial roots,',&
' but this approach has the crucial advantage that it does not require',&
' very much additional code.',&
'',&
' The elementary functions "exp", "log", "sqrt", "sin", "cos" and "atan"',&
' are applied to square matrices by diagonalizing the matrix, applying',&
' the functions to the individual eigenvalues and then transforming',&
' back. For example, "exp(A)" is computed by',&
'',&
' <X,D> = eig(A);',&
' for j = 1:n, D(j,j) = exp(D(j,j));',&
' X*D/X',&
'',&
' This is essentially method number 14 out of the 19 ''dubious''',&
' possibilities described in [8]. It is dubious because it doesn''t always',&
' work. The matrix of eigenvectors X can be arbitrarily badly conditioned',&
' and all accuracy lost in the computation of X*D/X. A warning message',&
' is printed if "rcond(X)" is very small, but this only catches the',&
' extreme cases. An example of a case not detected is when A has a double',&
' eigenvalue, but theoretically only one linearly independent eigenvector',&
' associated with it. The computed eigenvalues will be separated by',&
' something on the order of the square root of the roundoff level. This',&
' separation will be reflected in "rcond(X)" which will probably not',&
' be small enough to trigger the error message. The computed "exp(A)"',&
' will be accurate to only half precision. Better methods are known for',&
' computing "exp(A)", but they do not easily extend to the other five',&
' functions and would require a considerable amount of additional code.',&
'',&
' The expression A**p is evaluated by repeated multiplication if p is',&
' an integer greater than 1. Otherwise it is evaluated by',&
'',&
' <X,D> = eig(A);',&
' for j = 1:n, D(j,j) = exp(p*log(D(j,j)))',&
' X*D/X',&
'',&
' This suffers from the same potential loss of accuracy if X is',&
' badly conditioned. It was partly for this reason that the case p =',&
' 1 is included in the general case. Comparison of A**1 with A gives',&
' some idea of the loss of accuracy for other values of p and for the',&
' elementary functions.',&
'',&
' "rref", the reduced row echelon form, is of some interest in',&
' theoretical linear algebra, although it has little computational',&
' value. It is included in LALA for pedagogical reasons. The algorithm',&
' is essentially Gauss-Jordan elimination with detection of negligible',&
' columns applied to rectangular matrices.',&
'',&
' There are three separate places in LALA where the rank of a matrix',&
' is implicitly computed: in rref(A), in A\B for non-square A, and',&
' in the pseudoinverse pinv(A). Three different algorithms with three',&
' different criteria for negligibility are used and so it is possible',&
' that three different values could be produced for the same matrix. With',&
' rref(A), the rank of A is the number of nonzero rows. The elimination',&
' algorithm used for "rref" is the fastest of the three rank-determining',&
' algorithms, but it is the least sophisticated numerically and the',&
' least reliable. With A\B, the algorithm is essentially that used',&
' by example subroutine SQRST in chapter 9 of the LINPACK guide. With',&
' pinv(A), the algorithm is based on the singular value decomposition',&
' and is described in chapter 11 of the LINPACK guide. The SVD algorithm',&
' is the most time-consuming, but the most reliable and is therefore',&
' also used for rank(A).',&
'',&
' The uniformly distributed random numbers in "rand" are obtained from',&
' the machine-independent random number generator URAND described in',&
' [9]. It is possible to switch to normally distributed random numbers,',&
' which are obtained using a transformation also described in [9].',&
'',&
' The computation of',&
'',&
' 2 2',&
' sqrt(a + b )',&
'',&
' is required in many matrix algorithms, particularly those involving',&
' complex arithmetic. A new approach to carrying out this operation is',&
' described by Moler and Morrison [10]. It is a cubically convergent',&
' algorithm which starts with a and b , rather than with their squares,',&
' and thereby avoids destructive arithmetic underflows and overflows. In',&
' LALA, the algorithm is used for complex modulus, Euclidean vector',&
' norm, plane rotations, and the shift calculation in the eigenvalue',&
' and singular value iterations.',&
'',&
'================================================================================',&
'FLOP AND CHOP',&
'',&
' Detailed information about the amount of work involved in matrix',&
' calculations and the resulting accuracy is provided by "flop" and',&
' "chop". The basic unit of work is the "flop", or floating point',&
' operation. Roughly, one flop is one execution of a Fortran statement',&
' like',&
'',&
' S = S + X(I)*Y(I)',&
'',&
' or',&
'',&
' Y(I) = Y(I) + T*X(I)',&
'',&
' In other words, it consists of one floating point multiplication,',&
' together with one floating point addition and the associated',&
' indexing and storage reference operations.',&
'',&
' LALA will print the CPU time required for a particular',&
' statement when the statement is terminated by an extra comma. For',&
' example, the line',&
'',&
' n = 20; rand(n)*rand(n);,',&
'',&
' ends with an extra comma. Two 20 by 20 random matrices are generated',&
' and multiplied together. The result is assigned to "ans", but the',&
' semicolon suppresses its printing. The only output is',&
'',&
' 8800 flops',&
'',&
' This is n**3 + 2*n**2 flops, n**2 for each random matrix and n**3',&
' for the product.',&
'',&
' "flop" is a predefined vector with two components. "flop(1)" is',&
' the number of flops used by the most recently executed statement,',&
' except that statements with zero flops are ignored. For example,',&
' after executing the previous statement,',&
'',&
' flop(1)/n**3',&
'',&
' results in',&
'',&
' ans =',&
'',&
' 1.1000',&
'',&
' "flop(2)" is the cumulative total of all the flops used since',&
' the beginning of the LALA session. The statement',&
'',&
' flop = <0 0>',&
'',&
' resets the total.',&
'',&
' There are several difficulties associated with keeping a',&
' precise count of floating point operations.',&
'',&
' As the program generally uses complex values but only performs',&
' operations on the real matrices in many cases where all the imaginary',&
' values are zero it may not provide an accurate measure of the relative',&
' costs of real and complex arithmetic.',&
'',&
' The result of each floating point operation may also be "chopped"',&
' to simulate a computer with a shorter word length. The details',&
' of this chopping operation depend upon the format of the floating',&
' point word. Usually, the fraction in the floating point word can be',&
' regarded as consisting of several octal or hexadecimal digits. The',&
' least significant of these digits can be set to zero by a logical',&
' masking operation. Thus the statement',&
'',&
' chop(p)',&
'',&
' causes the p least significant octal or hexadecimal digits in',&
' the result of each floating point operation to be set to zero.',&
' For example, if the computer being used has an IBM 360 long floating',&
' point word with 14 hexadecimal digits in the fraction, then "chop(8)"',&
' results in simulation of a computer with only 6 hexadecimal digits',&
' in the fraction, i.e. a short floating point word. On a computer such',&
' as the CDC 6600 with 16 octal digits, "chop(8)" results in about the',&
' same accuracy because the remaining 8 octal digits represent the same',&
' number of bits as 6 hexadecimal digits.',&
'',&
' Some idea of the effect of "chop" on any particular system can',&
' be obtained by executing the following statements.',&
'',&
' long, t = 1/10',&
' long z, t = 1/10',&
' chop(8)',&
' long, t = 1/10',&
' long z, t = 1/10',&
'',&
' The following Fortran subprograms illustrate more details of',&
' "flop" and "chop". The first subprogram is a simplified example of a',&
' system-dependent function used within LALA itself. The common variable',&
' G_FLOP_COUNTER is essentially the first component of the variable',&
' FLOP. The common variable CHP is initially zero, but it is set to p',&
' by the statement "chop(p)". To shorten the DATA statement, we assume',&
' there are only 6 hexadecimal digits. We also assume an extension of',&
' Fortran that allows .AND. to be used as a binary operation between',&
' two real variables.',&
'',&
' REAL FUNCTION FLOP(X)',&
' REAL X',&
' INTEGER G_FLOP_COUNTER,CHP',&
' COMMON G_FLOP_COUNTER,CHP',&
' REAL MASK(5)',&
' DATA MASK/ZFFFFFFF0,ZFFFFFF00,ZFFFFF000,ZFFFF0000,ZFFF00000/',&
' G_FLOP_COUNTER = G_FLOP_COUNTER + 1',&
' IF (CHP .EQ. 0) FLOP = X',&
' IF (CHP .GE. 1 .AND. CHP .LE. 5) FLOP = X .AND. MASK(CHP)',&
' IF (CHP .GE. 6) FLOP = 0.0',&
' END REAL FUNCTION FLOP',&
'',&
' The following subroutine illustrates a typical use of the',&
' previous function within LALA. It is a simplified version of',&
' the Basic Linear Algebra Subprogram that adds a scalar multiple',&
' of one vector to another. We assume here that the vectors are',&
' stored with a memory increment of one.',&
'',&
' SUBROUTINE SAXPY(N,TR,TI,XR,XI,YR,YI)',&
' REAL TR,TI,XR(N),XI(N),YR(N),YI(N),FLOP',&
' IF (N .LE. 0) return',&
' IF (TR .EQ. 0.0 .AND. TI .EQ. 0.0) return',&
' DO I = 1, N',&
' YR(I) = FLOP(YR(I) + TR*XR(I) - TI*XI(I))',&
' YI(I) = YI(I) + TR*XI(I) + TI*XR(I)',&
' IF (YI(I) .NE. 0.0D0) YI(I) = FLOP(YI(I))',&
' enddo',&
' END SUBROUTINE SAXPY',&
'',&
' The saxpy operation is perhaps the most fundamental',&
' operation within LINPACK. It is used in the computation of the',&
' LU, the QR and the SVD factorizations, and in several other',&
' places. We see that adding a multiple of one vector with n',&
' components to another uses n flops if the vectors are real and',&
' between n and 2*n flops if the vectors have nonzero imaginary',&
' components.',&
'',&
' The permanent LALA variable "eps" is reset by the statement',&
' CHOP(p). Its new value is usually the smallest inverse power of',&
' two that satisfies the Fortran logical test',&
'',&
' FLOP(1.0+eps) .GT. 1.0',&
'',&
' However, if "eps" had been directly reset to a larger value, the',&
' old value is retained.',&
'',&
'================================================================================',&
'CENSUS EXAMPLE',&
'',&
' Our first extended example involves predicting the population of the',&
' United States in 1980 using extrapolation of various fits to the',&
' census data from 1900 through 1970. There are eight observations,',&
' so we begin with the LALA statement',&
'',&
' n = 8',&
'',&
' The values of the dependent variable, the population in millions,',&
' can be entered with',&
'',&
' y = < 75.995 91.972 105.711 123.203 ...',&
' 131.669 150.697 179.323 203.212>''',&
'',&
' In order to produce a reasonably scaled matrix, the independent',&
' variable, time, is transformed from the interval [1900,1970] to',&
' [-1.00,0.75]. This can be accomplished directly with',&
'',&
' t = -1.0:0.25:0.75',&
'',&
' or in a fancier, but perhaps clearer, way with',&
'',&
' t = 1900:10:1970; t = (t - 1940*ones(t))/40',&
'',&
' Either of these is equivalent to',&
'',&
' t = <-1 -.75 -.50 -.25 0 .25 .50 .75>',&
'',&
' The interpolating polynomial of degree n-1 involves an Vandermonde',&
' matrix of order n with elements that might be generated by',&
'',&
' for i = 1:n, for j = 1:n, a(i,j) = t(i)**(j-1);',&
'',&
' However, this results in an error caused by 0**0 when i = 5 and',&
' j = 1 . The preferable approach is',&
'',&
' A = ones(n,n);',&
' for i = 1:n, for j = 2:n, a(i,j) = t(i)*a(i,j-1);',&
'',&
' Now the statement',&
'',&
' cond(A)',&
'',&
' produces the output',&
'',&
' ans =',&
'',&
' 1.1819E+03',&
'',&
' which indicates that transformation of the time variable has resulted',&
' in a reasonably well conditioned matrix.',&
'',&
' The statement',&
'',&
' c = A\y',&
'',&
' results in',&
'',&
' C =',&
'',&
' 131.6690',&
' 41.0406',&
' 103.5396',&
' 262.4535',&
' -326.0658',&
' -662.0814',&
' 341.9022',&
' 533.6373',&
'',&
' These are the coefficients in the interpolating polynomial',&
'',&
' n-1',&
'',&
' c + c t + ... + c t',&
' 1 2 n',&
'',&
' Our transformation of the time variable has resulted in t = 1',&
' corresponding to the year 1980. Consequently, the extrapolated',&
' population is simply the sum of the coefficients. This can be',&
' computed by',&
'',&
' p = sum(c)',&
'',&
' The result is',&
'',&
' P =',&
'',&
' 426.0950',&
'',&
' which indicates a 1980 population of over 426 million. Clearly, using',&
' the seventh degree interpolating polynomial to extrapolate even a',&
' fairly short distance beyond the end of the data interval is not a',&
' good idea.',&
'',&
' The coefficients in least squares fits by polynomials of lower degree',&
' can be computed using fewer than n columns of the matrix.',&
'',&
' for k = 1:n, c = A(:,1:k)\y, p = sum(c)',&
'',&
' would produce the coefficients of these fits, as well as the',&
' resulting extrapolated population. If we do not want to print all the',&
' coefficients, we can simply generate a small table of populations',&
' predicted by polynomials of degrees zero through seven. We also',&
' compute the maximum deviation between the fitted and observed values.',&
'',&
' for k = 1:n, X = A(:,1:k); c = X\y; ...',&
' d(k) = k-1; p(k) = sum(c); e(k) = norm(X*c-y,''inf'');',&
' <d, p, e>',&
'',&
' The resulting output is',&
'',&
' 0 132.7227 70.4892',&
' 1 211.5101 9.8079',&
' 2 227.7744 5.0354',&
' 3 241.9574 3.8941',&
' 4 234.2814 4.0643',&
' 5 189.7310 2.5066',&
' 6 118.3025 1.6741',&
' 7 426.0950 0.0000',&
'',&
' The zeroth degree fit, 132.7 million, is the result of fitting a',&
' constant to the data and is simply the average. The results obtained',&
' with polynomials of degree one through four all appear reasonable. The',&
' maximum deviation of the degree four fit is slightly greater than the',&
' degree three, even though the sum of the squares of the deviations',&
' is less. The coefficients of the highest powers in the fits of degree',&
' five and six turn out to be negative and the predicted populations of',&
' less than 200 million are probably unrealistic. The hopefully absurd',&
' prediction of the interpolating polynomial concludes the table.',&
'',&
' We wish to emphasize that roundoff errors are not significant',&
' here. Nearly identical results would be obtained on other computers,',&
' or with other algorithms. The results simply indicate the difficulties',&
' associated with extrapolation of polynomial fits of even modest degree.',&
'',&
' A stabilized fit by a seventh degree polynomial can be obtained using',&
' the pseudoinverse, but it requires a fairly delicate choice of a',&
' tolerance. The statement',&
'',&
' s = svd(A)',&
'',&
' produces the singular values',&
'',&
' S =',&
'',&
' 3.4594',&
' 2.2121',&
' 1.0915',&
' 0.4879',&
' 0.1759',&
' 0.0617',&
' 0.0134',&
' 0.0029',&
'',&
' We see that the last three singular values are less than 0.1 ,',&
' consequently, A can be approximately by a matrix of rank five with an',&
' error less than 0.1 . The Moore-Penrose pseudoinverse of this rank',&
' five matrix is obtained from the singular value decomposition with',&
' the following statements',&
'',&
' c = pinv(A,0.1)*y, p = sum(c), e = norm(a*c-y,''inf'')',&
'',&
' The output is',&
'',&
' C =',&
'',&
' 134.7972',&
' 67.5055',&
' 23.5523',&
' 9.2834',&
' 3.0174',&
' 2.6503',&
' -2.8808',&
' 3.2467',&
'',&
' P =',&
'',&
' 241.1720',&
'',&
' E =',&
'',&
' 3.9469',&
'',&
' The resulting seventh degree polynomial has coefficients which are much',&
' smaller than those of the interpolating polynomial given earlier. The',&
' predicted population and the maximum deviation are reasonable. Any',&
' choice of the tolerance between the fifth and sixth singular values',&
' would produce the same results, but choices outside this range result',&
' in pseudoinverses of different rank and do not work as well.',&
'',&
' The one term exponential approximation',&
'',&
' y(t) = k exp(pt)',&
'',&
' can be transformed into a linear approximation by taking logarithms.',&
'',&
' log(y(t)) = log k + pt',&
'',&
' = c + c t',&
' 1 2',&
'',&
' The following segment makes use of the fact that a function of a',&
' vector is the function applied to the individual components.',&
'',&
' X = A(:,1:2);',&
' c = X\log(y)',&
' p = exp(sum(c))',&
' e = norm(exp(X*c)-y,''inf'')',&
'',&
' The resulting output is',&
'',&
' C =',&
'',&
' 4.9083',&
' 0.5407',&
'',&
' P =',&
'',&
' 232.5134',&
'',&
' E =',&
'',&
' 4.9141',&
'',&
' The predicted population and maximum deviation appear satisfactory and',&
' indicate that the exponential model is a reasonable one to consider.',&
'',&
' As a curiosity, we return to the degree six polynomial. Since the',&
' coefficient of the high order term is negative and the value of the',&
' polynomial at t = 1 is positive, it must have a root at some value',&
' of t greater than one. The statements',&
'',&
' X = A(:,1:7);',&
' c = X\y;',&
' c = c(7:-1:1); //reverse the order of the coefficients',&
' z = roots(c)',&
'',&
' produce',&
'',&
' Z =',&
'',&
' 1.1023- 0.0000*i',&
' 0.3021+ 0.7293*i',&
' -0.8790+ 0.6536*i',&
' -1.2939- 0.0000*i',&
' -0.8790- 0.6536*i',&
' 0.3021- 0.7293*i',&
'',&
' There is only one real, positive root. The corresponding time on the',&
' original scale is',&
'',&
' 1940 + 40*real(z(1))',&
'',&
' = 1984.091',&
'',&
' We conclude that the United States population should become zero',&
' early in February of 1984.',&
'',&
'================================================================================',&
'PARTIAL DIFFERENTIAL EQUATION EXAMPLE',&
'',&
' Our second extended example is a boundary value problem for Laplace''s',&
' equation. The underlying physical problem involves the conductivity',&
' of a medium with cylindrical inclusions and is considered by Keller',&
' and Sachs [7].',&
'',&
' Find a function u(x,y) satisfying Laplace''s equation',&
'',&
' u + u = 0',&
' xx yy',&
'',&
' The domain is a unit square with a quarter circle of radius rho removed',&
' from one corner. There are Neumann conditions on the top and bottom',&
' edges and Dirichlet conditions on the remainder of the boundary.',&
'',&
' u = 0',&
' n',&
'',&
' -------------',&
' | .',&
' | .',&
' | .',&
' | . u = 1',&
' | .',&
' | .',&
' | .',&
' u = 0 | |',&
' | |',&
' | |',&
' | | u = 1',&
' | |',&
' | |',&
' | |',&
' ------------------------',&
'',&
' u = 0',&
' n',&
'',&
' The effective conductivity of an medium is then given by the integral',&
' along the left edge,',&
'',&
' 1',&
' sigma = integral u (0,y) dy',&
' 0 n',&
'',&
' It is of interest to study the relation between the radius rho and',&
' the conductivity sigma. In particular, as rho approaches one, sigma',&
' becomes infinite.',&
'',&
' Keller and Sachs use a finite difference approximation. The following',&
' technique makes use of the fact that the equation is actually Laplace''s',&
' equation and leads to a much smaller matrix problem to solve.',&
'',&
' Consider an approximate solution of the form',&
'',&
' n 2j-1',&
' u = sum c r cos(2j-1)t',&
' j=1 j',&
'',&
' where r,t are polar coordinates (t is theta). The coefficients are',&
' to be determined. For any set of coefficients, this function already',&
' satisfies the differential equation because the basis functions are',&
' harmonic; it satisfies the normal derivative boundary condition on',&
' the bottom edge of the domain because we used cos t in preference to',&
' sin t ; and it satisfies the boundary condition on the left edge of',&
' the domain because we use only odd multiples of t .',&
'',&
' The computational task is to find coefficients so that the boundary',&
' conditions on the remaining edges are satisfied as well as possible. To',&
' accomplish this, pick m points (r,t) on the remaining edges. It is',&
' desirable to have m > n and in practice we usually choose m to be two',&
' or three times as large as n . Typical values of n are 10 or 20 and',&
' of m are 20 to 60. An m by n matrix A is generated. The i,j element',&
' is the j-th basis function, or its normal derivative, evaluated at',&
' the i-th boundary point. A right hand side with m components is also',&
' generated. In this example, the elements of the right hand side are',&
' either zero or one. The coefficients are then found by solving the',&
' overdetermined set of equations',&
'',&
' Ac = b',&
'',&
' in a least squares sense.',&
'',&
' Once the coefficients have been determined, the approximate solution',&
' is defined everywhere on the domain. It is then possible to compute the',&
' effective conductivity sigma . In fact, a very simple formula results,',&
'',&
' n j-1',&
' sigma = sum (-1) c',&
' j=1 j',&
'',&
' To use LALA for this problem, the following "program" is first stored',&
' in the local computer file system, say under the name "PDE".',&
'',&
' //Conductivity example.',&
' //Parameters ---',&
' rho //radius of cylindrical inclusion',&
' n //number of terms in solution',&
' m //number of boundary points',&
' //initialize operation counter',&
' flop = <0 0>;',&
' //initialize variables',&
' m1 = round(m/3); //number of points on each straight edge',&
' m2 = m - m1; //number of points with Dirichlet conditions',&
' pi = 4*atan(1);',&
' //generate points in Cartesian coordinates',&
' //right hand edge',&
' for i = 1:m1, x(i) = 1; y(i) = (1-rho)*(i-1)/(m1-1);',&
' //top edge',&
' for i = m2+1:m, x(i) = (1-rho)*(m-i)/(m-m2-1); y(i) = 1;',&
' //circular edge',&
' for i = m1+1:m2, t = pi/2*(i-m1)/(m2-m1+1); ...',&
' x(i) = 1-rho*sin(t); y(i) = 1-rho*cos(t);',&
' //convert to polar coordinates',&
' for i = 1:m-1, th(i) = atan(y(i)/x(i)); ...',&
' r(i) = sqrt(x(i)**2+y(i)**2);',&
' th(m) = pi/2; r(m) = 1;',&
' //generate matrix',&
' //Dirichlet conditions',&
' for i = 1:m2, for j = 1:n, k = 2*j-1; ...',&
' a(i,j) = r(i)**k*cos(k*th(i));',&
' //Neumann conditions',&
' for i = m2+1:m, for j = 1:n, k = 2*j-1; ...',&
' a(i,j) = k*r(i)**(k-1)*sin((k-1)*th(i));',&
' //generate right hand side',&
' for i = 1:m2, b(i) = 1;',&
' for i = m2+1:m, b(i) = 0;',&
' //solve for coefficients',&
' c = A\b',&
' //compute effective conductivity',&
' c(2:2:n) = -c(2:2:n);',&
' sigma = sum(c)',&
' //output total operation count',&
' ops = flop(2)',&
'',&
' The program can be used within LALA by setting the three parameters',&
' and then accessing the file. For example,',&
'',&
' rho = .9;',&
' n = 15;',&
' m = 30;',&
' exec(''PDE'')',&
'',&
' The resulting output is',&
'',&
' rho =',&
'',&
' .9000',&
'',&
' n =',&
'',&
' 15.',&
'',&
' m =',&
'',&
' 30.',&
'',&
' c =',&
'',&
' 2.2275',&
' -2.2724',&
' 1.1448',&
' 0.1455',&
' -0.1678',&
' -0.0005',&
' -0.3785',&
' 0.2299',&
' 0.3228',&
' -0.2242',&
' -0.1311',&
' 0.0924',&
' 0.0310',&
' -0.0154',&
' -0.0038',&
'',&
' sigm =',&
'',&
' 5.0895',&
'',&
' ops =',&
'',&
' 16204.',&
'',&
' A total of 16204 floating point operations were necessary to set up the',&
' matrix, solve for the coefficients and compute the conductivity. The',&
' operation count is roughly proportional to m*n**2. The results obtained',&
' for sigma as a function of rho by this approach are essentially the',&
' same as those obtained by the finite difference technique of Keller',&
' and Sachs, but the computational effort involved is much less.',&
'',&
'================================================================================',&
'EIGENVALUE SENSITIVITY EXAMPLE',&
'',&
' In this example, we construct a matrix whose eigenvalues are moderately',&
' sensitive to perturbations and then analyze that sensitivity. We',&
' begin with the statement',&
'',&
' B = <3 0 7; 0 2 0; 0 0 1>',&
'',&
' which produces',&
'',&
' B =',&
'',&
' 3. 0. 7.',&
' 0. 2. 0.',&
' 0. 0. 1.',&
'',&
' Obviously, the eigenvalues of B are 1, 2 and 3 . Moreover, since',&
' B is not symmetric, these eigenvalues are slightly sensitive to',&
' perturbation. (The value b(1,3) = 7 was chosen so that the elements',&
' of the matrix A below are less than 1000.)',&
'',&
' We now generate a similarity transformation to disguise the eigenvalues',&
' and make them more sensitive.',&
'',&
' L = <1 0 0; 2 1 0; -3 4 1>, M = L\L''',&
'',&
' L =',&
'',&
' 1. 0. 0.',&
' 2. 1. 0.',&
' -3. 4. 1.',&
'',&
' M =',&
'',&
' 1.0000 2.0000 -3.0000',&
' -2.0000 -3.0000 10.0000',&
' 11.0000 18.0000 -48.0000',&
'',&
' The matrix M has determinant equal to 1 and is moderately badly',&
' conditioned. The similarity transformation is',&
'',&
' A = M*B/M',&
'',&
' A =',&
'',&
' -64.0000 82.0000 21.0000',&
' 144.0000 -178.0000 -46.0000',&
' -771.0000 962.0000 248.0000',&
'',&
' Because det(M) = 1 , the elements of A would be exact integers',&
' if there were no roundoff. So,',&
'',&
' A = round(A)',&
'',&
' A =',&
'',&
' -64. 82. 21.',&
' 144. -178. -46.',&
' -771. 962. 248.',&
'',&
' This, then, is our test matrix. We can now forget how it',&
' was generated and analyze its eigenvalues.',&
'',&
' <X,D> = eig(A)',&
'',&
' D =',&
'',&
' 3.0000 0.0000 0.0000',&
' 0.0000 1.0000 0.0000',&
' 0.0000 0.0000 2.0000',&
'',&
' X =',&
'',&
' -.0891 3.4903 41.8091',&
' .1782 -9.1284 -62.7136',&
' -.9800 46.4473 376.2818',&
'',&
' Since A is similar to B, its eigenvalues are also 1, 2 and 3. They',&
' happen to be computed in another order by the EISPACK subroutines. The',&
' fact that the columns of X, which are the eigenvectors, are so far',&
' from being orthonormal is our first indication that the eigenvalues',&
' are sensitive. To see this sensitivity, we display more figures of',&
' the computed eigenvalues.',&
'',&
' long, diag(D)',&
'',&
' ans =',&
'',&
' 2.999999999973599',&
' 1.000000000015625',&
' 2.000000000011505',&
'',&
' We see that, on this computer, the last five significant figures are',&
' contaminated by roundoff error. A somewhat superficial explanation',&
' of this is provided by',&
'',&
' short, cond(X)',&
'',&
' ans =',&
'',&
' 3.2216e+05',&
'',&
' The condition number of X gives an upper bound for the relative',&
' error in the computed eigenvalues. However, this condition',&
' number is affected by scaling.',&
'',&
' X = X/diag(X(3,:)), cond(X)',&
'',&
' X =',&
'',&
' .0909 .0751 .1111',&
' -.1818 -.1965 -.1667',&
' 1.0000 1.0000 1.0000',&
'',&
' ans =',&
'',&
' 1.7692e+03',&
'',&
' Rescaling the eigenvectors so that their last components are all',&
' equal to one has two consequences. The condition of X is decreased',&
' by over two orders of magnitude. (This is about the minimum condition',&
' that can be obtained by such diagonal scaling.) Moreover, it is now',&
' apparent that the three eigenvectors are nearly parallel.',&
'',&
' More detailed information on the sensitivity of the individual',&
' eigenvalues involves the left eigenvectors.',&
'',&
' Y = inv(X''), Y''*A*X',&
'',&
' Y =',&
'',&
' -511.5000 259.5000 252.0000',&
' 616.0000 -346.0000 -270.0000',&
' 159.5000 -86.5000 -72.0000',&
'',&
' ans =',&
'',&
' 3.0000 .0000 .0000',&
' .0000 1.0000 .0000',&
' .0000 .0000 2.0000',&
'',&
' We are now in a position to compute the sensitivities of the individual',&
' eigenvalues.',&
'',&
' for j = 1:3, c(j) = norm(Y(:,j))*norm(X(:,j)); end, C',&
'',&
' C =',&
'',&
' 833.1092',&
' 450.7228',&
' 383.7564',&
'',&
' These three numbers are the reciprocals of the cosines of the',&
' angles between the left and right eigenvectors. It can be shown that',&
' perturbation of the elements of A can result in a perturbation of',&
' the j-th eigenvalue which is c(j) times as large. In this example,',&
' the first eigenvalue has the largest sensitivity.',&
'',&
' We now proceed to show that A is close to a matrix with a double',&
' eigenvalue. The direction of the required perturbation is given by',&
'',&
' E = -1.e-6*Y(:,1)*X(:,1)''',&
'',&
' E =',&
'',&
' 1.0e-03 *',&
'',&
' .0465 -.0930 .5115',&
' -.0560 .1120 -.6160',&
' -.0145 .0290 -.1595',&
'',&
' With some trial and error which we do not show, we bracket the',&
' point where two eigenvalues of a perturbed A coalesce and then',&
' become complex.',&
'',&
' eig(A + .4*E), eig(A + .5*E)',&
'',&
' ans =',&
'',&
' 1.1500',&
' 2.5996',&
' 2.2504',&
'',&
' ans =',&
'',&
' 2.4067 + .1753*i',&
' 2.4067 - .1753*i',&
' 1.1866 + 0.0000*i',&
'',&
' Now, a bisecting search, driven by the imaginary part of one of',&
' the eigenvalues, finds the point where two eigenvalues are nearly',&
' equal.',&
'',&
' r = .4; s = .5;',&
'',&
' while s-r > 1.e-14, t = (r+s)/2; d = eig(A+t*E); ...',&
' if imag(d(1))=0, r = t; else, s = t;',&
'',&
' long, T',&
'',&
' T =',&
'',&
' .450380734134507',&
'',&
' Finally, we display the perturbed matrix, which is obviously close',&
' to the original, and its pair of nearly equal eigenvalues. (We have',&
' dropped a few digits from the long output.)',&
'',&
' A+t*E, eig(A+t*E)',&
'',&
' A',&
'',&
' -63.999979057 81.999958114 21.000230369',&
'<PHONE_NUMBER>78 -177.999949557 -46.000277434',&
' -771.000006530<PHONE_NUMBER>61<PHONE_NUMBER>64',&
'',&
' ans =',&
'',&
' 2.415741150',&
' 2.415740621',&
' 1.168517777',&
'',&
' The first two eigenvectors of A + t*E are almost indistinguishable',&
' indicating that the perturbed matrix is almost defective.',&
'',&
' <X,D> = eig(A+t*E); X = X/diag(X(3,:))',&
'',&
' X =',&
'',&
' .096019578 .096019586 .071608466',&
' -.178329614 -.178329608 -.199190520',&
' 1.000000000 1.000000000 1.000000000',&
'',&
' short, cond(X)',&
'',&
' ans =',&
'',&
' 3.3997e+09',&
'',&
'================================================================================',&
'COMMUNICATING WITH OTHER PROGRAMS',&
'',&
' There are four different ways LALA can be used in',&
' conjunction with other programs:',&
'',&
' -- user() - a user-supplied subroutine',&
' -- exec() - reading commands from a file',&
' -- save() and load() -- reading specially formatted data files.',&
' -- lala() - call the interpreter with a CHARACTER array of',&
' commands or interactively.',&
'',&
' Let us illustrate each of these by equivalents of the following',&
' simple example.',&
'',&
' You can start the lala(1) program up and simply enter:',&
'',&
' n = 6',&
' for i = 1:n, for j = 1:n, a(i,j) = abs(i-j);',&
' a',&
' x = inv(a)',&
'',&
' An example user routine could be introduced into LALA that',&
' does the same thing as the "for" statement by compiling and',&
' linking the following subroutine into the calling program.',&
'',&
' program demo_user',&
' implicit none',&
' use M_matrix',&
' call set_usersub(lala_user)',&
' call lala()',&
' subroutine lala_user(a,m,n,s,t)',&
' implicit none',&
' doubleprecision a(:),s,t',&
' integer m,n',&
' n = int(a(1))',&
' m = n',&
' do j = 1, n',&
' do i = 1, n',&
' k = i + (j-1)*m',&
' a(k) = iabs(i-j)',&
' enddo',&
' enddo',&
' end subroutine lala_user',&
' end program demo_user',&
'',&
' A user-defined function can then be registered with the program',&
' with',&
'',&
' call set_usersub(SUBROUTINE_NAME)',&
'',&
' Note the routine must be defined with an explicit interface',&
' available in the calling unit.',&
'',&
' Then the LALA statements',&
'',&
' n = 6',&
' a = user(n)',&
' x = inv(a)',&
'',&
' do the job.',&
'',&
' The example procedure could be called by storing the following',&
' text in a file named, say, EXAMPLE.',&
'',&
' for i = 1:n, for j = 1:n, a(i,j) = abs(i-j);',&
'',&
' Then the LALA statements',&
'',&
' n = 6',&
' exec(''EXAMPLE'',0)',&
' x = inv(a)',&
'',&
' have the desired effect. The 0 as the optional second parameter',&
' of exec indicates that the text in the file should not be printed',&
' on the terminal.',&
'',&
' The matrices A and X could also be stored in files. Two',&
' separate main programs would be involved. The first is:',&
'',&
' program maina',&
' doubleprecision a(10,10)',&
' n = 6',&
' do j = 1, n',&
' do i = 1, n',&
' a(i,j) = iabs(i-j)',&
' enddo',&
' enddo',&
' OPEN(UNIT=1,FILE=''A'')',&
' write(1,''(a32,2i4)'') ''a'', n,n',&
' do j = 1, n',&
' write(1,102) (a(i,j),i=1,n)',&
' enddo',&
' 102 format(4z18)',&
' end program maina',&
'',&
' The OPEN statement may take different forms on different systems.',&
' It attaches Fortran logical unit number 1 to the file named A.',&
'',&
' The FORMAT number 102 may also be system dependent. This',&
' particular one is appropriate for hexadecimal computers with an 8',&
' byte double precision floating point word. Check, or modify,',&
' LALA subroutine SAVLOD.',&
'',&
' After this program is executed, enter LALA and give the',&
' following statements:',&
'',&
' load(''A'')',&
' X = inv(a)',&
' save(''X'',X)',&
'',&
' If all goes according to plan, this will read the matrix "a" from',&
' the file A, invert it, store the inverse in X and then write the',&
' matrix X on the file X. The following program can then access X.',&
'',&
' program mainx',&
' doubleprecision x(10,10)',&
' open(unit=1,file=''x'')',&
' rewind 1',&
' read (1, ''(a32,2i4)'') id,m,n',&
' do j = 1, n',&
' read(1,''(4z18)'') (x(i,j),i=1,m)',&
' enddo',&
' ...',&
' ...',&
'',&
'',&
' The most elaborate mechanism involves using LALA as a subroutine',&
' within another program. Communication with the LALA stack is',&
' accomplished using subroutine lala().',&
' The preamble of MATZ is:',&
'',&
' SUBROUTINE MATZ(A,LDA,M,N,ID,JOB,IERR)',&
' INTEGER LDA,M,N,JOB,IERR',&
' character(len=*) :: id',&
' DOUBLEPRECISION A(LDA,N)',&
'',&
' ! ACCESS LALA VARIABLE STACK',&
' !  LEADING DIMENSION LDA.',&
' ! ID IS THE NAME OF A. ID IS UP TO FOUR CHARACTERS.',&
' ! JOB = 0 GET REAL A FROM LALA,',&
' ! = 1 PUT REAL A INTO LALA,',&
' ! = 10 GET IMAG PART OF A FROM LALA,',&
' ! = 11 PUT IMAG PART OF A INTO LALA.',&
' ! RETURN WITH NONZERO IERR AFTER LALA ERROR MESSAGE.',&
' !',&
' ! USES LALA ROUTINES STACKG, STACKP AND ERROR',&
'',&
' The preamble of subroutine LALA is:',&
'',&
' SUBROUTINE LALA(INIT)',&
' ! INIT = 0 FOR FIRST ENTRY, NONZERO FOR SUBSEQUENT ENTRIES',&
'',&
' To do our example, write the following program:',&
'',&
' DOUBLEPRECISION A(10,10),X(10,10)',&
' DATA LDA/10/',&
' call M_88(0,'''')',&
' N = 6',&
' DO J = 1, N',&
' DO I = 1, N',&
' A(I,J) = IABS(I-J)',&
' enddo',&
' enddo',&
' call MATZ(A,LDA,N,N,''A'',1,IERR)',&
' IF (IERR .NE. 0) GO TO ...',&
' call LALA(1,'''')',&
' call MATZ(X,LDA,N,N,''X'',0,IERR)',&
' IF (IERR .NE. 0) GO TO ...',&
' ...',&
' ...',&
'',&
' When this program is executed, the call to LALA(0) produces the',&
' LALA greeting, then waits for input. The command',&
'',&
' quit',&
'',&
' sends control back to our example program. The matrix A is',&
' generated by the program and sent to the stack by the first call',&
' to MATZ. The call to LALA(1) produces the LALA(1) prompt. Then',&
' the statements',&
'',&
' X = inv(A)',&
' quit',&
'',&
' will invert our matrix, put the result on the stack and go back',&
' to our program. The second call to MATZ will retrieve X .',&
'',&
' By the way, this matrix X is interesting. Take a look at',&
' round(2*(n-1)*X).',&
'',&
'================================================================================',&
'ACKNOWLEDGEMENT',&
'',&
' LALA was inspired by the MATLAB subroutine. Most of the work on',&
' MATLAB was carried out at the University of New Mexico, where it was',&
' being supported by the National Science Foundation. Additional work',&
' has been done during visits to Stanford Linear Accelerator Center,',&
' Argonne National Laboratory and Los Alamos Scientific Laboratory,',&
' where support has been provided by NSF and the Department of Energy.',&
'',&
'================================================================================',&
'REFERENCES FOR THE MATLAB ROUTINE',&
'',&
' [1] J. J. Dongarra, J. R. Bunch, C. B. Moler and G. W. Stewart,',&
' LINPACK Users'' Guide, Society for Industrial and Applied',&
' Mathematics, Philadelphia, 1979.',&
'',&
' [2] B. T. Smith, J. M. Boyle, J. J. Dongarra, B. S. Garbow, Y.',&
' Ikebe, V. C. Klema, C. B. Moler, Matrix Eigensystem Routines',&
' -- EISPACK Guide, Lecture Notes in Computer Science, volume',&
' 6, second edition, Springer-Verlag, 1976.',&
'',&
' [3] B. S. Garbow, J. M. Boyle, J. J. Dongarra, C. B. Moler,',&
' Matrix Eigensystem Routines -- EISPACK Guide Extension,',&
' Lecture Notes in Computer Science, volume 51, Springer-',&
' Verlag, 1977.',&
'',&
' [4] S. Cohen and S. Piper, SPEAKEASY III Reference Manual,',&
' Speakeasy Computing Corp., Chicago, Ill., 1979.',&
'',&
' [5] J. H. Wilkinson and C. Reinsch, Handbook for Automatic',&
' Computation, volume II, Linear Algebra, Springer-Verlag,',&
' 1971.',&
'',&
' [6] Niklaus Wirth, Algorithms + Data Structures = Programs,',&
' Prentice-Hall, 1976.',&
'',&
' [7] H. B. Keller and D. Sachs, "Calculations of the Conductivity',&
' of a Medium Containing Cylindrical Inclusions", J. Applied',&
' Physics 35, 537-538, 1964.',&
'',&
' [8] C. B. Moler and C. F. Van Loan, Nineteen Dubious Ways to',&
' Compute the Exponential of a Matrix, SIAM Review 20, 801-',&
' 836, 1979.',&
'',&
' [9] G. E. Forsythe, M. A. Malcolm and C. B. Moler, Computer',&
' Methods for Mathematical Computations, Prentice-Hall, 1977.',&
'',&
' [10] C. B. Moler and D. R. Morrison, "Replacing square roots by',&
' Pythagorean sums", University of New Mexico, Computer',&
' Science Department, technical report, submitted for',&
' publication, 1980.',&
'',&
'================================================================================',&
'SUMMARY A list of basic (case-sensitive) section and topic names',&
' .______________._________________________________________________________.',&
' |SYNTAX | [ ] < > ( ) = . , ! ; \ / '''' + - * : semi ? |',&
' |______________._________________________________________________________|',&
' |VARIABLES | ans clear who |',&
' |______________._________________________________________________________|',&
' |BASIC | atan cos exp log sin sqrt |',&
' |______________._________________________________________________________|',&
' |HIGH | abs base chol chop cond conjg det diag |',&
' | | eig eye hess invh imag inv kron lu |',&
' | | magic norm ones orth pinv poly prod qr |',&
' | | rand rank rcond rat real rref roots round |',&
' | | schur shape sum svd tril triu user zeros |',&
' |______________._________________________________________________________|',&
' |FLOW control | else end if for while exit quit |',&
' |______________._________________________________________________________|',&
' |FILE access | exec load print save delete |',&
' |______________._________________________________________________________|',&
' |OUTPUT options| lines long short diary display plot |',&
' |______________._________________________________________________________|',&
' |ENVIRONMENT | getenv |',&
' |______________._________________________________________________________|',&
' |DOCUMENTATION | help fhelp NEWS |',&
' |______________._________________________________________________________|',&
' |MISCELLANEOUS | eps debug flops sh MACROS EDIT CHARS |',&
' |______________._________________________________________________________|',&
'================================================================================',&
'SAMPLE',&
' Here are a few sample statements:',&
'',&
' A = <1 2; 3 4>',&
' b = <5 6>''',&
' x = A\b',&
' <V,D> = eig(A), norm(A-V*D/V)',&
' help \ , help eig',&
' exec(''demo'',7)',&
'',&
' For more information, generate the LALA Users'' Guide',&
' using',&
'',&
' help manual',&
' w help.txt',&
' q',&
'================================================================================',&
'DOCUMENTATION',&
'fhelp topic|SECTION_NAME',&
'',&
' "fhelp" is identical in usage to "help" except that it searches a',&
' collection of descriptions of Fortran intrinsics.',&
'',&
' fhelp verify',&
' fhelp pack',&
'',&
' See "help"',&
'',&
'help topic|SECTION_NAME',&
'',&
' "help" gives assistance. It is equivalent to "help SUMMARY"',&
' by default.',&
'',&
' o "help" with no options lists common topic and section names.',&
' o The special topic "topics" shows all topic lines.',&
' o The special topic "manual" displays all the help text.',&
' o The special topic "search" shows lines from the manual',&
' containing the subsequent string',&
'',&
' Enter "h" at the "continue ..." prompt for additional options.',&
'',&
' For example:',&
'',&
' help // a list of common topics and section names',&
' help topics // a list of topics including the first line of',&
' // the topic.',&
' help abs // produces help on the function "abs".',&
' help FLOW // the entire section on flow control is displayed.',&
' help manual // show all the help text',&
' help help // obviously prints this message.',&
' help search factor // show all lines containing "factor".',&
'',&
' Alternatively, To place all the documenation in a file, use',&
' "help manual" and enter "w help.txt" at the "continue .." prompt.',&
'NEWS',&
' LALA is intended to be used primarily by families of FORTRAN',&
' programs that wish to add a consistent interactive "calculator"',&
' mode for interactively inspecting and modifying data.',&
'',&
' May, 1981.',&
'',&
' This is a port of the Argonne National Lab. FORTRAN 77 MATLAB',&
' routine circa 1981.',&
'',&
' Mar, 1990.',&
'',&
' Input lines can now be recalled and edited. A "??" on a line by',&
' itself calls the command history mode. Enter "?" after entering',&
' the mode for details.',&
'',&
' Apr, 2021.',&
'',&
' Rewritten but largely true to the original documentation.',&
'',&
'what does nothing for now',&
'',&
'sh Starts the command shell interactively, using the command defined by',&
' the environment variable SHELL. Note that in addition any line',&
' starting with an exclamation (!) is passed to the system for',&
' execution.',&
'================================================================================',&
'SYNTAX',&
'[ See "<"',&
'] See "<"',&
'> See "<" . Also see MACROS.',&
'< < > or [ ] are brackets used in forming vectors and matrices.',&
' "<6.9 9.64 sqrt(-1)>" is a vector with three elements separated by',&
' blanks. "[1+I 2-I 3]" and "[1 +I 2 -I 3]" are not the same. The',&
' first has three elements, the second has five. <11 12 13; 21 22',&
' 23> is a 2 by 3 matrix. The semicolon ends the first row.',&
'',&
' Vectors and matrices can be used inside < > brackets. <A B; C>',&
' is allowed if the number of rows of A equals the number of rows',&
' of B and the number of columns of A plus the number of columns of',&
' B equals the number of columns of C. This rule generalizes in a',&
' hopefully obvious way to allow fairly complicated constructions.',&
'',&
' A = < > stores an empty matrix in A, thereby removing it from the',&
' list of current variables.',&
'',&
' For the use of < and > on the left of the = in multiple assignment',&
' statements, See "lu", "eig", "svd" and so on.',&
'',&
' In "while" and "if" clauses, "<>" means less than or greater than,',&
' i.e. not equal, "<" means less than, ">" means greater than,',&
' "<=" means less than or equal, ">=" means greater than or equal.',&
'',&
' For the use of ">" and "<" to delineate macros, see MACROS.',&
'',&
'{ see "(".',&
'} see "(".',&
') See "(" .',&
'( ( ) or { } are used to indicate precedence in arithmetic expressions',&
' and to enclose arguments of functions in the usual way. They are',&
' also used to enclose subscripts of vectors and matrices in a manner',&
' somewhat more general than the usual way. If X and V are vectors,',&
' then X(V) is <X(V(1)), X(V(2)), ..., X(V(N))>. The components of V',&
' are rounded to nearest integers and used as subscripts. An error',&
' occurs if any such subscript is less than 1 or greater than the',&
' dimension of X. Some examples:',&
'',&
' X(3) is the third element of X .',&
' X([1 2 3]) is the first three elements of X. So is',&
' X([sqrt(2), sqrt(3), 4*atan(1)]) .',&
' If X has N components, X(N:-1:1) reverses them.',&
'',&
' The same indirect subscripting is used in matrices. If V has',&
' M components and W has N components, then A(V,W) is the M by N',&
' matrix formed from the elements of A whose subscripts are the',&
' elements of V and W. For example... A(<1,5>,:) = A(<5,1>,:)',&
' interchanges rows 1 and 5 of A.',&
'',&
'= Used in assignment statements and to mean equality in "while"',&
' and "if" clauses.',&
'',&
'. Decimal point. 314/100, 3.14 and .314E1 are all the',&
' same.',&
'',&
' Element-by-element multiplicative operations are obtained',&
' using .* , ./ , or .\ . For example, C = A ./ B is the',&
' matrix with elements c(i,j) = a(i,j)/b(i,j) .',&
'',&
' Kronecker tensor products and quotients are obtained with',&
' .*. , ./. and .\. . See "kron".',&
'',&
' Two or more points at the end of the line indicate',&
' continuation. The total line length limit is 1024',&
' characters.',&
'',&
', Used to separate matrix subscripts and function arguments.',&
' Used at the end of "for", "while" and "if" clauses. Used to',&
' separate statements in multi-statement lines. In this',&
' situation, it may be replaced by a semicolon to suppress',&
' printing.',&
'',&
'! If an exclamation is the first character of a line the',&
' rest of the line is passed to the system to be executed.',&
'',&
' Examples:',&
'',&
' // enter command history mode and change all occurrences of',&
' // "abc" to "123" on the last command entered.',&
' !!c/abc/123',&
'',&
' // pass command to system',&
' !ls -ltrasd',&
'',&
' see "EDIT"',&
'',&
'; Used inside brackets to end rows.',&
'',&
' Used after an expression or statement to suppress printing.',&
' See "semi".',&
'',&
'\ Backslash or matrix left division. A\B is roughly the',&
' same as "inv(A)*B", except it is computed in a different',&
' way. If A is an N by N matrix and B is a column vector',&
' with N components, or a matrix with several such columns,',&
' then X = A\B is the solution to the equation A*X = B',&
' computed by Gaussian elimination. A warning message is',&
' printed if A is badly scaled or nearly singular.',&
' A\eye produces the inverse of A .',&
'',&
' If A is an M by N matrix with M < or > N and B is a',&
' column vector with M components, or a matrix with several',&
' such columns, then X = A\B is the solution in the least',&
' squares sense to the under- or overdetermined system of',&
' equations A*X = B. The effective rank, K, of A is',&
' determined from the QR decomposition with pivoting. A',&
' solution X is computed which has at most K nonzero',&
' components per column. If K < N this will usually not be',&
' the same solution as pinv(A)*B .',&
' A\eye produces a generalized inverse of A.',&
'',&
' If A and B have the same dimensions, then A .\ B has',&
' elements a(i,j)\b(i,j) .',&
'',&
' Also, see "edit".',&
'',&
'/ Slash or matrix right division. B/A is roughly the same',&
' as B*inv(A) . More precisely, B/A = (A''\B'')'' . See \ .',&
'',&
' IF A and B have the same dimensions, then A ./ B has',&
' elements a(i,j)/b(i,j) .',&
'',&
' Two or more slashes together on a line indicate a logical end of',&
' line. Any following text is ignored.',&
'',&
''' Transpose. X'' is the complex conjugate transpose of X .',&
'',&
' A quote is also use to delmit text. ''ANY TEXT'' is a vector whose',&
' components are the LALA internal codes for the characters. A',&
' quote within the text is indicated by two quotes. See "display"',&
' and "FILE" .',&
'',&
'+ Addition. X + Y . X and Y must have the same dimensions.',&
'',&
'- Subtraction. X - Y . X and Y must have the same',&
' dimensions.',&
'',&
'* Matrix multiplication, X*Y . Any scalar (1 by 1 matrix)',&
' may multiply anything. Otherwise, the number of columns of',&
' X must equal the number of rows of Y .',&
'',&
' Element-by-element multiplication is obtained with X .* Y .',&
'',&
' The Kronecker tensor product is denoted by X .*. Y .',&
'',&
' Powers. X**p is X to the p power. p must be a',&
' scalar. If X is a matrix, see "HIGH" .',&
'',&
': Colon. Used in subscripts, "for" iterations and possibly',&
' elsewhere.',&
'',&
' j:k is the same as <j, j+1, ..., k>',&
' is empty if j > k .',&
' j:i:k is the same as [j, j+i,j+2*i, ..., k]',&
' (Fortran DO loop users beware of the unusual order!)',&
'',&
' j:i:k is the same as <j, j+i, j+2i, ..., k>',&
' j:i:k is empty if i > 0 and j > k or if i < 0 and j < k .',&
'',&
' The colon notation can be used to pick out selected rows,',&
' columns and elements of vectors and matrices.',&
'',&
' A(:) is all the elements of A, regarded as a single column.',&
' However, used on the left side of an assignment, A(:)',&
' fills A, but preserves its shape.',&
' A(:,j) is the j-th column of A',&
' A(j:k) is A(j), A(j+1), ... , A(k)',&
' A(:,j:k) is A(:,j), A(:,j+1), ... ,A(:,k) and so on.',&
' A(:,:) is the same as A.',&
'',&
' For the use of the colon in the "for" statement, See "for" .',&
'',&
'semi "semi" toggles the action of semicolons at the end of lines.',&
' It will make semicolons cause rather than suppress printing.',&
' A second "semi" restores the initial interpretation.',&
'================================================================================',&
'VARIABLES',&
'',&
'ans Variable created automatically when expressions are not',&
' assigned to anything else.',&
'',&
'clear Erases all variables, except "eps", "flop", "eye" and "rand".',&
' X = <> erases only variable X . So does "clear X".',&
'',&
'who Lists current variables.',&
'================================================================================',&
'MACROS',&
'',&
' The macro facility involves text and inward pointing angle',&
' brackets. If "STRING" is the source text for any LALA',&
' expression or statement, then',&
'',&
' t = ''STRING'';',&
' encodes the text as a vector of integers and stores that',&
' vector in t. "display(t)" will print the text and',&
'',&
' >t<',&
' causes the text to be interpreted, either as a statement or',&
' as a factor in an expression. For example',&
'',&
' t = ''1/(i+j-1)'';',&
' display(t)',&
' for i = 1:n, for j = 1:n, a(i,j) = >t<;',&
'',&
' generates the Hilbert matrix of order n.',&
' Another example showing indexed text,',&
'',&
' S = <''x = 3 ''',&
' ''y = 4 ''',&
' ''z = sqrt(x*x+y*y)''>',&
' for k = 1:3, >S(k,:)<',&
'',&
' It is necessary that the strings making up the "rows" of',&
' the "matrix" S have the same lengths.',&
'',&
'================================================================================',&
'BASIC FUNCTIONS',&
'',&
' For matrix arguments X , the functions "sin", "cos", "atan",',&
' "sqrt", "log", "exp" and X**p are computed using eigenvalues D',&
' and eigenvectors V . If <V,D> = eig(X) then f(X) = V*f(D)/V . This',&
' method may give inaccurate results if V is badly conditioned. Some',&
' idea of the accuracy can be obtained by comparing X**1 with X .',&
' For vector arguments, the function is applied to each component.',&
'',&
'atan atan(X) is the arctangent of X . See "BASIC" .',&
'',&
'cos cos(X) is the cosine of X . See "BASIC" .',&
'',&
'exp exp(X) is the exponential of X , e to the X . See "BASIC".',&
'',&
'log log(X) is the natural logarithm of X.',&
'',&
' Complex results are produced if X is not positive, or has',&
' nonpositive eigenvalues.',&
'',&
' See "BASIC".',&
'',&
'sin sin(X) is the sine of X. See "BASIC".',&
'',&
'sqrt sqrt(X) is the square root of X. See "BASIC". Complex',&
' results are produced if X is not positive, or has',&
' nonpositive eigenvalues.',&
'================================================================================',&
'HIGH LEVEL FUNCTIONS',&
'',&
'abs abs(X) is the absolute value, or complex modulus,',&
' of the elements of X .',&
'',&
'base base(X,B) is a vector containing the base B representation',&
' of X. This is often used in conjunction with "display".',&
' "display(X,B)" is the same as "display(base(X,B))". For example,',&
' "display(4*atan(1),16)" prints the hexadecimal representation of pi.',&
'',&
'chol Cholesky factorization. "chol(X)" uses only the diagonal',&
' and upper triangle of X. The lower triangular is assumed to be',&
' the (complex conjugate) transpose of the upper. If X is positive',&
' definite, then "R = chol(X)" produces an upper triangular R so',&
' that R''*R = X . If X is not positive definite, an error message',&
' is printed.',&
'',&
'chop Truncate arithmetic. "chop(P)" causes P places to be chopped',&
' off after each arithmetic operation in subsequent computations. This',&
' means P hexadecimal digits on some computers and P octal digits',&
' on others. "chop(0)" restores full precision.',&
'',&
'cond Condition number in 2-norm. "cond(X)" is the ratio of the',&
' largest singular value of X to the smallest.',&
'',&
'conjg "conjg(X)" is the complex conjugate of X .',&
'',&
'det "det(X)" is the determinant of the square matrix X .',&
'',&
'diag If V is a row or column vector with N components,',&
' "diag(V,K)" is a square matrix of order "N+abs(K)" with the',&
' elements of V on the K-th diagonal. K = 0 is the main diagonal,',&
' K > 0 is above the main diagonal and K < 0 is below the main',&
' diagonal. "diag(V)" simply puts V on the main diagonal. eg.',&
'',&
' diag(-M:M) + diag(ones(2*M,1),1) + diag(ones(2*M,1),-1)',&
'',&
' produces a tridiagonal matrix of order 2*M+1 .',&
'',&
' If X is a matrix, "diag(X,K)" is a column vector formed from the',&
' elements of the K-th diagonal of X. "diag(X)" is the main diagonal',&
' of X. "diag(diag(X))" is a diagonal matrix .',&
'',&
'eig Eigenvalues and eigenvectors.',&
' "eig(X)" is a vector containing the eigenvalues of a square',&
' matrix X.',&
' "<V,D> = eig(X)" produces a diagonal matrix D of',&
' eigenvalues and a full matrix V whose columns are the',&
' corresponding eigenvectors so that X*V = V*D .',&
'',&
'eye Identity matrix. "eye(N)" is the N by N identity matrix.',&
' "eye(M,N)" is an M by N matrix with 1''s on the diagonal and',&
' zeros elsewhere. "eye(A)" is the same size as A. "eye"',&
' with no arguments is an identity matrix of whatever order',&
' is appropriate in the context. For example "A + 3*eye"',&
' adds 3 to each diagonal element of A.',&
'',&
'hess Hessenberg form. The Hessenberg form of a matrix is zero',&
' below the first subdiagonal. If the matrix is symmetric or',&
' Hermitian, the form is tridiagonal. <P,H> = "hess(A)" produces a',&
' unitary matrix P and a Hessenberg matrix H so that A = P*H*P''. By',&
' itself, "hess(A)" returns H.',&
'',&
'invh Inverse Hilbert matrix. "invh(N)" is the inverse of a N_by_N',&
' Hilbert matrix (which is a famous example of a badly conditioned',&
' matrix). The result is exact for N less than about 15, depending',&
' upon the computer.',&
'',&
' for i = 1:N, for j = 1:N, A(i,j) = 1/(i+j-1);',&
'',&
' generates the NxN Hilbert matrix.',&
'',&
' "invh" has an alias of "inverse_hilbert" and "invhilb".',&
'',&
'aimag see "imag"',&
'imag "imag(X)" is the imaginary part of X .',&
'',&
'inv "inv(X)" is the inverse of the square matrix X . A warning',&
' message is printed if X is badly scaled or nearly',&
' singular.',&
'',&
'kron "kron(X,Y)" is the Kronecker tensor product of X and Y. It',&
' is also denoted by X .*. Y . The result is a large matrix',&
' formed by taking all possible products between the elements',&
' of X and those of Y . For example, if X is 2 by 3, then',&
' X .*. Y is',&
'',&
' < x(1,1)*Y x(1,2)*Y x(1,3)*Y',&
' x(2,1)*Y x(2,2)*Y x(2,3)*Y >',&
'',&
' The five-point discrete Laplacian for an n-by-n grid can be',&
' generated by',&
'',&
' T = diag(ones(n-1,1),1); T = T + T''; I = eye(T);',&
' A = T.*.I + I.*.T - 4*eye;',&
'',&
' Just in case they might be useful, LALA includes',&
' constructions called Kronecker tensor quotients, denoted by',&
' X ./. Y and X .\. Y . They are obtained by replacing the',&
' element-wise multiplications in X .*. Y with divisions.',&
'',&
'lu Factors from Gaussian elimination. <L,U> = LU(X) stores a',&
' upper triangular matrix in U and a ''psychologically lower',&
' triangular matrix'', i.e. a product of lower triangular and',&
' permutation matrices, in L , so that X = L*U . By itself,',&
' "lu(X)" returns the output from CGEFA .',&
'',&
'magic Magic square. "magic(N)" is an N by N matrix constructed',&
' from the integers 1 through N**2 with equal row, column and',&
' diagonal sums. N must be a positive whole number not equal to two.',&
'',&
'norm computes the norm or P-norm of X',&
'',&
' norm(X,P) computes the P-norm of X. P=2 is the default, which defines',&
' the standard norm.',&
'',&
' For matrices..',&
' norm(X,1) is the 1-norm of X; ie. the largest column sum',&
' of X.',&
'',&
' norm(X,2) the largest singular value of X.',&
' or norm(X)',&
'',&
' norm(X,''inf'') is the infinity norm of X; ie. the largest row',&
' sum of X.',&
'',&
' norm(X,''fro'') is the F-norm, i.e. "sqrt(sum(diag(X''*X)))" .',&
'',&
' For vectors..',&
' norm(V,P) the same as sum(V(I)**P)**(1/P) .',&
' ??? what about negative values of (I) and odd P? abs() or not',&
'',&
' norm(V,2) the square root of the sum of the squares of',&
' or norm(V) the entries of V.',&
'',&
' norm(V,''inf'') the value is max(abs(V)) .',&
'<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<',&
'!! If X is a vector, then',&
'!!',&
'!! norm(x,p) = sum(abs(x) .^ p) ^ (1/p)',&
'!! norm(x,1) is the sum of the absolute values of X.',&
'!! norm(x)/sqrt(n) is the root-mean-square value.',&
'!! norm(x,-inf)=min(abs(x))',&
'<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<',&
'',&
'ones All ones. "ones(N)" is an N by N matrix of ones. "ones(M,N)"',&
' is an M by N matrix of ones . "ones(A)" is the same size as A and',&
' all ones .',&
'',&
' a=magic(4)',&
' a=a+ones(a)*3 // Add 3 to each element of "a"',&
'',&
'orth Orthogonalization. "Q = orth(X)" is a matrix with',&
' orthonormal columns, i.e. Q''*Q = eye, which span the same',&
' space as the columns of X .',&
'',&
'pinv Pseudoinverse.',&
'',&
' "X = pinv(A)" produces a matrix X of the same dimensions as A''',&
' so that A*X*A = A , X*A*X = X and AX and XA are Hermitian . The',&
' computation is based on "svd(A)" and any singular values less',&
' than a tolerance are treated as zero. The default tolerance is',&
' "norm(shape(A),''inf'')*norM(A)*eps". This tolerance may be overridden',&
' with "X = pinv(A,tol)". See "rank".',&
'',&
'poly Characteristic polynomial.',&
'',&
' If A is an N by N matrix, "poly(A)" is a column vector with',&
' N+1 elements which are the coefficients of the characteristic',&
' polynomial, "det(lambda*eye - A)" .',&
'',&
' If V is a vector, "poly(V)" is a vector whose elements are the',&
' coefficients of the polynomial whose roots are the elements of V',&
' . For vectors, "roots" and "poly" are inverse functions of each',&
' other, up to ordering, scaling, and roundoff error.',&
'',&
' "roots(poly(1:20))" generates Wilkinson''s famous example.',&
'',&
'prod "prod(X)" is the product of all the elements of X .',&
'',&
'qr Orthogonal-triangular decomposition. "<Q,R> = qr(X)" produces an',&
' upper triangular matrix R of the same',&
' dimension as X and a unitary matrix Q so that X = Q*R .',&
'',&
' "<Q,R,E> = qr(X)" produces a permutation matrix E, an upper',&
' triangular R with decreasing diagonal elements and a unitary Q',&
' so that X*E = Q*R . By itself, "qr(X)" returns the output of',&
' "cqrdc". "triu(qr(X))" is R .',&
'',&
'rand Random numbers and matrices. "rand(N)" is an N by N matrix',&
' with random entries. "rand(M,N)" is an M by N matrix with',&
' random entries. "rand(A)" is the same size as A. "rand"',&
' with no arguments is a scalar whose value changes each time',&
' it is referenced.',&
'',&
' Ordinarily, random numbers are uniformly distributed in',&
' the interval "(0.0,1.0). rand(''normal'')" switches to a',&
' normal distribution with mean 0.0 and variance 1.0.',&
' "rand(''uniform'')" switches back to the uniform distribution.',&
' "rand(''seed'')" returns the current value of the seed for the',&
' generator. "rand(''seed'',n)" sets the seed to n.',&
' "rand(''seed'',0)" resets the seed to 0, its value when LALA',&
' is first entered.',&
'',&
'rank Rank. "K = rank(X)" is the number of singular values of X',&
' that are larger than "norm(shape(X),''inf'')*norm(X)*eps".',&
' "K = rank(X,tol)" is the number of singular values of X that',&
' are larger than tol.',&
'',&
'rcond "rcond(X)" is an estimate for the reciprocal of the',&
' condition of X in the 1-norm obtained by the LINPACK',&
' condition estimator. If X is well conditioned, rcond(X)',&
' is near 1.0. If X is badly conditioned, rcond(X) is',&
' near 0.0.',&
' <R, Z> = rcond(A) sets R to rcond(A) and also produces a',&
' vector Z so that',&
'',&
' norm(A*Z,1) = R*norm(A,1)*norm(Z,1)',&
'',&
' So, if rcond(A) is small, then Z is an approximate null',&
' vector.',&
'',&
'rat An experimental function which attempts to remove the',&
' roundoff error from results that should be "simple"',&
' rational numbers.',&
' "rat(X)" approximates each element of X by a continued',&
' fraction of the form',&
'',&
' a/b = d1 + 1/(d2 + 1/(d3 + ... + 1/dk))',&
'',&
' with k <= len, integer di and abs(di) <= max . The default',&
' values of the parameters are len = 5 and max = 100.',&
' "rat(len,max)" changes the default values. Increasing either',&
' len or max increases the number of possible fractions.',&
' "<A,B> = rat(X)" produces integer matrices A and B so that',&
'',&
' A ./ B = rat(X)',&
'',&
' Some examples:',&
'',&
' long',&
' T = invh(6), X = inv(T)',&
' <A,B> = rat(X)',&
' H = A ./ B, S = inv(H)',&
'',&
' short e',&
' d = 1:8, e = ones(d), A = abs(d''*e - e''*d)',&
' X = inv(A)',&
' rat(X)',&
' display(ans)',&
'',&
'real "real(X)" is the real part of X.',&
'',&
'rref "rref(A)" is the reduced row echelon form of the rectangular',&
' matrix. rref(A,B) is the same as rref(<A,B>) .',&
'',&
'roots Find polynomial roots. "roots(C)" computes the roots of the',&
' polynomial whose coefficients are the elements of the vector C.',&
' If C has N+1 components, the polynomial is',&
'',&
' C(1)*X**N + ... + C(N)*X + C(N+1)',&
'',&
' See "poly".',&
'',&
'round "round(X)" rounds the elements of X to the nearest integers.',&
'',&
'schur Schur decomposition. "<U,T> = schur(X)" produces an upper',&
' triangular matrix T , with the eigenvalues of X on the',&
' diagonal, and a unitary matrix U so that X = U*T*U'' and',&
' U''*U = eye . By itself, "schur(X)" returns T .',&
'',&
'shape If X is an M by N matrix, then shape(X) is <M, N> .',&
' Can also be used with a multiple assignment,',&
' <M, N> = shape(X) .',&
'',&
'sum "sum(X)" is the sum of all the elements of X.',&
' "sum(diag(X))" is the trace of X.',&
'',&
'svd Singular value decomposition. "<U,S,V> = svd(X)" produces a',&
' diagonal matrix S , of the same dimension as X and with',&
' nonnegative diagonal elements in decreasing order, and',&
' unitary matrices U and V so that X = U*S*V'' .',&
'',&
' By itself, "svd(X)" returns a vector containing the singular',&
' values.',&
'',&
' "<U,S,V> = svd(X,0)" produces the "economy size"',&
' decomposition. If X is m by n with m > n, then only the',&
' first n columns of U are computed and S is n by n .',&
'',&
'tril Lower triangle. "tril(X)" is the lower triangular part of X.',&
' "tril(X,K)" is the elements on and below the K-th diagonal of',&
' X. K = 0 is the main diagonal, K > 0 is above the main',&
' diagonal and K < 0 is below the main diagonal.',&
'',&
'triu Upper triangle. "triu(X)" is the upper triangular part of X.',&
' "triu(X,K)" is the elements on and above the K-th diagonal of X. K',&
' = 0 is the main diagonal, K > 0 is above the main diagonal and K <',&
' 0 is below the main diagonal.',&
'',&
'user Allows personal Fortran subroutines to be linked into',&
' LALA. The subroutine should have the heading',&
'',&
' SUBROUTINE USER(A,M,N,S,T)',&
' REAL or DOUBLEPRECISION A(M,N),S,T',&
'',&
' The LALA statement "Y = user(X,s,t)" results in a call to the',&
' subroutine with a copy of the matrix X stored in the argument A,',&
' its column and row dimensions in M and N, and the scalar parameters',&
' s and t stored in S and T. If s and t are omitted, they are set',&
' to 0.0. After the return, A is stored in Y. The dimensions M and',&
' N may be reset within the subroutine. The statement Y = "user(K)"',&
' results in a call with M = 1, N = 1 and A(1,1) = "float(K)". After',&
' the subroutine has been written, it must be compiled and linked',&
' to the LALA object code within the local operating system.',&
'',&
'zeros',&
' Returns a matrix of all zeros.',&
'',&
' zeros(N) returns an N by N matrix of zeroes.',&
' zeros(M,N) returns an M by N matrix of zeroes.',&
' zeros(X) returns a matrix of zeroes of the same order as X.',&
'================================================================================',&
'FLOW CONTROL',&
'',&
'else Used with "if".',&
'',&
'end Terminates the scope of "for", "while" and "if" statements.',&
' Without "end"s, "for" and "while" repeat all statements up to',&
' the end of the line. Each "end" is paired with the closest',&
' previous unpaired "for" or "while" and serves to terminate its',&
' scope. The line',&
'',&
' for I=1:N, for J=1:N, A(I,J)=1/(I+J-1); A',&
'',&
' would cause A to be printed N**2 times, once for each new',&
' element. On the other hand, the line',&
'',&
' for I=1:N, for J=1:N, A(I,J)=1/(I+J-1); end, end, A',&
'',&
' will lead to only the final printing of A.',&
' Similar considerations apply to "while".',&
'',&
' See "exit" (terminates execution of loops or of LALA itself).',&
'',&
'if Conditionally execute statements',&
'',&
' SIMPLE FORM',&
' Enter',&
'',&
' if expression rop expression, statements',&
'',&
' where rop is =, <, >, <=, >=, or <> (not equal). The',&
' statements are executed once if the indicated comparison',&
' between the real parts of the first components of the two',&
' expressions is true, otherwise the statements are skipped.',&
'',&
' EXAMPLE',&
' Enter',&
'',&
' if abs(i-j) = 1, a(i,j) = -1;',&
'',&
' More complicated forms use "end" in the same way it is used with',&
' "for" and "while" and use "else" as an abbreviation for "end",',&
'',&
' if expression not rop expression',&
'',&
' EXAMPLE',&
' Enter',&
'',&
' for i = 1:n, for j = 1:n, ...',&
' if i = j, a(i,j) = 2; else if abs(i-j) = 1, a(i,j) = -1; ...',&
' else a(i,j) = 0;',&
'',&
' An easier way to accomplish the same thing is',&
'',&
' a = 2*eye(n);',&
' for i = 1:n-1, a(i,i+1) = -1; a(i+1,i) = -1;',&
'',&
'for Repeat statements a specific number of times.',&
'',&
' for variable = expr, statement, ..., statement, end',&
'',&
' The "end" at the end of a line may be omitted. The comma before the',&
' "end" may also be omitted. The columns of the expression are stored',&
' one at a time in the variable and then the following statements,',&
' up to the "end", are executed. The expression is often of the form',&
' X:Y, in which case its columns are simply scalars. Some examples',&
' (assume N has already been assigned a value).',&
'',&
' for I = 1:N, for J = 1:N, A(I,J) = 1/(I+J-1);',&
' for J = 2:N-1, A(J,J) = J; end; A',&
' for S = 1.0: -0.1: 0.0, ... steps S with increments of -0.1 .',&
' for E = eye(N), ... sets E to the unit N-vectors.',&
' for V = A, ... has the same effect as',&
' for J = 1:N, V = A(:,J); ... except J is also set here.',&
'',&
'while Repeat statements an indefinite number of times.',&
'',&
' while expr rop expr, statement, ..., statement, end',&
'',&
' where rop is =, <, >, <=, >=, or <> (not equal). The "end"',&
' at the end of a line may be omitted. The comma before the',&
' "end" may also be omitted. The commas may be replaced by',&
' semicolons to avoid printing. The statements are',&
' repeatedly executed as long as the indicated comparison',&
' between the real parts of the first components of the two',&
' expressions is true.',&
'',&
' EXAMPLE',&
' (assume a matrix A is already defined).',&
'',&
' E = 0*A; F = E + eye; N = 1;',&
' while norm(E+F-E,1) > 0, E = E + F; F = A*F/N; N = N + 1;',&
' E',&
'',&
'exit Causes termination of a "for" or "while" loop.',&
' If not in a loop, terminates execution of LALA.',&
' Also see "quit".',&
'',&
'quit From the terminal, causes return to the operating system',&
' or other program which invoked LALA. From inside an',&
' "exec", causes return to the invoking "exec", or to the',&
' terminal.',&
'================================================================================',&
'FILE ACCESS',&
'',&
' The "exec", "save", "delete", "load", "diary", and "print"',&
' functions access files. The ''file'' parameter takes different',&
' forms for different operating systems. On most systems, ''file''',&
' may be a string of up to 1024 characters in quotes. For example,',&
' "save(''A'')" or "exec(''LALA/demo.exec'')" . The string will be used',&
' as the name of a file in the local operating system.',&
'',&
' Check your local installation for details. The filename must be',&
' composed of recognized characters. See "char".',&
'',&
' Also see "quit" and "exit".',&
'',&
'delete "delete(''filename'')" deletes the given file.',&
'',&
'exec "exec(''file'',k)" obtains subsequent LALA input from an',&
' external file. The printing of input is controlled by the',&
' optional parameter k .',&
'',&
' Files are searched for by the given name. If not found, it is searched',&
' for in the colon-separated directory names in the environment variable',&
' LALA_PATH. It is looked for first literally by the given name, and then',&
' by the name suffixed with ".la".',&
'',&
' "include" is an alias for "exec".',&
'',&
' If k = 0 , there is no echo, prompt or pause. This is the',&
' default if the exec command is followed by a semicolon.',&
' If k = 1 , the input is echoed.',&
' If k = 2 , the LALA prompt <> is printed.',&
' If k = 3 , there will be echos and prompts, but no pauses.',&
' This is the the default if the exec command is not',&
' followed by a semicolon.',&
' If k = 4 , LALA pauses before each prompt and waits for a',&
' null line to continue.',&
' If k = 7 , there will be echos, prompts and pauses. This is',&
' useful for demonstrations on video terminals.',&
'',&
' "exec"''s may be nested, i.e. the text in the file may contain',&
' "exec" of another file.',&
'',&
' "exec" may not be recursive, as Fortran (currently) does not allow',&
' for multiple opens of a file.',&
'',&
' "exec"s may also be driven by "for" and "while" loops.',&
'',&
'include "include" is an alias for "exec".',&
'',&
'load "load(''file'')" retrieves all the variables from the file .',&
' See FILE and "save" for more details. To prepare your own',&
' file for "load"ing, change the "read" to "write" in the code',&
' given under "save".',&
'',&
'print "print(''file'',X)" prints X on the file using the current',&
' format determined by "short", "long z", etc. See FILE.',&
'',&
'doc does nothing at the moment',&
'',&
'save "save(''file'')" stores all the current variables in a file.',&
' "save(''file'',X)" saves only X . See FILE .',&
'',&
' The variables may be retrieved later by "load(''file'')" or by your',&
' own program using the following code for each matrix. The lines',&
' involving "ximag" may be eliminated if everything is known to',&
' be real.',&
'',&
' > ! attach LUN to ''file'', then ...',&
' > doubleprecision :: xreal(mmax,nmax)',&
' > doubleprecision :: ximag(mmax,nmax)',&
' > character(len=32) :: id',&
' > read(LUN,''(a32,3i9)'') id,m,n,img',&
' > do j = 1, n',&
' > read(LUN,''(4z18)'') (xreal(i,j), i=1,m)',&
' > if (img .ne. 0) read(LUN,102) (ximag(i,j),i=1,m)',&
' > enddo',&
' > ! The formats used are system dependent. These are typical.',&
' > ! See SUBROUTINE mat_savlod(3f) in your local implementation',&
' > ! of LALA.',&
'',&
'================================================================================',&
'OUTPUT OPTIONS',&
' ( Also see "FILE" ("exec", "load", "print", "save" ))',&
'',&
'lines An internal count is kept of the number of lines of output',&
' since the last input. Whenever this count approaches a',&
' limit, the user is asked whether or not to suppress',&
' printing until the next input. Initially the limit is 21.',&
' "lines(N)" resets the limit to N .',&
'',&
'long See "short" also.',&
'',&
' Determine output format. All computations are done in',&
' complex arithmetic and double precision if it is available.',&
' "short" and "long" merely switch between different output',&
' formats.',&
'',&
' long // Scaled fixed point format with about 15 digits.',&
' long e // Floating point format with about 15 digits.',&
' long z // System dependent format, often hexadecimal.',&
'',&
'short See "long" also.',&
' Determine output format. All computations are done in',&
' complex arithmetic and double precision if it is available.',&
' "short" and "long" merely switch between different output',&
' formats.',&
'',&
' short // Scaled fixed point format with about 5 digits.',&
' short e // Floating point format with about 5 digits.',&
'',&
'diary "diary(''file'')" causes a copy of all subsequent terminal input and',&
' most of the resulting output to be written on the file. "diary(0)"',&
' turns it off. See "FILE".',&
'',&
'display "display(X)" prints X in a compact format.',&
'',&
' If base >= 2 is specified the values are printed as numeric',&
' values in the specified base.',&
'',&
' display(0:10,10 ) // display values in base 10',&
' display(0:10,16 ) // display values as hexadecimal values',&
' display(0:10,2 ) // display values as binary numbers',&
'',&
' If no base is specified and all the elements of X are integers',&
' between 0 and 255, then X is interpreted as LALA text and printed',&
' accordingly.',&
'',&
' <>display(''the analysis is complete'')',&
' the analysis is complete',&
' display(32:126) // print the printable default LALA characters',&
'',&
' Otherwise or if the base is one, + , - and blank are printed for',&
' positive, negative and zero elements.',&
'',&
' display(rand(24,80)-rand(24,80))',&
'',&
' Imaginary parts are ignored.',&
'',&
' Note that "display(X,B)" is the same as "display(base(X,B))" except',&
' for base 1 except it forces "display" to display numeric values.',&
'',&
' "display" has an alias of "disp".',&
'',&
'plot "plot(X,Y)" produces a plot of the elements of Y against',&
' those of X. plot(Y) is the same as plot(1:n,Y) where n is the number',&
' of elements in Y. plot(X,Y,P) or "plot(X,Y,p1,...,pk)" passes the',&
' optional parameter vector P or scalars p1 through pk to the plot',&
' routine. The default plot routine is a crude printer-plot. This',&
' version writes the data as a simple X Y table into a scratch file',&
' called "scratch.dat" and then calls the command',&
'',&
' xy scratch.dat [options]',&
'',&
' Hopefully, you have the command xy(1) in your search path.',&
' If not, you can make one by creating a script that calls',&
' a plotting utility.',&
'',&
' t = 0:50;',&
' plot( t.*cos(t), t.*sin(t) )',&
' opts='' -m -1 -title test plot -d pdf''',&
' plot( t.*cos(t), t.*sin(t),opts)',&
'================================================================================',&
'ENVIRONMENT',&
'',&
'getenv get environment variable or return a space',&
'',&
' // read commands from a file if an environment variable is set.',&
' MATRC=getenv(''MATRC'');',&
' if MATRC <> '' '', exec(''MATRC'');',&
'================================================================================',&
'PERFORMANCE INFORMATION',&
'',&
'flops Count of floating point operations.',&
'',&
' "flops" is a permanently defined row vector with two elements.',&
' "flops(1)" is the cpu time consumed by the the previous',&
' statement. "flops(2)" is a cumulative total. "flops" can be used',&
' in the same way as any other vector. "flops(2) = 0" resets the',&
' cumulative total. In addition, "flops(1)" will be printed whenever',&
' a statement is terminated by an extra comma. For example,',&
'',&
' X = inv(A);,',&
'',&
' or',&
'',&
' cond(A), (as the last statement on the line).',&
'================================================================================',&
'MISCELLANEOUS',&
'',&
'CHAR special issues regarding strings',&
'',&
' LALA has a limited facility for handling text. Any string of',&
' characters delineated by quotes (with two quotes used to allow one',&
' quote within the string) is saved as a vector of integer values that',&
' are the ADE (Ascii Decimal Equivalent) value of the character.',&
'',&
' In commands { and } are equivalent to ( and )',&
'',&
' When defining an array [ and ] or < and > may be used as the delimiters.',&
'',&
' lala(3f) is too flexible about that and lets them be interchanged freely',&
' instead of being matched but that will probably change to be more strictly',&
' enforced.',&
'',&
' Currently " is not a special character but will probably be allowed as a',&
' string quoting character in the future.',&
'',&
' For example',&
'',&
' ''2*A + 3'' // is the same as < 50 42 65 32 43 32 51 >.',&
' display(32:126) // display the basic visible ASCII characters',&
'',&
'',&
' So if you wanted to home the cursor and clear the screen on an',&
' ANSI-compatible terminal and entered',&
'',&
' display(<27,''[H'',27,''[2J''>)',&
'',&
' The terminal screen would clear. More usefully, if you define the',&
' string',&
'',&
' clr=''display([27,91,''''H'''',27,91,''''2J''''])''',&
'',&
' Then entering',&
'',&
' >clr',&
'',&
' will clear the screen on ANSI terminals and emulators.',&
'',&
'DECIMAL ADE TABLE',&
' The ASCII Decimal Equivalents',&
' *-------*-------*-------*-------*-------*-------*-------*-------*',&
' | 00 nul| 01 soh| 02 stx| 03 etx| 04 eot| 05 enq| 06 ack| 07 bel|',&
' | 08 bs | 09 ht | 10 nl | 11 vt | 12 np | 13 cr | 14 so | 15 si |',&
' | 16 dle| 17 dc1| 18 dc2| 19 dc3| 20 dc4| 21 nak| 22 syn| 23 etb|',&
' | 24 can| 25 em | 26 sub| 27 esc| 28 fs | 29 gs | 30 rs | 31 us |',&
' | 32 sp | 33 ! | 34 " | 35 # | 36 $ | 37 % | 38 & | 39 '' |',&
' | 40 ( | 41 ) | 42 * | 43 + | 44 , | 45 - | 46 . | 47 / |',&
' | 48 0 | 49 1 | 50 2 | 51 3 | 52 4 | 53 5 | 54 6 | 55 7 |',&
' | 56 8 | 57 9 | 58 : | 59 ; | 60 < | 61 = | 62 > | 63 ? |',&
' | 64 @ | 65 A | 66 B | 67 C | 68 D | 69 E | 70 F | 71 G |',&
' | 72 H | 73 I | 74 J | 75 K | 76 L | 77 M | 78 N | 79 O |',&
' | 80 P | 81 Q | 82 R | 83 S | 84 T | 85 U | 86 V | 87 W |',&
' | 88 X | 89 Y | 90 Z | 91 [ | 92 \ | 93 ] | 94 ^ | 95 _ |',&
' | 96 ` | 97 a | 98 b | 99 c |100 d |101 e |102 f |103 g |',&
' |104 h |105 i |106 j |107 k |108 l |109 m |110 n |111 o |',&
' |112 p |113 q |114 r |115 s |116 t |117 u |118 v |119 w |',&
' |120 x |121 y |122 z |123 { |124 | |125 } |126 ~ |127 del|',&
' *-------*-------*-------*-------*-------*-------*-------*-------*',&
'',&
'?? Two exclamation marks beginning a line enters command history mode.',&
' The rest of the line is treated as an optional initial history',&
' edit command. Enter "???" to enter history mode and then display',&
' additional instructions.',&
' see "EDIT" for further details.',&
'',&
'EDIT',&
' A command line consisting of two question marks("??") will cause a',&
' small line-based editor to be called (very similar to the CDC NOS',&
' editor "xedit") with a copy of the previous input lines. When the',&
' editor returns control to LALA, it will execute the edited command',&
' (by default).',&
'',&
' In editor mode the command to be edited is shifted over one and the',&
' first character of input determines the edit mode. The letter "c"',&
' begins a string change (ie. "c/oldstring/newstring/"). The "l"',&
' command lists the lines. A number goes to that command line as',&
' listed by the "l" command. If the change command begins with a',&
' space a letter replaces the one above it with the exception of',&
' the special characters # (delete) & (blank out) and ^ (insert the',&
' following string here).',&
'',&
' An editing loop is entered until a carriage return on an empty',&
' line is entered to accept the new line or a period is entered to',&
' cancel the editing.',&
'',&
' For example, if you had entered a line such as:',&
'',&
' <M,N>=shape(A);for I = 1:M, for J = 1:N, A(I,J) = A(I,J)+3.6;',&
'',&
' Then to repeat the command changing "3.6" to "5.1" enter',&
'',&
' ??',&
' the previous command is then displayed. Now enter',&
'',&
' c/3.6/5.1',&
'',&
' and then enter a carriage return and the edited line will be',&
' executed.',&
'',&
' The first command can appear on the same line if the line starts',&
' with "?? " (two question marks followed by a space). For example',&
'',&
' ?? /rand',&
'',&
' would take you into edit mode on the last command containing the',&
' string "rand"',&
'',&
' Enter "?" in edit mode to display further help on editor mode.',&
'',&
'eps Floating point relative accuracy. A permanent variable',&
' whose value is initially the distance from 1.0 to the next largest',&
' floating point number. The value is changed by "chop", and other',&
' values may be assigned. "eps" is used as a default tolerance by "pinv"',&
' and "rank".',&
'',&
'lala A placeholder for a new command.',&
'',&
'debug "debu(1)" turns on verbose low-level debugging for the developer,',&
' "debu(0)" turns it back off.',&
'',&
'================================================================================',&
'']
end subroutine mat_help_text