Return to doc CVS log

Up to [Research Unix] / researchv10dc / cmd / matlab

Annotation of researchv10dc/cmd/matlab/doc, revision 1.1

1.1     ! root        1: 
        !             2: 
        !             3: 6/24/81
        !             4: 
        !             5: 
        !             6: 
        !             7: 
        !             8: 
        !             9: 
        !            10: 
        !            11: 
        !            12: 
        !            13: 
        !            14: 
        !            15:                        MATLAB Users' Guide
        !            16:                             May, 1981
        !            17: 
        !            18: 
        !            19:                            Cleve Moler
        !            20:                  Department of Computer Science
        !            21:                     University of New Mexico
        !            22: 
        !            23: 
        !            24:      ABSTRACT.  MATLAB is an  interactive  computer  program
        !            25:      that   serves   as   a   convenient   "laboratory"  for
        !            26:      computations  involving  matrices.   It  provides  easy
        !            27:      access  to matrix software developed by the LINPACK and
        !            28:      EISPACK projects.  The program is  written  in  Fortran
        !            29:      and  is  designed  to  be  readily  installed under any
        !            30:      operating system which permits interactive execution of
        !            31:      Fortran programs.
        !            32: 
        !            33: 
        !            34: 
        !            35:                             CONTENTS
        !            36: 
        !            37:           1.  Elementary operations              page  2
        !            38:           2.  MATLAB functions                         8
        !            39:           3.  Rows, columns and submatrices            9
        !            40:           4.  FOR, WHILE and IF                       10
        !            41:           5.  Commands, text, files and macros        12
        !            42:           6.  Census example                          13
        !            43:           7.  Partial differential equation           19
        !            44:           8.  Eigenvalue sensitivity example          23
        !            45:           9.  Syntax diagrams                         27
        !            46:          10.  The parser-interpreter                  31
        !            47:          11.  The numerical algorithms                34
        !            48:          12.  FLOP and CHOP                           37
        !            49:          13.  Communicating with other programs       41
        !            50:          Appendix.  The HELP file                     46
        !            51: 
        !            52: 
        !            53: 
        !            54: 
        !            55: 
        !            56: 
        !            57: 
        !            58: 
        !            59: 
        !            60: 
        !            61: 
        !            62: 
        !            63: 
        !            64: 
        !            65: 
        !            66: 
        !            67: 
        !            68: 
        !            69: 6/24/81
        !            70: 
        !            71: 
        !            72: 
        !            73: 
        !            74: 
        !            75: 
        !            76: 
        !            77: 
        !            78: 
        !            79: 
        !            80: 
        !            81:                        MATLAB Users' Guide
        !            82:                          November, 1980
        !            83: 
        !            84: 
        !            85:                            Cleve Moler
        !            86:                  Department of Computer Science
        !            87:                     University of New Mexico
        !            88: 
        !            89: 
        !            90: 
        !            91:      MATLAB is an interactive computer program that serves  as  a
        !            92: convenient  "laboratory" for computations involving matrices.  It
        !            93: provides easy access to matrix software developed by the  LINPACK
        !            94: and EISPACK projects [1-3].  The capabilities range from standard
        !            95: tasks such as solving simultaneous linear equations and inverting
        !            96: matrices, through symmetric and nonsymmetric eigenvalue problems,
        !            97: to fairly sophisticated matrix tools such as the  singular  value
        !            98: decomposition.
        !            99: 
        !           100:      It is expected that one of MATLAB's primary uses will be  in
        !           101: the  classroom.   It  should be useful in introductory courses in
        !           102: applied linear algebra, as  well  as  more  advanced  courses  in
        !           103: numerical analysis, matrix theory, statistics and applications of
        !           104: matrices to other disciplines.  In nonacademic  settings,  MATLAB
        !           105: can  serve as a "desk calculator" for the quick solution of small
        !           106: problems involving matrices.
        !           107: 
        !           108:      The program is written in Fortran  and  is  designed  to  be
        !           109: readily  installed  under  any  operating  system  which  permits
        !           110: interactive  execution  of  Fortran  programs.    The   resources
        !           111: required  are  fairly  modest.  There are less than 7000 lines of
        !           112: Fortran  source  code,  including   the   LINPACK   and   EISPACK
        !           113: subroutines  used.   With  proper use of overlays, it is possible
        !           114: run the system on a minicomputer with only 32K bytes of memory.
        !           115: 
        !           116:      The size of the matrices  that  can  be  handled  in  MATLAB
        !           117: depends  upon  the  amount  of storage that is set aside when the
        !           118: system is compiled on a particular machine.  We have  found  that
        !           119: an  allocation of 5000 words for matrix elements is usually quite
        !           120: satisfactory.  This provides room for several 20 by 20  matrices,
        !           121: for  example.   One  implementation  on  a  virtual memory system
        !           122: provides 100,000 elements.  Since most  of  the  algorithms  used
        !           123: access  memory  in  a  sequential  fashion,  the  large amount of
        !           124: allocated storage causes no difficulties.
        !           125: 
        !           126: 
        !           127: 
        !           128: 
        !           129: 
        !           130: 
        !           131: 
        !           132: 
        !           133: 
        !           134: 
        !           135: MATLAB, page 2
        !           136: 
        !           137: 
        !           138: 
        !           139:      In some ways, MATLAB  resembles  SPEAKEASY  [4]  and,  to  a
        !           140: lesser  extent, APL.  All are interactive terminal languages that
        !           141: ordinarily accept single-line  commands  or  statements,  process
        !           142: them  immediately,  and  print  the  results.  All have arrays or
        !           143: matrices as principal data types.  But for MATLAB, the matrix  is
        !           144: the  only  data  type  (although  scalars,  vectors  and text are
        !           145: special cases), the underlying system is  portable  and  requires
        !           146: fewer resources, and the supporting subroutines are more powerful
        !           147: and, in some cases, have better numerical properties.
        !           148: 
        !           149:      Together, LINPACK and EISPACK represent the state of the art
        !           150: in software for matrix computation.  EISPACK is a package of over
        !           151: 70 Fortran subroutines for various matrix eigenvalue computations
        !           152: that are based for the most part on Algol procedures published by
        !           153: Wilkinson, Reinsch  and  their  colleagues  [5].   LINPACK  is  a
        !           154: package  of  40  Fortran subroutines (in each of four data types)
        !           155: for solving  and  analyzing  simultaneous  linear  equations  and
        !           156: related matrix problems.  Since MATLAB is not primarily concerned
        !           157: with either execution time  efficiency  or  storage  savings,  it
        !           158: ignores  most  of  the special matrix properties that LINPACK and
        !           159: EISPACK subroutines  use  to  advantage.   Consequently,  only  8
        !           160: subroutines   from  LINPACK  and  5  from  EISPACK  are  actually
        !           161: involved.
        !           162: 
        !           163:      In  more  advanced  applications,  MATLAB  can  be  used  in
        !           164: conjunction  with other programs in several ways.  It is possible
        !           165: to define new MATLAB functions and add them to the system.   With
        !           166: most  operating  systems,  it  is  possible to use the local file
        !           167: system to  pass  matrices  between  MATLAB  and  other  programs.
        !           168: MATLAB  command  and statement input can be obtained from a local
        !           169: file  instead  of  from  the  terminal.   The  most   power   and
        !           170: flexibility  is obtained by using MATLAB as a subroutine which is
        !           171: called by other programs.
        !           172: 
        !           173:      This document first gives an overview  of  MATLAB  from  the
        !           174: user's  point  of  view. Several extended examples involving data
        !           175: fitting, partial differential equations,  eigenvalue  sensitivity
        !           176: and other topics are included.  A formal definition of the MATLAB
        !           177: language and an brief description of the parser  and  interpreter
        !           178: are   given.   The  system  was  designed  and  programmed  using
        !           179: techniques described by Wirth [6], implemented  in  nonrecursive,
        !           180: portable  Fortran.   There  is  a brief discussion of some of the
        !           181: matrix algorithms and of their numerical properties.   The  final
        !           182: section  describes  how  MATLAB  can be used with other programs.
        !           183: The appendix includes the HELP documentation available on-line.
        !           184: 
        !           185: 
        !           186: 1.  Elementary operations
        !           187: 
        !           188: 
        !           189:      MATLAB works with essentially only one  kind  of  object,  a
        !           190: rectangular matrix with complex elements.  If the imaginary parts
        !           191: of the elements are all zero, they  are  not  printed,  but  they
        !           192: 
        !           193: 
        !           194: 
        !           195: 
        !           196: 
        !           197: 
        !           198: 
        !           199: 
        !           200: 
        !           201: MATLAB, page 3
        !           202: 
        !           203: 
        !           204: 
        !           205: still  occupy  storage.   In  some situations, special meaning is
        !           206: attached to 1 by 1 matrices, that is scalars, and to 1 by n and m
        !           207: by 1 matrices, that is row and column vectors.
        !           208: 
        !           209:      Matrices can be introduced into  MATLAB  in  four  different
        !           210: ways:
        !           211:         --  Explicit list of elements,
        !           212:         --  Use of FOR and WHILE statements,
        !           213:         --  Read from an external file,
        !           214:         --  Execute an external Fortran program.
        !           215: 
        !           216:      The explicit list is surrounded by angle brackets,  '<'  and
        !           217: '>', and uses the semicolon ';' to indicate the ends of the rows.
        !           218: For example, the input line
        !           219: 
        !           220:    A = <1 2 3; 4 5 6; 7 8 9>
        !           221: 
        !           222: will result in the output
        !           223: 
        !           224:    A     =
        !           225: 
        !           226:        1.    2.   3.
        !           227:        4.    5.   6.
        !           228:        7.    8.   9.
        !           229: 
        !           230: The matrix A  will  be  saved  for  later  use.   The  individual
        !           231: elements  are separated by commas or blanks and can be any MATLAB
        !           232: expressions, for example
        !           233: 
        !           234:    x = < -1.3, 4/5, 4*atan(1) >
        !           235: 
        !           236: results in
        !           237: 
        !           238:    X     =
        !           239: 
        !           240:      -1.3000   0.8000   3.1416
        !           241: 
        !           242: The elementary functions available include sqrt, log,  exp,  sin,
        !           243: cos, atan, abs, round, real, imag, and conjg.
        !           244: 
        !           245:      Large matrices can be spread  across  several  input  lines,
        !           246: with  the  carriage  returns replacing the semicolons.  The above
        !           247: matrix could also have been produced by
        !           248: 
        !           249:    A = < 1 2 3
        !           250:          4 5 6
        !           251:          7 8 9 >
        !           252: 
        !           253: 
        !           254:      Matrices can be input from the local  file  system.   Say  a
        !           255: file named 'xyz' contains five lines of text,
        !           256: 
        !           257: 
        !           258: 
        !           259: 
        !           260: 
        !           261: 
        !           262: 
        !           263: 
        !           264: 
        !           265: 
        !           266: 
        !           267: MATLAB, page 4
        !           268: 
        !           269: 
        !           270: 
        !           271:    A = <
        !           272:    1 2 3
        !           273:    4 5 6
        !           274:    7 8 9
        !           275:    >;
        !           276: 
        !           277: then the  MATLAB  statement  EXEC('xyz')  reads  the  matrix  and
        !           278: assigns it to A .
        !           279: 
        !           280:      The FOR statement allows the generation  of  matrices  whose
        !           281: elements  are  given  by  simple formulas.  Our example matrix  A
        !           282: could also have been produced by
        !           283: 
        !           284:    for i = 1:3, for j = 1:3, a(i,j) = 3*(i-1)+j;
        !           285: 
        !           286: The semicolon at the end of the  line  suppresses  the  printing,
        !           287: which  in  this  case  would  have  been  nine versions of A with
        !           288: changing elements.
        !           289: 
        !           290:      Several statements may be given  on  a  line,  separated  by
        !           291: semicolons or commas.
        !           292: 
        !           293:      Two  consecutive  periods  anywhere  on  a   line   indicate
        !           294: continuation.   The  periods  and  any  following  characters are
        !           295: deleted, then another line is input  and  concatenated  onto  the
        !           296: previous line.
        !           297: 
        !           298:      Two  consecutive  slashes  anywhere  on  a  line  cause  the
        !           299: remainder  of  the  line  to  be  ignored.   This  is  useful for
        !           300: inserting comments.
        !           301: 
        !           302:      Names of variables are formed by a letter, followed  by  any
        !           303: number of letters and digits, but only the first 4 characters are
        !           304: remembered.
        !           305: 
        !           306:      The special character  prime  (')  is  used  to  denote  the
        !           307: transpose of a matrix, so
        !           308: 
        !           309:    x = x'
        !           310: 
        !           311: changes the row vector above into the column vector
        !           312: 
        !           313:    X     =
        !           314: 
        !           315:      -1.3000
        !           316:       0.8000
        !           317:       3.1416
        !           318: 
        !           319: 
        !           320:      Individual matrix elements may be  referenced  by  enclosing
        !           321: their  subscripts  in  parentheses.  When any element is changed,
        !           322: the entire matrix is reprinted.  For  example,  using  the  above
        !           323: matrix,
        !           324: 
        !           325: 
        !           326: 
        !           327: 
        !           328: 
        !           329: 
        !           330: 
        !           331: 
        !           332: 
        !           333: MATLAB, page 5
        !           334: 
        !           335: 
        !           336: 
        !           337:    a(3,3) = a(1,3) + a(3,1)
        !           338: 
        !           339: results in
        !           340: 
        !           341:    A     =
        !           342: 
        !           343:        1.    2.    3.
        !           344:        4.    5.    6.
        !           345:        7.    8.   10.
        !           346: 
        !           347: 
        !           348:      Addition, subtraction and  multiplication  of  matrices  are
        !           349: denoted  by  +, -, and * .  The operations are performed whenever
        !           350: the matrices have the proper dimensions.  For example,  with  the
        !           351: above  A  and  x,  the  expressions  A  + x and x*A are incorrect
        !           352: because A is 3 by 3 and x is now 3 by 1.  However,
        !           353: 
        !           354:    b = A*x
        !           355: 
        !           356: is correct and results in the output
        !           357: 
        !           358:    B     =
        !           359: 
        !           360:       9.7248
        !           361:      17.6496
        !           362:      28.7159
        !           363: 
        !           364: Note that both upper and lower case letters are allowed for input
        !           365: (on  those  systems  which  have  both),  but  that lower case is
        !           366: converted to upper case.
        !           367: 
        !           368:      There are two "matrix division" symbols in MATLAB, \ and / .
        !           369: (If  your  terminal  does not have a backslash, use $ instead, or
        !           370: see CHAR.) If A and B are matrices, then A\B and  B/A  correspond
        !           371: formally  to left and right multiplication of B by the inverse of
        !           372: A , that is inv(A)*B and B*inv(A), but  the  result  is  obtained
        !           373: directly  without  the computation of the inverse.  In the scalar
        !           374: case, 3\1 and 1/3 have the  same  value,  namely  one-third.   In
        !           375: general,  A\B  denotes the solution X to the equation A*X = B and
        !           376: B/A denotes the solution to X*A = B.
        !           377: 
        !           378:      Left division, A\B, is defined whenever B has as  many  rows
        !           379: as   A  .   If  A  is  square,  it  is  factored  using  Gaussian
        !           380: elimination.   The  factors  are  used  to  solve  the  equations
        !           381: A*X(:,j) = B(:,j) where B(:,j) denotes the j-th column of B.  The
        !           382: result is a matrix X with the same dimensions  as  B.   If  A  is
        !           383: nearly  singular  (according  to the LINPACK condition estimator,
        !           384: RCOND), a warning message is printed.  If A is not square, it  is
        !           385: factored   using   Householder   orthogonalization   with  column
        !           386: pivoting.   The  factors  are  used  to  solve  the   under-   or
        !           387: overdetermined equations in a least squares sense.  The result is
        !           388: an m by n matrix X where m is the number of columns of A and n is
        !           389: the  number  of  columns  of  B .  Each column of X has at most k
        !           390: 
        !           391: 
        !           392: 
        !           393: 
        !           394: 
        !           395: 
        !           396: 
        !           397: 
        !           398: 
        !           399: MATLAB, page 6
        !           400: 
        !           401: 
        !           402: 
        !           403: nonzero components, where k is the effective rank of A .
        !           404: 
        !           405:      Right division,  B/A,  can  be  defined  in  terms  of  left
        !           406: division by  B/A = (A'\B')'.
        !           407: 
        !           408:      For example, since our vector  b   was computed as  A*x, the
        !           409: statement
        !           410: 
        !           411:    y = A\b
        !           412: 
        !           413: results in
        !           414: 
        !           415:    Y     =
        !           416: 
        !           417:      -1.3000
        !           418:       0.8000
        !           419:       3.1416
        !           420: 
        !           421: Of course,  y  is  not  exactly  equal  to   x   because  of  the
        !           422: roundoff  errors involved in both  A*x  and  A\b , but we are not
        !           423: printing enough digits to see the difference.  The result of  the
        !           424: statement
        !           425: 
        !           426:    e = x - y
        !           427: 
        !           428: depends upon the particular computer being used.  In one case  it
        !           429: produces
        !           430: 
        !           431:    E     =
        !           432: 
        !           433:       1.0e-15 *
        !           434: 
        !           435:         .3053
        !           436:        -.2498
        !           437:         .0000
        !           438: 
        !           439: The quantity 1.0e-15 is a scale factor which multiplies  all  the
        !           440: components  which  follow.  Thus our vectors  x  and  y  actually
        !           441: agree to about 15 decimal places on this computer.
        !           442: 
        !           443:      It   is   also   possible   to   obtain   element-by-element
        !           444: multiplicative  operations.  If A and B have the same dimensions,
        !           445: then A .* B denotes the matrix  whose  elements  are  simply  the
        !           446: products  of the individual elements of A and B . The expressions
        !           447: A ./ B and A .\ B give the quotients of the individual elements.
        !           448: 
        !           449:      There are several possible output formats.  The statement
        !           450: 
        !           451:    long, x
        !           452: 
        !           453: results in
        !           454: 
        !           455:    X     =
        !           456: 
        !           457: 
        !           458: 
        !           459: 
        !           460: 
        !           461: 
        !           462: 
        !           463: 
        !           464: 
        !           465: MATLAB, page 7
        !           466: 
        !           467: 
        !           468: 
        !           469:       -1.300000000000000
        !           470:         .800000000000000
        !           471:        3.141592653589793
        !           472: 
        !           473: The statement
        !           474: 
        !           475:    short
        !           476: 
        !           477: restores the original format.
        !           478: 
        !           479:      The expression A**p means  A  to  the  p-th  power.   It  is
        !           480: defined  if  A  is a square matrix and p is a scalar.  If p is an
        !           481: integer greater than one,  the  power  is  computed  by  repeated
        !           482: multiplication.   For  other values of p the calculation involves
        !           483: the eigenvalues and eigenvectors of A.
        !           484: 
        !           485:      Previously defined matrices and matrix  expressions  can  be
        !           486: used inside brackets to generate larger matrices, for example
        !           487: 
        !           488:    C = <A, b; <4 2 0>*x, x'>
        !           489: 
        !           490: results in
        !           491: 
        !           492: 
        !           493:    C     =
        !           494: 
        !           495:       1.0000   2.0000   3.0000   9.7248
        !           496:       4.0000   5.0000   6.0000  17.6496
        !           497:       7.0000   8.0000  10.0000  28.7159
        !           498:      -3.6000  -1.3000   0.8000   3.1416
        !           499: 
        !           500: 
        !           501:      There are four predefined variables,  EPS,  FLOP,  RAND  and
        !           502: EYE.  The variable EPS is used as a tolerance is determining such
        !           503: things as near singularity and rank.  Its initial  value  is  the
        !           504: distance  from  1.0  to the next largest floating point number on
        !           505: the particular computer being used.  The user may reset  this  to
        !           506: any other value, including zero. EPS is changed by CHOP, which is
        !           507: described in section 12.
        !           508: 
        !           509:      The value of RAND is a random variable, with a choice  of  a
        !           510: uniform or a normal distribution.
        !           511: 
        !           512:      The name EYE is used  in  place  of  I  to  denote  identity
        !           513: matrices  because  I is often used as a subscript or as sqrt(-1).
        !           514: The dimensions of EYE are determined by context.  For example,
        !           515: 
        !           516:    B = A + 3*EYE
        !           517: 
        !           518: adds 3 to the diagonal elements of A and
        !           519: 
        !           520:    X = EYE/A
        !           521: 
        !           522: 
        !           523: 
        !           524: 
        !           525: 
        !           526: 
        !           527: 
        !           528: 
        !           529: 
        !           530: 
        !           531: MATLAB, page 8
        !           532: 
        !           533: 
        !           534: 
        !           535: is one of several ways in MATLAB to invert a matrix.
        !           536: 
        !           537:      FLOP provides a  count  of  the  number  of  floating  point
        !           538: operations, or "flops", required for each calculation.
        !           539: 
        !           540:      A statement may consist of an  expression  alone,  in  which
        !           541: case a variable named ANS is created and the result stored in ANS
        !           542: for possible future use.  Thus
        !           543: 
        !           544:    A\A - EYE
        !           545: 
        !           546: is the same as
        !           547: 
        !           548:    ANS = A\A - EYE
        !           549: 
        !           550: (Roundoff error usually causes this result  to  be  a  matrix  of
        !           551: "small" numbers, rather than all zeros.)
        !           552: 
        !           553:      All computations are done  using  either  single  or  double
        !           554: precision  real  arithmetic,  whichever  is  appropriate  for the
        !           555: particular computer.  There  is  no  mixed-precision  arithmetic.
        !           556: The  Fortran  COMPLEX  data type is not used because many systems
        !           557: create  unnecessary  underflows  and   overflows   with   complex
        !           558: operations and because some systems do not allow double precision
        !           559: complex arithmetic.
        !           560: 
        !           561: 
        !           562: 2.  MATLAB functions
        !           563: 
        !           564:      Much of MATLAB's computational power comes from the  various
        !           565: matrix functions available.  The current list includes:
        !           566: 
        !           567:    INV(A)          - Inverse.
        !           568:    DET(A)          - Determinant.
        !           569:    COND(A)         - Condition number.
        !           570:    RCOND(A)        - A measure of nearness to singularity.
        !           571:    EIG(A)          - Eigenvalues and eigenvectors.
        !           572:    SCHUR(A)        - Schur triangular form.
        !           573:    HESS(A)         - Hessenberg or tridiagonal form.
        !           574:    POLY(A)         - Characteristic polynomial.
        !           575:    SVD(A)          - Singular value decomposition.
        !           576:    PINV(A,eps)     - Pseudoinverse with optional tolerance.
        !           577:    RANK(A,eps)     - Matrix rank with optional tolerance.
        !           578:    LU(A)           - Factors from Gaussian elimination.
        !           579:    CHOL(A)         - Factor from Cholesky factorization.
        !           580:    QR(A)           - Factors from Householder orthogonalization.
        !           581:    RREF(A)         - Reduced row echelon form.
        !           582:    ORTH(A)         - Orthogonal vectors spanning range of A.
        !           583:    EXP(A)          - e to the A.
        !           584:    LOG(A)          - Natural logarithm.
        !           585:    SQRT(A)         - Square root.
        !           586:    SIN(A)          - Trigonometric sine.
        !           587:    COS(A)          - Cosine.
        !           588: 
        !           589: 
        !           590: 
        !           591: 
        !           592: 
        !           593: 
        !           594: 
        !           595: 
        !           596: 
        !           597: MATLAB, page 9
        !           598: 
        !           599: 
        !           600: 
        !           601:    ATAN(A)         - Arctangent.
        !           602:    ROUND(A)        - Round the elements to nearest integers.
        !           603:    ABS(A)          - Absolute value of the elements.
        !           604:    REAL(A)         - Real parts of the elements.
        !           605:    IMAG(A)         - Imaginary parts of the elements.
        !           606:    CONJG(A)        - Complex conjugate.
        !           607:    SUM(A)          - Sum of the elements.
        !           608:    PROD(A)         - Product of the elements.
        !           609:    DIAG(A)         - Extract or create diagonal matrices.
        !           610:    TRIL(A)         - Lower triangular part of A.
        !           611:    TRIU(A)         - Upper triangular part of A.
        !           612:    NORM(A,p)       - Norm with p = 1, 2 or 'Infinity'.
        !           613:    EYE(m,n)        - Portion of identity matrix.
        !           614:    RAND(m,n)       - Matrix with random elements.
        !           615:    ONES(m,n)       - Matrix of all ones.
        !           616:    MAGIC(n)        - Interesting test matrices.
        !           617:    HILBERT(n)      - Inverse Hilbert matrices.
        !           618:    ROOTS(C)        - Roots of polynomial with coefficients C.
        !           619:    DISPLAY(A,p)    - Print base p representation of A.
        !           620:    KRON(A,B)       - Kronecker tensor product of A and B.
        !           621:    PLOT(X,Y)       - Plot Y as a function of X .
        !           622:    RAT(A)          - Find "simple" rational approximation to A.
        !           623:    USER(A)         - Function defined by external program.
        !           624: 
        !           625:      Some of these functions have different interpretations  when
        !           626: the  argument  is  a  matrix  or  a  vector and some of them have
        !           627: additional optional arguments.  Details are  given  in  the  HELP
        !           628: document in the appendix.
        !           629: 
        !           630:      Several of these functions can  be  used  in  a  generalized
        !           631: assignment statement with two or three variables on the left hand
        !           632: side.  For example
        !           633: 
        !           634:    <X,D> = EIG(A)
        !           635: 
        !           636: stores the eigenvectors of A in  the  matrix  X  and  a  diagonal
        !           637: matrix containing the eigenvalues in the matrix D.  The statement
        !           638: 
        !           639:    EIG(A)
        !           640: 
        !           641: simply computes the eigenvalues and stores them in ANS.
        !           642: 
        !           643:      Future versions of MATLAB will probably  include  additional
        !           644: functions, since they can easily be added to the system.
        !           645: 
        !           646: 
        !           647: 
        !           648: 3.  Rows, columns and submatrices
        !           649: 
        !           650: 
        !           651:      Individual elements of a matrix can be  accessed  by  giving
        !           652: their subscripts in parentheses, eg. A(1,2), x(i), TAB(ind(k)+1).
        !           653: An expression used as a  subscript  is  rounded  to  the  nearest
        !           654: 
        !           655: 
        !           656: 
        !           657: 
        !           658: 
        !           659: 
        !           660: 
        !           661: 
        !           662: 
        !           663: MATLAB, page 10
        !           664: 
        !           665: 
        !           666: 
        !           667: integer.
        !           668: 
        !           669:      Individual rows and columns can be accessed  using  a  colon
        !           670: ':' (or a '|') for the free subscript. For example, A(1,:) is the
        !           671: first row of A and A(:,j) is the j-th column.  Thus
        !           672: 
        !           673:    A(i,:) = A(i,:) + c*A(k,:)
        !           674: 
        !           675: adds c times the k-th row of A to the i-th row.
        !           676: 
        !           677:      The colon is used in several other ways in MATLAB,  but  all
        !           678: of the uses are based on the following definition.
        !           679: 
        !           680:    j:k    is the same as  <j, j+1, ..., k>
        !           681:    j:k    is empty if  j > k .
        !           682:    j:i:k  is the same as  <j, j+i, j+2i, ..., k>
        !           683:    j:i:k  is empty if  i > 0 and j > k or if i < 0 and j < k .
        !           684: 
        !           685: The colon is usually used with integers, but it  is  possible  to
        !           686: use arbitrary real scalars as well.  Thus
        !           687: 
        !           688:    1:4  is the same as  <1, 2, 3, 4>
        !           689:    0: 0.1: 0.5 is the same as <0.0, 0.1, 0.2, 0.3, 0.4, 0.5>
        !           690: 
        !           691: 
        !           692:      In general, a subscript can be a vector.  If  X  and  V  are
        !           693: vectors, then X(V) is <X(V(1)), X(V(2)), ..., X(V(n))> . This can
        !           694: also be used with matrices.  If V has m components and  W  has  n
        !           695: components,  then  A(V,W)  is  the  m by n matrix formed from the
        !           696: elements of A whose subscripts are  the  elements  of  V  and  W.
        !           697: Combinations  of the colon notation and the indirect subscripting
        !           698: allow manipulation of various submatrices. For example,
        !           699: 
        !           700:    A(<1,5>,:) = A(<5,1>,:)  interchanges rows 1 and 5 of A.
        !           701:    A(2:k,1:n)  is the submatrix formed from rows 2 through k
        !           702:       and columns 1 through n of A .
        !           703:    A(:,<3 1 2>)  is a permutation of the first three columns.
        !           704: 
        !           705: 
        !           706:      The notation A(:) has a special meaning.  On the right  hand
        !           707: side  of  an assignment statement, it denotes all the elements of
        !           708: A, regarded as a single column.  When an expression  is  assigned
        !           709: to  A(:),  the  current  dimensions  of  A,  rather  than  of the
        !           710: expression, are used.
        !           711: 
        !           712: 
        !           713: 4.  FOR, WHILE and IF
        !           714: 
        !           715: 
        !           716:      The FOR clause allows statements to be repeated  a  specific
        !           717: number of times.  The general form is
        !           718: 
        !           719:    FOR variable = expr,  statement, ..., statement, END
        !           720: 
        !           721: 
        !           722: 
        !           723: 
        !           724: 
        !           725: 
        !           726: 
        !           727: 
        !           728: 
        !           729: MATLAB, page 11
        !           730: 
        !           731: 
        !           732: 
        !           733: The END and the comma before it may be omitted.  In general,  the
        !           734: expression  may be a matrix, in which case the columns are stored
        !           735: one at a time in the variable and the following statements, up to
        !           736: the  END or the end of the line, are executed.  The expression is
        !           737: often of the form j:k, and its "columns" are simply  the  scalars
        !           738: from j to k.  Some examples (assume n has already been assigned a
        !           739: value):
        !           740: 
        !           741:    for i = 1:n, for j = 1:n, A(i,j) = 1/(i+j-1);
        !           742: 
        !           743: generates the Hilbert matrix.
        !           744: 
        !           745:    for j = 2:n-1, for i = j:n-1, ...
        !           746:       A(i,j) = 0; end; A(j,j) = j; end; A
        !           747: 
        !           748: changes all but the "outer edge" of the lower triangle  and  then
        !           749: prints the final matrix.
        !           750: 
        !           751:    for h = 1.0: -0.1: -1.0, (<h, cos(pi*h)>)
        !           752: 
        !           753: prints a table of cosines.
        !           754: 
        !           755:    <X,D> = EIG(A); for v = X, v, A*v
        !           756: 
        !           757: displays eigenvectors, one at a time.
        !           758: 
        !           759:      The  WHILE  clause  allows  statements  to  be  repeated  an
        !           760: indefinite number of times.  The general form is
        !           761: 
        !           762:    WHILE expr relop expr,   statement,..., statement, END
        !           763: 
        !           764: where relop is =, <,  >,  <=,  >=,  or  <>  (not  equal)  .   The
        !           765: statements  are  repeatedly  executed  as  long  as the indicated
        !           766: comparison between the real parts of the first components of  the
        !           767: two  expressions  is true.  Here are two examples.  (Exercise for
        !           768: the reader: What do these segments do?)
        !           769: 
        !           770:    eps = 1;
        !           771:    while 1 + eps > 1, eps = eps/2;
        !           772:    eps = 2*eps
        !           773: 
        !           774:    E = 0*A;  F = E + EYE; n = 1;
        !           775:    while NORM(E+F-E,1) > 0, E = E + F; F = A*F/n; n = n + 1;
        !           776:    E
        !           777: 
        !           778: 
        !           779:      The IF clause allows conditional  execution  of  statements.
        !           780: The general form is
        !           781: 
        !           782:    IF expr relop expr,   statement, ..., statement,
        !           783:       ELSE statement, ..., statement
        !           784: 
        !           785: The first group of statements are executed  if  the  relation  is
        !           786: 
        !           787: 
        !           788: 
        !           789: 
        !           790: 
        !           791: 
        !           792: 
        !           793: 
        !           794: 
        !           795: MATLAB, page 12
        !           796: 
        !           797: 
        !           798: 
        !           799: true  and the second group are executed if the relation is false.
        !           800: The ELSE and the statements following it  may  be  omitted.   For
        !           801: example,
        !           802: 
        !           803:    if abs(i-j) = 2, A(i,j) = 0;
        !           804: 
        !           805: 
        !           806: 5.  Commands, text, files and macros.
        !           807: 
        !           808: 
        !           809:      MATLAB has several commands which control the output  format
        !           810: and the overall execution of the system.
        !           811: 
        !           812:      The HELP command allows on-line access to short portions  of
        !           813: text   describing   various  operations,  functions  and  special
        !           814: characters.   The  entire  HELP  document  is  reproduced  in  an
        !           815: appendix.
        !           816: 
        !           817:      Results are usually printed in a scaled fixed  point  format
        !           818: that shows 4 or 5 significant figures.  The commands SHORT, LONG,
        !           819: SHORT E, LONG E and LONG Z alter the output format,  but  do  not
        !           820: alter the precision of the computations or the internal storage.
        !           821: 
        !           822:      The WHO, WHAT and WHY commands provide information about the
        !           823: functions and variables that are currently defined.
        !           824: 
        !           825:      The CLEAR command erases all variables,  except  EPS,  FLOP,
        !           826: RAND  and  EYE.  The  statement  A = <> indicates that a "0 by 0"
        !           827: matrix is to be stored in A.  This causes A to be erased so  that
        !           828: its storage can be used for other variables.
        !           829: 
        !           830:      The RETURN and EXIT commands cause return to the  underlying
        !           831: operating system through the Fortran RETURN statement.
        !           832: 
        !           833:      MATLAB has a limited facility for handling text.  Any string
        !           834: of characters delineated by quotes (with two quotes used to allow
        !           835: one quote within the string) is saved  as  a  vector  of  integer
        !           836: values  with '1' = 1, 'A' = 10, ' ' = 36, etc. (The complete list
        !           837: is in the appendix under CHAR.) For example
        !           838: 
        !           839:    '2*A + 3'  is the same as  <2 43 10 36 41 36 3>
        !           840: 
        !           841: It is possible,  though  seldom  very  meaningful,  to  use  such
        !           842: strings  in matrix operations.  More frequently, the text is used
        !           843: as a special argument to various functions.
        !           844: 
        !           845:    NORM(A,'inf')    computes the infinity norm of A .
        !           846:    DISPLAY(T)       prints the text stored in T .
        !           847:    EXEC('file')     obtains MATLAB input from an external file.
        !           848:    SAVE('file')     stores all the current variables in a file.
        !           849:    LOAD('file')     retrieves all the variables from a file.
        !           850:    PRINT('file',X)  prints X on a file.
        !           851:    DIARY('file')    makes a copy of the complete MATLAB session.
        !           852: 
        !           853: 
        !           854: 
        !           855: 
        !           856: 
        !           857: 
        !           858: 
        !           859: 
        !           860: 
        !           861: MATLAB, page 13
        !           862: 
        !           863: 
        !           864: 
        !           865: 
        !           866:      The text can also be used in a limited  string  substitution
        !           867: macro  facility.   If a variable, say T, contains the source text
        !           868: for a MATLAB statement or expression, then the construction
        !           869: 
        !           870:    > T <
        !           871: 
        !           872: causes T to be executed or evaluated.  For example
        !           873: 
        !           874:    T = '2*A + 3';
        !           875:    S = 'B = >T< + 5'
        !           876:    A = 4;
        !           877:    > S <
        !           878: 
        !           879: produces
        !           880: 
        !           881:    B     =
        !           882: 
        !           883:       16.
        !           884: 
        !           885: Some other examples are given under MACRO in the appendix.   This
        !           886: facility  is  useful for fairly short statements and expressions.
        !           887: More complicated MATLAB "programs" should use the EXEC facility.
        !           888: 
        !           889:      The operations which access external files cannot be handled
        !           890: in  a  completely  machine-independent manner by portable Fortran
        !           891: code.  It  is  necessary  for  each  particular  installation  to
        !           892: provide  a  subroutine  which associates external text files with
        !           893: Fortran logical unit numbers.
        !           894: 
        !           895: 
        !           896: 6.  Census example
        !           897: 
        !           898: 
        !           899:      Our  first  extended   example   involves   predicting   the
        !           900: population  of  the  United States in 1980 using extrapolation of
        !           901: various fits to the census data from 1900  through  1970.   There
        !           902: are eight observations, so we begin with the MATLAB statement
        !           903: 
        !           904:    n = 8
        !           905: 
        !           906: The values of the dependent variable, the population in millions,
        !           907: can be entered with
        !           908: 
        !           909:    y = < 75.995   91.972  105.711  123.203   ...
        !           910:         131.669  150.697  179.323  203.212>'
        !           911: 
        !           912: In order to produce a reasonably scaled matrix,  the  independent
        !           913: variable,  time,  is transformed from the interval [1900,1970] to
        !           914: [-1.00,0.75].  This can be accomplished directly with
        !           915: 
        !           916:    t = -1.0:0.25:0.75
        !           917: 
        !           918: 
        !           919: 
        !           920: 
        !           921: 
        !           922: 
        !           923: 
        !           924: 
        !           925: 
        !           926: 
        !           927: MATLAB, page 14
        !           928: 
        !           929: 
        !           930: 
        !           931: or in a fancier, but perhaps clearer, way with
        !           932: 
        !           933:    t = 1900:10:1970;   t = (t - 1940*ones(t))/40
        !           934: 
        !           935: Either of these is equivalent to
        !           936: 
        !           937:    t = <-1 -.75 -.50 -.25 0 .25 .50 .75>
        !           938: 
        !           939:      The interpolating polynomial of  degree   n-1   involves  an
        !           940: Vandermonde  matrix  of  order   n   with  elements that might be
        !           941: generated by
        !           942: 
        !           943:    for i = 1:n, for j = 1:n, a(i,j) = t(i)**(j-1);
        !           944: 
        !           945: However, this results in an error caused by 0**0  when  i = 5 and
        !           946: j = 1 .  The preferable approach is
        !           947: 
        !           948:    A = ones(n,n);
        !           949:    for i = 1:n, for j = 2:n, a(i,j) = t(i)*a(i,j-1);
        !           950: 
        !           951: Now the statement
        !           952: 
        !           953:    cond(A)
        !           954: 
        !           955: produces the output
        !           956: 
        !           957:    ANS   =
        !           958: 
        !           959:       1.1819E+03
        !           960: 
        !           961: which indicates that transformation  of  the  time  variable  has
        !           962: resulted in a reasonably well conditioned matrix.
        !           963: 
        !           964:      The statement
        !           965: 
        !           966:    c = A\y
        !           967: 
        !           968: results in
        !           969: 
        !           970:    C     =
        !           971: 
        !           972:      131.6690
        !           973:       41.0406
        !           974:      103.5396
        !           975:      262.4535
        !           976:     -326.0658
        !           977:     -662.0814
        !           978:      341.9022
        !           979:      533.6373
        !           980: 
        !           981: These are the coefficients in the interpolating polynomial
        !           982: 
        !           983:                           n-1
        !           984: 
        !           985: 
        !           986: 
        !           987: 
        !           988: 
        !           989: 
        !           990: 
        !           991: 
        !           992: 
        !           993: MATLAB, page 15
        !           994: 
        !           995: 
        !           996: 
        !           997:       c  + c t + ... + c t
        !           998:        1    2           n
        !           999: 
        !          1000: Our transformation of the time variable has resulted in   t  =  1
        !          1001: corresponding  to  the year 1980.  Consequently, the extrapolated
        !          1002: population is simply the sum of the coefficients.   This  can  be
        !          1003: computed by
        !          1004: 
        !          1005:    p = sum(c)
        !          1006: 
        !          1007: The result is
        !          1008: 
        !          1009:    P     =
        !          1010: 
        !          1011:      426.0950
        !          1012: 
        !          1013: which indicates a 1980 population of over 426 million.   Clearly,
        !          1014: using  the seventh degree interpolating polynomial to extrapolate
        !          1015: even a fairly short distance beyond the end of the data  interval
        !          1016: is not a good idea.
        !          1017: 
        !          1018:      The coefficients in least squares  fits  by  polynomials  of
        !          1019: lower  degree can be computed using fewer than  n  columns of the
        !          1020: matrix.
        !          1021: 
        !          1022:    for k = 1:n, c = A(:,1:k)\y,  p = sum(c)
        !          1023: 
        !          1024: would produce the coefficients of these  fits,  as  well  as  the
        !          1025: resulting  extrapolated  population.   If we do not want to print
        !          1026: all the coefficients, we can simply generate  a  small  table  of
        !          1027: populations  predicted  by  polynomials  of  degrees zero through
        !          1028: seven.  We also compute the maximum deviation between the  fitted
        !          1029: and observed values.
        !          1030: 
        !          1031:    for k = 1:n, X = A(:,1:k);  c = X\y;  ...
        !          1032:       d(k) = k-1;  p(k) = sum(c);  e(k) = norm(X*c-y,'inf');
        !          1033:    <d, p, e>
        !          1034: 
        !          1035: The resulting output is
        !          1036: 
        !          1037:       0   132.7227  70.4892
        !          1038:       1   211.5101   9.8079
        !          1039:       2   227.7744   5.0354
        !          1040:       3   241.9574   3.8941
        !          1041:       4   234.2814   4.0643
        !          1042:       5   189.7310   2.5066
        !          1043:       6   118.3025   1.6741
        !          1044:       7   426.0950   0.0000
        !          1045: 
        !          1046: The zeroth degree fit, 132.7 million, is the result of fitting  a
        !          1047: constant  to  the  data  and  is simply the average.  The results
        !          1048: obtained with polynomials of degree one through four  all  appear
        !          1049: reasonable.   The  maximum  deviation  of  the degree four fit is
        !          1050: 
        !          1051: 
        !          1052: 
        !          1053: 
        !          1054: 
        !          1055: 
        !          1056: 
        !          1057: 
        !          1058: 
        !          1059: MATLAB, page 16
        !          1060: 
        !          1061: 
        !          1062: 
        !          1063: slightly greater than the degree three, even though  the  sum  of
        !          1064: the  squares  of the deviations is less.  The coefficients of the
        !          1065: highest powers in the fits of degree five and six turn out to  be
        !          1066: negative  and  the predicted populations of less than 200 million
        !          1067: are probably unrealistic.  The hopefully absurd prediction of the
        !          1068: interpolating polynomial concludes the table.
        !          1069: 
        !          1070:      We  wish  to  emphasize  that  roundoff   errors   are   not
        !          1071: significant  here.  Nearly identical results would be obtained on
        !          1072: other computers, or with other algorithms.   The  results  simply
        !          1073: indicate   the  difficulties  associated  with  extrapolation  of
        !          1074: polynomial fits of even modest degree.
        !          1075: 
        !          1076:      A stabilized fit by  a  seventh  degree  polynomial  can  be
        !          1077: obtained  using  the  pseudoinverse,  but  it  requires  a fairly
        !          1078: delicate choice of a tolerance. The statement
        !          1079: 
        !          1080:    s = svd(A)
        !          1081: 
        !          1082: produces the singular values
        !          1083: 
        !          1084:    S     =
        !          1085: 
        !          1086:       3.4594
        !          1087:       2.2121
        !          1088:       1.0915
        !          1089:       0.4879
        !          1090:       0.1759
        !          1091:       0.0617
        !          1092:       0.0134
        !          1093:       0.0029
        !          1094: 
        !          1095: We see that the last three singular values are less  than  0.1  ,
        !          1096: consequently,   A   can be approximately by a matrix of rank five
        !          1097: with an error less than 0.1 .  The Moore-Penrose pseudoinverse of
        !          1098: this  rank  five  matrix  is  obtained  from  the  singular value
        !          1099: decomposition with the following statements
        !          1100: 
        !          1101:    c = pinv(A,0.1)*y, p = sum(c), e = norm(a*c-y,'inf')
        !          1102: 
        !          1103: The output is
        !          1104: 
        !          1105: 
        !          1106: 
        !          1107: 
        !          1108: 
        !          1109: 
        !          1110: 
        !          1111: 
        !          1112: 
        !          1113: 
        !          1114: 
        !          1115: 
        !          1116: 
        !          1117: 
        !          1118: 
        !          1119: 
        !          1120: 
        !          1121: 
        !          1122: 
        !          1123: 
        !          1124: 
        !          1125: MATLAB, page 17
        !          1126: 
        !          1127: 
        !          1128: 
        !          1129:    C     =
        !          1130: 
        !          1131:     134.7972
        !          1132:      67.5055
        !          1133:      23.5523
        !          1134:       9.2834
        !          1135:       3.0174
        !          1136:       2.6503
        !          1137:      -2.8808
        !          1138:       3.2467
        !          1139: 
        !          1140:    P     =
        !          1141: 
        !          1142:     241.1720
        !          1143: 
        !          1144:    E     =
        !          1145: 
        !          1146:       3.9469
        !          1147: 
        !          1148: The resulting seventh degree polynomial  has  coefficients  which
        !          1149: are much smaller than those of the interpolating polynomial given
        !          1150: earlier.  The predicted population and the maximum deviation  are
        !          1151: reasonable.   Any  choice  of the tolerance between the fifth and
        !          1152: sixth singular values would produce the same results, but choices
        !          1153: outside this range result in pseudoinverses of different rank and
        !          1154: do not work as well.
        !          1155: 
        !          1156:      The one term exponential approximation
        !          1157: 
        !          1158:      y(t) = k exp(pt)
        !          1159: 
        !          1160: can  be  transformed  into  a  linear  approximation  by   taking
        !          1161: logarithms.
        !          1162: 
        !          1163:      log(y(t)) = log k + pt
        !          1164: 
        !          1165:                = c  + c t
        !          1166:                   1    2
        !          1167: 
        !          1168: The following segment makes use of the fact that a function of  a
        !          1169: vector is the function applied to the individual components.
        !          1170: 
        !          1171:    X = A(:,1:2);
        !          1172:    c = X\log(y)
        !          1173:    p = exp(sum(c))
        !          1174:    e = norm(exp(X*c)-y,'inf')
        !          1175: 
        !          1176: The resulting output is
        !          1177: 
        !          1178: 
        !          1179: 
        !          1180: 
        !          1181: 
        !          1182: 
        !          1183: 
        !          1184: 
        !          1185: 
        !          1186: 
        !          1187: 
        !          1188: 
        !          1189: 
        !          1190: 
        !          1191: MATLAB, page 18
        !          1192: 
        !          1193: 
        !          1194: 
        !          1195:    C     =
        !          1196: 
        !          1197:       4.9083
        !          1198:       0.5407
        !          1199: 
        !          1200:    P     =
        !          1201: 
        !          1202:     232.5134
        !          1203: 
        !          1204:    E     =
        !          1205: 
        !          1206:       4.9141
        !          1207: 
        !          1208: The   predicted   population   and   maximum   deviation   appear
        !          1209: satisfactory  and  indicate  that  the  exponential  model  is  a
        !          1210: reasonable one to consider.
        !          1211: 
        !          1212:      As a curiousity, we return to  the  degree  six  polynomial.
        !          1213: Since  the coefficient of the high order term is negative and the
        !          1214: value of the polynomial at t = 1 is positive, it must have a root
        !          1215: at some value of  t  greater than one.  The statements
        !          1216: 
        !          1217:    X = A(:,1:7);
        !          1218:    c = X\y;
        !          1219:    c = c(7:-1:1);  //reverse the order of the coefficients
        !          1220:    z = roots(c)
        !          1221: 
        !          1222: produce
        !          1223: 
        !          1224:    Z     =
        !          1225: 
        !          1226:       1.1023-  0.0000*i
        !          1227:       0.3021+  0.7293*i
        !          1228:      -0.8790+  0.6536*i
        !          1229:      -1.2939-  0.0000*i
        !          1230:      -0.8790-  0.6536*i
        !          1231:       0.3021-  0.7293*i
        !          1232: 
        !          1233: There is only one real, positive root.  The corresponding time on
        !          1234: the original scale is
        !          1235: 
        !          1236:    1940 + 40*real(z(1))
        !          1237: 
        !          1238:      =  1984.091
        !          1239: 
        !          1240: We conclude that the United States population should become  zero
        !          1241: early in February of 1984.
        !          1242: 
        !          1243: 
        !          1244: 
        !          1245: 
        !          1246: 
        !          1247: 
        !          1248: 
        !          1249: 
        !          1250: 
        !          1251: 
        !          1252: 
        !          1253: 
        !          1254: 
        !          1255: 
        !          1256: 
        !          1257: MATLAB, page 19
        !          1258: 
        !          1259: 
        !          1260: 
        !          1261: 7.  Partial differential equation example
        !          1262: 
        !          1263: 
        !          1264:      Our second extended example is a boundary value problem  for
        !          1265: Laplace's equation.  The underlying physical problem involves the
        !          1266: conductivity of a  medium  with  cylindrical  inclusions  and  is
        !          1267: considered by Keller and Sachs [7].
        !          1268: 
        !          1269:      Find a function  u(x,y)  satisfying Laplace's equation
        !          1270: 
        !          1271:                u   + u   = 0
        !          1272:                 xx    yy
        !          1273: 
        !          1274: The domain is a unit square with a quarter circle of  radius  rho
        !          1275: removed from one corner.  There are Neumann conditions on the top
        !          1276: and bottom edges and Dirichlet conditions on the remainder of the
        !          1277: boundary.
        !          1278: 
        !          1279: 
        !          1280:                          u  = 0
        !          1281:                           n
        !          1282: 
        !          1283:                      -------------
        !          1284:                     |             .
        !          1285:                     |             .
        !          1286:                     |              .
        !          1287:                     |               .  u = 1
        !          1288:                     |                 .
        !          1289:                     |                    .
        !          1290:                     |                       .
        !          1291:              u = 0  |                        |
        !          1292:                     |                        |
        !          1293:                     |                        |
        !          1294:                     |                        |  u = 1
        !          1295:                     |                        |
        !          1296:                     |                        |
        !          1297:                     |                        |
        !          1298:                      ------------------------
        !          1299: 
        !          1300:                               u  = 0
        !          1301:                                n
        !          1302: 
        !          1303: 
        !          1304: The effective conductivity of an medium  is  then  given  by  the
        !          1305: integral along the left edge,
        !          1306: 
        !          1307:                             1
        !          1308:                  sigma = integral  u (0,y) dy
        !          1309:                            0        n
        !          1310: 
        !          1311: It is of interest to study the relation between  the  radius  rho
        !          1312: and  the  conductivity  sigma.   In particular, as rho approaches
        !          1313: one, sigma becomes infinite.
        !          1314: 
        !          1315: 
        !          1316: 
        !          1317: 
        !          1318: 
        !          1319: 
        !          1320: 
        !          1321: 
        !          1322: 
        !          1323: MATLAB, page 20
        !          1324: 
        !          1325: 
        !          1326: 
        !          1327:      Keller and Sachs use a finite difference approximation.  The
        !          1328: following  technique  makes  use of the fact that the equation is
        !          1329: actually Laplace's equation and leads to a  much  smaller  matrix
        !          1330: problem to solve.
        !          1331: 
        !          1332:      Consider an approximate solution of the form
        !          1333: 
        !          1334:                  n      2j-1
        !          1335:            u =  sum  c r    cos(2j-1)t
        !          1336:                 j=1   j
        !          1337: 
        !          1338: where  r,t  are polar coordinates (t is theta).  The coefficients
        !          1339: are to be determined.  For any set of coefficients, this function
        !          1340: already satisfies the differential  equation  because  the  basis
        !          1341: functions  are  harmonic;  it  satisfies  the  normal  derivative
        !          1342: boundary condition on the bottom edge of the  domain  because  we
        !          1343: used   cos  t   in  preference  to   sin t ; and it satisfies the
        !          1344: boundary condition on the left edge of the domain because we  use
        !          1345: only odd multiples of  t .
        !          1346: 
        !          1347:      The computational task is to find coefficients  so that  the
        !          1348: boundary  conditions on the remaining edges are satisfied as well
        !          1349: as possible.  To accomplish this, pick  m  points  (r,t)  on  the
        !          1350: remaining edges.  It is desirable to have  m > n  and in practice
        !          1351: we usually choose m  to be two or three times as large  as   n  .
        !          1352: Typical  values  of  n  are 10 or 20 and of  m  are 20 to 60.  An
        !          1353: m  by  n  matrix  A  is generated.  The  i,j  element is the j-th
        !          1354: basis  function,  or its normal derivative, evaluated at the i-th
        !          1355: boundary point.  A right hand side with  m   components  is  also
        !          1356: generated.   In this example, the elements of the right hand side
        !          1357: are either zero or one.   The  coefficients  are  then  found  by
        !          1358: solving the overdetermined set of equations
        !          1359: 
        !          1360:             Ac = b
        !          1361: 
        !          1362: in a least squares sense.
        !          1363: 
        !          1364:      Once the coefficients have been determined, the  approximate
        !          1365: solution  is  defined  everywhere  on  the  domain.   It  is then
        !          1366: possible to compute the effective conductivity sigma .  In  fact,
        !          1367: a very simple formula results,
        !          1368: 
        !          1369:                      n       j-1
        !          1370:            sigma =  sum  (-1)   c
        !          1371:                     j=1          j
        !          1372: 
        !          1373:      To use MATLAB for this problem, the following  "program"  is
        !          1374: first  stored  in  the  local computer file system, say under the
        !          1375: name "PDE".
        !          1376: 
        !          1377: 
        !          1378: 
        !          1379: 
        !          1380: 
        !          1381: 
        !          1382: 
        !          1383: 
        !          1384: 
        !          1385: 
        !          1386: 
        !          1387: 
        !          1388: 
        !          1389: MATLAB, page 21
        !          1390: 
        !          1391: 
        !          1392: 
        !          1393: //Conductivity example.
        !          1394: //Parameters ---
        !          1395:    rho       //radius of cylindrical inclusion
        !          1396:    n         //number of terms in solution
        !          1397:    m         //number of boundary points
        !          1398: //initialize operation counter
        !          1399:    flop = <0 0>;
        !          1400: //initialize variables
        !          1401:    m1 = round(m/3);   //number of points on each straight edge
        !          1402:    m2 = m - m1;       //number of points with Dirichlet conditions
        !          1403:    pi = 4*atan(1);
        !          1404: //generate points in Cartesian coordinates
        !          1405:    //right hand edge
        !          1406:    for i = 1:m1, x(i) = 1; y(i) = (1-rho)*(i-1)/(m1-1);
        !          1407:    //top edge
        !          1408:    for i = m2+1:m, x(i) = (1-rho)*(m-i)/(m-m2-1); y(i) = 1;
        !          1409:    //circular edge
        !          1410:    for i = m1+1:m2, t = pi/2*(i-m1)/(m2-m1+1); ...
        !          1411:       x(i) = 1-rho*sin(t);  y(i) = 1-rho*cos(t);
        !          1412: //convert to polar coordinates
        !          1413:    for i = 1:m-1, th(i) = atan(y(i)/x(i));  ...
        !          1414:       r(i) = sqrt(x(i)**2+y(i)**2);
        !          1415:    th(m) = pi/2;  r(m) = 1;
        !          1416: //generate matrix
        !          1417:    //Dirichlet conditions
        !          1418:    for i = 1:m2, for j = 1:n, k = 2*j-1; ...
        !          1419:       a(i,j) = r(i)**k*cos(k*th(i));
        !          1420:    //Neumann conditions
        !          1421:    for i = m2+1:m, for j = 1:n, k = 2*j-1; ...
        !          1422:       a(i,j) = k*r(i)**(k-1)*sin((k-1)*th(i));
        !          1423: //generate right hand side
        !          1424:    for i = 1:m2, b(i) = 1;
        !          1425:    for i = m2+1:m, b(i) = 0;
        !          1426: //solve for coefficients
        !          1427:    c = A\b
        !          1428: //compute effective conductivity
        !          1429:    c(2:2:n) = -c(2:2:n);
        !          1430:    sigma = sum(c)
        !          1431: //output total operation count
        !          1432:    ops = flop(2)
        !          1433: 
        !          1434: 
        !          1435: 
        !          1436: 
        !          1437:      The program can be used within MATLAB by setting  the  three
        !          1438: parameters and then accessing the file.  For example,
        !          1439: 
        !          1440:    rho = .9;
        !          1441:    n = 15;
        !          1442:    m = 30;
        !          1443:    exec('PDE')
        !          1444: 
        !          1445: The resulting output is
        !          1446: 
        !          1447: 
        !          1448: 
        !          1449: 
        !          1450: 
        !          1451: 
        !          1452: 
        !          1453: 
        !          1454: 
        !          1455: MATLAB, page 22
        !          1456: 
        !          1457: 
        !          1458: 
        !          1459:    RHO   =
        !          1460: 
        !          1461:       .9000
        !          1462: 
        !          1463:    N     =
        !          1464: 
        !          1465:     15.
        !          1466: 
        !          1467:    M     =
        !          1468: 
        !          1469:     30.
        !          1470: 
        !          1471:    C     =
        !          1472: 
        !          1473:       2.2275
        !          1474:      -2.2724
        !          1475:       1.1448
        !          1476:       0.1455
        !          1477:      -0.1678
        !          1478:      -0.0005
        !          1479:      -0.3785
        !          1480:       0.2299
        !          1481:       0.3228
        !          1482:      -0.2242
        !          1483:      -0.1311
        !          1484:       0.0924
        !          1485:       0.0310
        !          1486:      -0.0154
        !          1487:      -0.0038
        !          1488: 
        !          1489:    SIGM  =
        !          1490: 
        !          1491:       5.0895
        !          1492: 
        !          1493:    OPS   =
        !          1494: 
        !          1495:       16204.
        !          1496: 
        !          1497: 
        !          1498:      A total of 16204 floating point operations were necessary to
        !          1499: set  up  the  matrix,  solve for the coefficients and compute the
        !          1500: conductivity.  The operation count  is  roughly  proportional  to
        !          1501: m*n**2.   The  results obtained for sigma as a function of rho by
        !          1502: this approach are essentially the same as those obtained  by  the
        !          1503: finite   difference  technique  of  Keller  and  Sachs,  but  the
        !          1504: computational effort involved is much less.
        !          1505: 
        !          1506: 
        !          1507: 
        !          1508: 
        !          1509: 
        !          1510: 
        !          1511: 
        !          1512: 
        !          1513: 
        !          1514: 
        !          1515: 
        !          1516: 
        !          1517: 
        !          1518: 
        !          1519: 
        !          1520: 
        !          1521: MATLAB, page 23
        !          1522: 
        !          1523: 
        !          1524: 
        !          1525: 8.  Eigenvalue sensitivity example
        !          1526: 
        !          1527: 
        !          1528:      In this example, we construct a matrix whose eigenvalues are
        !          1529: moderately  sensitive  to  perturbations  and  then  analyze that
        !          1530: sensitivity. We begin with the statement
        !          1531: 
        !          1532:    B = <3 0 7; 0 2 0; 0 0 1>
        !          1533: 
        !          1534: which produces
        !          1535: 
        !          1536:    B     =
        !          1537: 
        !          1538:        3.    0.    7.
        !          1539:        0.    2.    0.
        !          1540:        0.    0.    1.
        !          1541: 
        !          1542: 
        !          1543:      Obviously, the eigenvalues of B are 1, 2 and 3 .   Moreover,
        !          1544: since   B  is  not  symmetric,  these  eigenvalues  are  slightly
        !          1545: sensitive to perturbation.  (The value b(1,3) = 7 was  chosen  so
        !          1546: that the elements of the matrix A below are less than 1000.)
        !          1547: 
        !          1548:      We now generate a similarity transformation to disguise  the
        !          1549: eigenvalues and make them more sensitive.
        !          1550: 
        !          1551:    L = <1 0 0; 2 1 0; -3 4 1>, M = L\L'
        !          1552: 
        !          1553:    L     =
        !          1554: 
        !          1555:        1.    0.    0.
        !          1556:        2.    1.    0.
        !          1557:       -3.    4.    1.
        !          1558: 
        !          1559:    M     =
        !          1560: 
        !          1561:        1.0000    2.0000   -3.0000
        !          1562:       -2.0000   -3.0000   10.0000
        !          1563:       11.0000   18.0000  -48.0000
        !          1564: 
        !          1565: The matrix M has determinant equal to 1 and is  moderately  badly
        !          1566: conditioned.  The similarity transformation is
        !          1567: 
        !          1568:    A = M*B/M
        !          1569: 
        !          1570:    A     =
        !          1571: 
        !          1572:      -64.0000   82.0000   21.0000
        !          1573:      144.0000 -178.0000  -46.0000
        !          1574:     -771.0000  962.0000  248.0000
        !          1575: 
        !          1576: Because  det(M) = 1 , the elements of  A  would be exact integers
        !          1577: if there were no roundoff.  So,
        !          1578: 
        !          1579: 
        !          1580: 
        !          1581: 
        !          1582: 
        !          1583: 
        !          1584: 
        !          1585: 
        !          1586: 
        !          1587: MATLAB, page 24
        !          1588: 
        !          1589: 
        !          1590: 
        !          1591:    A = round(A)
        !          1592: 
        !          1593:    A     =
        !          1594: 
        !          1595:      -64.   82.   21.
        !          1596:      144. -178.  -46.
        !          1597:     -771.  962.  248.
        !          1598: 
        !          1599: 
        !          1600:      This, then, is our test matrix.  We can now  forget  how  it
        !          1601: was generated and analyze its eigenvalues.
        !          1602: 
        !          1603:    <X,D> = eig(A)
        !          1604: 
        !          1605:    D     =
        !          1606: 
        !          1607:        3.0000    0.0000    0.0000
        !          1608:        0.0000    1.0000    0.0000
        !          1609:        0.0000    0.0000    2.0000
        !          1610: 
        !          1611:    X     =
        !          1612: 
        !          1613:        -.0891    3.4903   41.8091
        !          1614:         .1782   -9.1284  -62.7136
        !          1615:        -.9800   46.4473  376.2818
        !          1616: 
        !          1617: Since A is similar to B, its eigenvalues are also  1,  2  and  3.
        !          1618: They  happen  to  be  computed  in  another  order by the EISPACK
        !          1619: subroutines.  The fact that the  columns  of  X,  which  are  the
        !          1620: eigenvectors,  are  so  far  from  being orthonormal is our first
        !          1621: indication that  the  eigenvalues  are  sensitive.  To  see  this
        !          1622: sensitivity, we display more figures of the computed eigenvalues.
        !          1623: 
        !          1624:    long, diag(D)
        !          1625: 
        !          1626:    ANS   =
        !          1627: 
        !          1628:       2.999999999973599
        !          1629:       1.000000000015625
        !          1630:       2.000000000011505
        !          1631: 
        !          1632: We see that, on this computer, the last five significant  figures
        !          1633: are  contaminated  by  roundoff  error.  A  somewhat  superficial
        !          1634: explanation of this is provided by
        !          1635: 
        !          1636:    short,  cond(X)
        !          1637: 
        !          1638:    ANS   =
        !          1639: 
        !          1640:       3.2216e+05
        !          1641: 
        !          1642: The condition number of X gives an upper bound for  the  relative
        !          1643: error  in  the  computed  eigenvalues.   However,  this condition
        !          1644: 
        !          1645: 
        !          1646: 
        !          1647: 
        !          1648: 
        !          1649: 
        !          1650: 
        !          1651: 
        !          1652: 
        !          1653: MATLAB, page 25
        !          1654: 
        !          1655: 
        !          1656: 
        !          1657: number is affected by scaling.
        !          1658: 
        !          1659:    X = X/diag(X(3,:)),  cond(X)
        !          1660: 
        !          1661:    X     =
        !          1662: 
        !          1663:         .0909     .0751     .1111
        !          1664:        -.1818    -.1965    -.1667
        !          1665:        1.0000    1.0000    1.0000
        !          1666: 
        !          1667:    ANS   =
        !          1668: 
        !          1669:       1.7692e+03
        !          1670: 
        !          1671: 
        !          1672:      Rescaling the eigenvectors so that their last components are
        !          1673: all  equal  to  one  has  two consequences. The condition of X is
        !          1674: decreased by over two orders of magnitude.  (This  is  about  the
        !          1675: minimum condition that can be obtained by such diagonal scaling.)
        !          1676: Moreover, it is now apparent  that  the  three  eigenvectors  are
        !          1677: nearly parallel.
        !          1678: 
        !          1679:      More  detailed  information  on  the  sensitivity   of   the
        !          1680: individual eigenvalues involves the left eigenvectors.
        !          1681: 
        !          1682:    Y = inv(X'),  Y'*A*X
        !          1683: 
        !          1684:    Y     =
        !          1685: 
        !          1686:     -511.5000  259.5000  252.0000
        !          1687:      616.0000 -346.0000 -270.0000
        !          1688:      159.5000  -86.5000  -72.0000
        !          1689: 
        !          1690:    ANS   =
        !          1691: 
        !          1692:        3.0000     .0000     .0000
        !          1693:         .0000    1.0000     .0000
        !          1694:         .0000     .0000    2.0000
        !          1695: 
        !          1696: We are now in a position to  compute  the  sensitivities  of  the
        !          1697: individual eigenvalues.
        !          1698: 
        !          1699:    for j = 1:3, c(j) = norm(Y(:,j))*norm(X(:,j)); end,  C
        !          1700: 
        !          1701:    C     =
        !          1702: 
        !          1703:      833.1092
        !          1704:      450.7228
        !          1705:      383.7564
        !          1706: 
        !          1707: These three numbers are the reciprocals of  the  cosines  of  the
        !          1708: angles  between the left and right eigenvectors.  It can be shown
        !          1709: that  perturbation  of  the  elements  of  A  can  result  in   a
        !          1710: 
        !          1711: 
        !          1712: 
        !          1713: 
        !          1714: 
        !          1715: 
        !          1716: 
        !          1717: 
        !          1718: 
        !          1719: MATLAB, page 26
        !          1720: 
        !          1721: 
        !          1722: 
        !          1723: perturbation of the j-th eigenvalue which is c(j) times as large.
        !          1724: In  this  example,  the  first   eigenvalue   has   the   largest
        !          1725: sensitivity.
        !          1726: 
        !          1727:      We now proceed to show that A is close to a  matrix  with  a
        !          1728: double eigenvalue.  The direction of the required perturbation is
        !          1729: given by
        !          1730: 
        !          1731:    E = -1.e-6*Y(:,1)*X(:,1)'
        !          1732: 
        !          1733:    E     =
        !          1734: 
        !          1735:       1.0e-03 *
        !          1736: 
        !          1737:         .0465    -.0930     .5115
        !          1738:        -.0560     .1120    -.6160
        !          1739:        -.0145     .0290    -.1595
        !          1740: 
        !          1741: With some trial and error which we do not show,  we  bracket  the
        !          1742: point  where  two  eigenvalues of a perturbed A coalesce and then
        !          1743: become complex.
        !          1744: 
        !          1745:    eig(A + .4*E),  eig(A + .5*E)
        !          1746: 
        !          1747:    ANS   =
        !          1748: 
        !          1749:        1.1500
        !          1750:        2.5996
        !          1751:        2.2504
        !          1752: 
        !          1753:    ANS   =
        !          1754: 
        !          1755:       2.4067 +  .1753*i
        !          1756:       2.4067 -  .1753*i
        !          1757:       1.1866 + 0.0000*i
        !          1758: 
        !          1759: Now, a bisecting search, driven by the imaginary part of  one  of
        !          1760: the eigenvalues, finds the point where two eigenvalues are nearly
        !          1761: equal.
        !          1762: 
        !          1763:    r = .4;  s = .5;
        !          1764: 
        !          1765:    while s-r > 1.e-14, t = (r+s)/2; d = eig(A+t*E); ...
        !          1766:      if imag(d(1))=0, r = t; else, s = t;
        !          1767: 
        !          1768:    long,  T
        !          1769: 
        !          1770:    T     =
        !          1771: 
        !          1772:         .450380734134507
        !          1773: 
        !          1774: 
        !          1775:      Finally, we display the perturbed matrix, which is obviously
        !          1776: 
        !          1777: 
        !          1778: 
        !          1779: 
        !          1780: 
        !          1781: 
        !          1782: 
        !          1783: 
        !          1784: 
        !          1785: MATLAB, page 27
        !          1786: 
        !          1787: 
        !          1788: 
        !          1789: close  to the original, and its pair of nearly equal eigenvalues.
        !          1790: (We have dropped a few digits from the long output.)
        !          1791: 
        !          1792:    A+t*E,  eig(A+t*E)
        !          1793: 
        !          1794:    A
        !          1795: 
        !          1796:     -63.999979057   81.999958114   21.000230369
        !          1797:     143.999974778 -177.999949557  -46.000277434
        !          1798:    -771.000006530  962.000013061  247.999928164
        !          1799: 
        !          1800:    ANS   =
        !          1801: 
        !          1802:       2.415741150
        !          1803:       2.415740621
        !          1804:       1.168517777
        !          1805: 
        !          1806: 
        !          1807:      The  first  two  eigenvectors  of  A  +   t*E   are   almost
        !          1808: indistinguishable  indicating that the perturbed matrix is almost
        !          1809: defective.
        !          1810: 
        !          1811:    <X,D> = eig(A+t*E);  X = X/diag(X(3,:))
        !          1812: 
        !          1813:    X     =
        !          1814: 
        !          1815:        .096019578     .096019586    .071608466
        !          1816:       -.178329614    -.178329608   -.199190520
        !          1817:       1.000000000    1.000000000   1.000000000
        !          1818: 
        !          1819:    short,  cond(X)
        !          1820: 
        !          1821:    ANS   =
        !          1822: 
        !          1823:       3.3997e+09
        !          1824: 
        !          1825: 
        !          1826: 9.  Syntax diagrams
        !          1827: 
        !          1828: 
        !          1829:      A formal description of the language acceptable  to  MATLAB,
        !          1830: as well as a flow chart of the MATLAB program, is provided by the
        !          1831: syntax diagrams or syntax graphs of Wirth [6].  There are  eleven
        !          1832: non-terminal symbols in the language:
        !          1833: 
        !          1834:    line, statement, clause, expression, term,
        !          1835:    factor, number, integer, name, command, text .
        !          1836: 
        !          1837: The diagrams define each of the non-terminal  symbols  using  the
        !          1838: others and the terminal symbols:
        !          1839: 
        !          1840:    letter -- A through Z,
        !          1841:    digit  -- 0 through 9,
        !          1842: 
        !          1843: 
        !          1844: 
        !          1845: 
        !          1846: 
        !          1847: 
        !          1848: 
        !          1849: 
        !          1850: 
        !          1851: MATLAB, page 28
        !          1852: 
        !          1853: 
        !          1854: 
        !          1855:    char   -- ( ) ; : + - * / \ = . , < >
        !          1856:    quote  -- '
        !          1857: 
        !          1858: 
        !          1859: line
        !          1860: 
        !          1861:        |-----> statement >----|
        !          1862:        |                      |
        !          1863:        |-----> clause >-------|
        !          1864:        |                      |
        !          1865: -------|-----> expr >---------|------>
        !          1866:      | |                      | |
        !          1867:      | |-----> command >------| |
        !          1868:      | |                      | |
        !          1869:      | |-> > >-> expr >-> < >-| |
        !          1870:      | |                      | |
        !          1871:      | |----------------------| |
        !          1872:      |                          |
        !          1873:      |        |-< ; <-|         |
        !          1874:      |--------|       |---------|
        !          1875:               |-< , <-|
        !          1876: 
        !          1877: 
        !          1878: 
        !          1879: 
        !          1880: statement
        !          1881: 
        !          1882:      |-> name >--------------------------------|
        !          1883:      |          |                              |
        !          1884:      |          |         |--> : >---|         |
        !          1885:      |          |         |          |         |
        !          1886:      |          |-> ( >---|-> expr >-|---> ) >-|
        !          1887:      |                  |              |       |
        !          1888: -----|                  |-----< , <----|       |--> = >--> expr >--->
        !          1889:      |                                         |
        !          1890:      |       |--< , <---|                      |
        !          1891:      |       |          |                      |
        !          1892:      |-> < >---> name >---> > >----------------|
        !          1893: 
        !          1894: 
        !          1895: 
        !          1896: 
        !          1897: 
        !          1898: 
        !          1899: 
        !          1900: 
        !          1901: 
        !          1902: 
        !          1903: 
        !          1904: 
        !          1905: 
        !          1906: 
        !          1907: 
        !          1908: 
        !          1909: 
        !          1910: 
        !          1911: 
        !          1912: 
        !          1913: 
        !          1914: 
        !          1915: 
        !          1916: 
        !          1917: MATLAB, page 29
        !          1918: 
        !          1919: 
        !          1920: 
        !          1921: clause
        !          1922: 
        !          1923:      |---> FOR   >---> name >---> = >---> expr >--------------|
        !          1924:      |                                                        |
        !          1925:      | |-> WHILE >-|                                          |
        !          1926:      |-|           |-> expr >----------------------           |
        !          1927:      | |-> IF    >-|          |   |   |   |   |   |           |
        !          1928: -----|                        <   <=  =   <>  >=  >           |---->
        !          1929:      |                        |   |   |   |   |   |           |
        !          1930:      |                        ----------------------> expr >--|
        !          1931:      |                                                        |
        !          1932:      |---> ELSE  >--------------------------------------------|
        !          1933:      |                                                        |
        !          1934:      |---> END   >--------------------------------------------|
        !          1935: 
        !          1936: 
        !          1937: 
        !          1938: 
        !          1939: expr
        !          1940: 
        !          1941:        |-> + >-|
        !          1942:        |       |
        !          1943: -------|-------|-------> term >---------->
        !          1944:        |       |    |             |
        !          1945:        |-> - >-|    |  |-< + <-|  |
        !          1946:                     |  |       |  |
        !          1947:                     |--|-< - <-|--|
        !          1948:                        |       |
        !          1949:                        |-< : <-|
        !          1950: 
        !          1951: 
        !          1952: 
        !          1953: 
        !          1954: term
        !          1955: 
        !          1956: ---------------------> factor >---------------------->
        !          1957:         |                                   |
        !          1958:         |             |-< * <-|             |
        !          1959:         |  |-------|  |       |  |-------|  |
        !          1960:         |--|       |--|-< / <-|--|       |--|
        !          1961:            |-< . <-|  |       |  |-< . <-|
        !          1962:                       |-< \ <-|
        !          1963: 
        !          1964: 
        !          1965: 
        !          1966: 
        !          1967: 
        !          1968: 
        !          1969: 
        !          1970: 
        !          1971: 
        !          1972: 
        !          1973: 
        !          1974: 
        !          1975: 
        !          1976: 
        !          1977: 
        !          1978: 
        !          1979: 
        !          1980: 
        !          1981: 
        !          1982: 
        !          1983: MATLAB, page 30
        !          1984: 
        !          1985: 
        !          1986: 
        !          1987: factor
        !          1988: 
        !          1989:      |----------------> number >---------------|
        !          1990:      |                                         |
        !          1991:      |-> name >--------------------------------|
        !          1992:      |          |                              |
        !          1993:      |          |         |--> : >---|         |
        !          1994:      |          |         |          |         |
        !          1995:      |          |-> ( >---|-> expr >-|---> ) >-|
        !          1996:      |                  |              |       |
        !          1997:      |                  |-----< , <----|       |
        !          1998:      |                                         |
        !          1999: -----|------------> ( >-----> expr >-----> ) >-|-|-------|----->
        !          2000:      |                                         | |       | |
        !          2001:      |                  |--------------|       | |-> ' >-| |
        !          2002:      |                  |              |       |           |
        !          2003:      |------------> < >-|---> expr >---|-> > >-|           |
        !          2004:      |                    |          |         |           |
        !          2005:      |                    |--<   <---|         |           |
        !          2006:      |                    |          |         |           |
        !          2007:      |                    |--< ; <---|         |           |
        !          2008:      |                    |          |         |           |
        !          2009:      |                    |--< , <---|         |           |
        !          2010:      |                                         |           |
        !          2011:      |------------> > >-----> expr >-----> < >-|           |
        !          2012:      |                                         |           |
        !          2013:      |-----> factor >---> ** >---> factor >----|           |
        !          2014:      |                                                     |
        !          2015:      |------------> ' >-----> text >-----> ' >-------------|
        !          2016: 
        !          2017: 
        !          2018: 
        !          2019: 
        !          2020: number
        !          2021: 
        !          2022:     |----------|                          |-> + >-|
        !          2023:     |          |                          |       |
        !          2024: -----> int >-----> . >---> int >-----> E >---------> int >---->
        !          2025:              |                   | |      |       |        |
        !          2026:              |                   | |      |-> - >-|        |
        !          2027:              |                   | |                       |
        !          2028:              |---------------------------------------------|
        !          2029: 
        !          2030: 
        !          2031: 
        !          2032: 
        !          2033: int
        !          2034: 
        !          2035: ------------> digit >----------->
        !          2036:           |           |
        !          2037:           |-----------|
        !          2038: 
        !          2039: 
        !          2040: 
        !          2041: 
        !          2042: 
        !          2043: 
        !          2044: 
        !          2045: 
        !          2046: 
        !          2047: 
        !          2048: 
        !          2049: MATLAB, page 31
        !          2050: 
        !          2051: 
        !          2052: 
        !          2053: 
        !          2054: 
        !          2055: name
        !          2056: 
        !          2057:                   |--< letter <--|
        !          2058:                   |              |
        !          2059: ------> letter >--|--------------|----->
        !          2060:                   |              |
        !          2061:                   |--< digit  <--|
        !          2062: 
        !          2063: 
        !          2064: 
        !          2065: 
        !          2066: command
        !          2067: 
        !          2068:                         |--> name >--|
        !          2069:                         |            |
        !          2070: --------> name >--------|------------|---->
        !          2071:                         |            |
        !          2072:                         |--> char >--|
        !          2073:                         |            |
        !          2074:                         |---> ' >----|
        !          2075: 
        !          2076: text
        !          2077: 
        !          2078:                 |-> letter >--|
        !          2079:                 |             |
        !          2080:                 |-> digit >---|
        !          2081: ----------------|             |-------------->
        !          2082:             |   |-> char >----|   |
        !          2083:             |   |             |   |
        !          2084:             |   |-> ' >-> ' >-|   |
        !          2085:             |                     |
        !          2086:             |---------------------|
        !          2087: 
        !          2088: 
        !          2089: 10.  The parser-interpreter
        !          2090: 
        !          2091: 
        !          2092:      The structure of the parser-interpreter is similar  to  that
        !          2093: of  Wirth's  compiler  [6] for his simple language, PL/0 , except
        !          2094: that MATLAB  is  programmed  in  Fortran,  which  does  not  have
        !          2095: explicit recursion.  The interrelation of the primary subroutines
        !          2096: is shown in the following diagram.
        !          2097: 
        !          2098: 
        !          2099: 
        !          2100: 
        !          2101: 
        !          2102: 
        !          2103: 
        !          2104: 
        !          2105: 
        !          2106: 
        !          2107: 
        !          2108: 
        !          2109: 
        !          2110: 
        !          2111: 
        !          2112: 
        !          2113: 
        !          2114: 
        !          2115: MATLAB, page 32
        !          2116: 
        !          2117: 
        !          2118: 
        !          2119:       MAIN
        !          2120:         |
        !          2121:       MATLAB    |--CLAUSE
        !          2122:         |       |    |
        !          2123:       PARSE-----|--EXPR----TERM----FACTOR
        !          2124:                 |    |       |       |
        !          2125:                 |    |-------|-------|
        !          2126:                 |    |       |       |
        !          2127:                 |  STACK1  STACK2  STACKG
        !          2128:                 |
        !          2129:                 |--STACKP--PRINT
        !          2130:                 |
        !          2131:                 |--COMAND
        !          2132:                 |
        !          2133:                 |
        !          2134:                 |          |--CGECO
        !          2135:                 |          |
        !          2136:                 |          |--CGEFA
        !          2137:                 |          |
        !          2138:                 |--MATFN1--|--CGESL
        !          2139:                 |          |
        !          2140:                 |          |--CGEDI
        !          2141:                 |          |
        !          2142:                 |          |--CPOFA
        !          2143:                 |
        !          2144:                 |
        !          2145:                 |          |--IMTQL2
        !          2146:                 |          |
        !          2147:                 |          |--HTRIDI
        !          2148:                 |          |
        !          2149:                 |--MATFN2--|--HTRIBK
        !          2150:                 |          |
        !          2151:                 |          |--CORTH
        !          2152:                 |          |
        !          2153:                 |          |--COMQR3
        !          2154:                 |
        !          2155:                 |
        !          2156:                 |--MATFN3-----CSVDC
        !          2157:                 |
        !          2158:                 |
        !          2159:                 |          |--CQRDC
        !          2160:                 |--MATFN4--|
        !          2161:                 |          |--CQRSL
        !          2162:                 |
        !          2163:                 |
        !          2164:                 |          |--FILES
        !          2165:                 |--MATFN5--|
        !          2166:                            |--SAVLOD
        !          2167: 
        !          2168:      Subroutine  PARSE  controls  the  interpretation   of   each
        !          2169: statement.    It  calls  subroutines  that  process  the  various
        !          2170: syntactic  quantities  such  as  command,  expression,  term  and
        !          2171: factor.   A  fairly  simple  program stack mechanism allows these
        !          2172: 
        !          2173: 
        !          2174: 
        !          2175: 
        !          2176: 
        !          2177: 
        !          2178: 
        !          2179: 
        !          2180: 
        !          2181: MATLAB, page 33
        !          2182: 
        !          2183: 
        !          2184: 
        !          2185: subroutines to recursively "call"  each  other  along  the  lines
        !          2186: allowed  by  the  syntax  diagrams.   The  four STACK subroutines
        !          2187: manage the variable memory  and  perform  elementary  operations,
        !          2188: such as matrix addition and transposition.
        !          2189: 
        !          2190:      The  four  subroutines  MATFN1  though  MATFN4  are   called
        !          2191: whenever  "serious"  matrix  computations are required.  They are
        !          2192: interface routines which call the  various  LINPACK  and  EISPACK
        !          2193: subroutines.  MATFN5 primarily handles the file access tasks.
        !          2194: 
        !          2195:      Two large real arrays, STKR and STKI, are used to store  all
        !          2196: the  matrices.   Four integer arrays are used to store the names,
        !          2197: the row and column dimensions, and the  pointers  into  the  real
        !          2198: stacks.  The following diagram illustrates this storage scheme.
        !          2199: 
        !          2200: TOP         IDSTK     MSTK NSTK LSTK               STKR       STKI
        !          2201:  --      -- -- -- --   --   --   --              --------   --------
        !          2202: |  |--->|  |  |  |  | |  | |  | |  |----------->|        | |        |
        !          2203:  --      -- -- -- --   --   --   --              --------   --------
        !          2204:         |  |  |  |  | |  | |  | |  |            |        | |        |
        !          2205:          -- -- -- --   --   --   --              --------   --------
        !          2206:               .         .    .    .                  .          .
        !          2207:               .         .    .    .                  .          .
        !          2208:               .         .    .    .                  .          .
        !          2209:          -- -- -- --   --   --   --              --------   --------
        !          2210: BOT     |  |  |  |  | |  | |  | |  |            |        | |        |
        !          2211:  --      -- -- -- --   --   --   --              --------   --------
        !          2212: |  |--->| X|  |  |  | | 2| | 1| |  |----------->|  3.14  | |  0.00  |
        !          2213:  --      -- -- -- --   --   --   --              --------   --------
        !          2214:         | A|  |  |  | | 2| | 2| |  |---------   |  0.00  | |  1.00  |
        !          2215:          -- -- -- --   --   --   --          \   --------   --------
        !          2216:         | E| P| S|  | | 1| | 1| |  |-------   ->| 11.00  | |  0.00  |
        !          2217:          -- -- -- --   --   --   --        \     --------   --------
        !          2218:         | F| L| O| P| | 1| | 2| |  |------  \   | 21.00  | |  0.00  |
        !          2219:          -- -- -- --   --   --   --       \  \   --------   --------
        !          2220:         | E| Y| E|  | |-1| |-1| |  |---    \ |  | 12.00  | |  0.00  |
        !          2221:          -- -- -- --   --   --   --    \   | |   --------   --------
        !          2222:         | R| A| N| D| | 1| | 1| |  |-   \  | |  | 22.00  | |  0.00  |
        !          2223:          -- -- -- --   --   --   --  \  |  \ \   --------   --------
        !          2224:                                      |  \   \ ->| 1.E-15 | |  0.00  |
        !          2225:                                      \   \   \   --------   --------
        !          2226:                                       \   \   ->|  0.00  | |  0.00  |
        !          2227:                                        \   \     --------   --------
        !          2228:                                         \   \   |  0.00  | |  0.00  |
        !          2229:                                          \   \   --------   --------
        !          2230:                                           \   ->|  1.00  | |  0.00  |
        !          2231:                                            \     --------   --------
        !          2232:                                             --->| URAND  | |  0.00  |
        !          2233:                                                  --------   --------
        !          2234: 
        !          2235:      The top portion of the stack is used for temporary variables
        !          2236: and the bottom portion for saved variables.  The figure shows the
        !          2237: situation after the line
        !          2238: 
        !          2239: 
        !          2240: 
        !          2241: 
        !          2242: 
        !          2243: 
        !          2244: 
        !          2245: 
        !          2246: 
        !          2247: MATLAB, page 34
        !          2248: 
        !          2249: 
        !          2250: 
        !          2251:    A = <11,12; 21,22>,  x = <3.14, sqrt(-1)>'
        !          2252: 
        !          2253: has been processed.  The four permanent names,  EPS,  FLOP,  RAND
        !          2254: and  EYE,  occupy the last four positions of the variable stacks.
        !          2255: RAND has dimensions 1 by 1, but whenever its value is  requested,
        !          2256: a random number generator is used instead.  EYE has dimensions -1
        !          2257: by -1 to indicate that the actual dimensions must  be  determined
        !          2258: later by context.  The two saved variables have dimensions 2 by 2
        !          2259: and 2 by 1 and so take up a total of 6 locations.
        !          2260: 
        !          2261:      Subsequent statements involving  A  and  x  will  result  in
        !          2262: temporary  copies  being  made in the top of the stack for use in
        !          2263: the actual calculations.  Whenever the top of the  stack  reaches
        !          2264: the  bottom,  a  message  indicating  memory has been exceeded is
        !          2265: printed, but the current variables are not affected.
        !          2266: 
        !          2267:      This modular structure makes it possible to implement MATLAB
        !          2268: on a system with a limited amount of memory.  The object code for
        !          2269: the MATFN's and the LINPACK-EISPACK subroutines is rarely needed.
        !          2270: Although  it  is  not  standard,  many  Fortran operating systems
        !          2271: provide some overlay mechanism so that this code is brought  into
        !          2272: the  main memory only when required.  The variables, which occupy
        !          2273: a relatively small portion of the memory, remain in place,  while
        !          2274: the subroutines which process them are loaded a few at a time.
        !          2275: 
        !          2276: 
        !          2277: 11.  The numerical algorithms
        !          2278: 
        !          2279: 
        !          2280:      The algorithms underlying the  basic  MATLAB  functions  are
        !          2281: described  in the LINPACK and EISPACK guides [1-3]. The following
        !          2282: list gives the subroutines used by these functions.
        !          2283: 
        !          2284:    INV(A)          - CGECO,CGEDI
        !          2285:    DET(A)          - CGECO,CGEDI
        !          2286:    LU(A)           - CGEFA
        !          2287:    RCOND(A)        - CGECO
        !          2288:    CHOL(A)         - CPOFA
        !          2289:    SVD(A)          - CSVDC
        !          2290:    COND(A)         - CSVDC
        !          2291:    NORM(A,2)       - CSVDC
        !          2292:    PINV(A,eps)     - CSVDC
        !          2293:    RANK(A,eps)     - CSVDC
        !          2294:    QR(A)           - CQRDC,CQRSL
        !          2295:    ORTH(A)         - CQRDC,CSQSL
        !          2296:    A\B and B/A     - CGECO,CGESL if A is square.
        !          2297:                    - CQRDC,CQRSL if A is not square.
        !          2298:    EIG(A)          - HTRIDI,IMTQL2,HTRIBK if A is Hermitian.
        !          2299:                    - CORTH,COMQR2         if A is not Hermitian.
        !          2300:    SCHUR(A)        - same as EIG.
        !          2301:    HESS(A)         - same as EIG.
        !          2302: 
        !          2303: 
        !          2304: 
        !          2305: 
        !          2306: 
        !          2307: 
        !          2308: 
        !          2309: 
        !          2310: 
        !          2311: 
        !          2312: 
        !          2313: MATLAB, page 35
        !          2314: 
        !          2315: 
        !          2316: 
        !          2317:      Minor modifications were made to all these subroutines.  The
        !          2318: LINPACK  routines  were  changed  to  replace the Fortran complex
        !          2319: arithmetic with explicit references to real and imaginary  parts.
        !          2320: Since  most of the floating point arithmetic is concentrated in a
        !          2321: few low-level subroutines which perform  vector  operations  (the
        !          2322: Basic  Linear  Algebra  Subprograms),  this  was not an extensive
        !          2323: change.  It also facilitated implementation of the FLOP and  CHOP
        !          2324: features  which count and optionally truncate each floating point
        !          2325: operation.
        !          2326: 
        !          2327:      The EISPACK subroutine COMQR2 was modified to  allow  access
        !          2328: to  the  Schur  triangular  form, ordinarily just an intermediate
        !          2329: result.   IMTQL2  was  modified  to  make  computation   of   the
        !          2330: eigenvectors   optional.    Both  subroutines  were  modified  to
        !          2331: eliminate the machine-dependent accuracy parameter  and  all  the
        !          2332: EISPACK subroutines were changed to include FLOP and CHOP.
        !          2333: 
        !          2334:      The algorithms employed for the  POLY  and  ROOTS  functions
        !          2335: illustrate  an  interesting  aspect  of  the  modern  approach to
        !          2336: eigenvalue computation.   POLY(A)  generates  the  characteristic
        !          2337: polynomial  of  A  and  ROOTS(POLY(A))  finds  the  roots of that
        !          2338: polynomial, which are, of course, the eigenvalues of A . But both
        !          2339: POLY  and  ROOTS  use  EISPACK eigenvalues subroutines, which are
        !          2340: based on similarity transformations.  So the  classical  approach
        !          2341: which  characterizes  eigenvalues  as roots of the characteristic
        !          2342: polynomial is actually reversed.
        !          2343: 
        !          2344:      If A is an n by n matrix, POLY(A) produces the  coefficients
        !          2345: C(1) through C(n+1), with C(1) = 1, in
        !          2346: 
        !          2347:       DET(z*EYE-A) = C(1)*z**n + ... + C(n)*z + C(n+1) .
        !          2348: 
        !          2349: The algorithm can be expressed compactly using MATLAB:
        !          2350: 
        !          2351:       Z = EIG(A);
        !          2352:       C = 0*ONES(n+1,1);  C(1) = 1;
        !          2353:       for j = 1:n, C(2:j+1) = C(2:j+1) - Z(j)*C(1:j);
        !          2354:       C
        !          2355: 
        !          2356: This recursion is easily derived by expanding the product
        !          2357: 
        !          2358:       (z - z(1))*(z - z(2))* ... * (z-z(n)) .
        !          2359: 
        !          2360: It is possible to prove that POLY(A) produces the coefficients in
        !          2361: the  characteristic  polynomial of a matrix within roundoff error
        !          2362: of  A .  This is true even if the  eigenvalues  of  A  are  badly
        !          2363: conditioned.    The  traditional  algorithms  for  obtaining  the
        !          2364: characteristic polynomial which do not use the eigenvalues do not
        !          2365: have such satisfactory numerical properties.
        !          2366: 
        !          2367:      If C is a vector with n+1  components,  ROOTS(C)  finds  the
        !          2368: roots of the polynomial of degree n ,
        !          2369: 
        !          2370: 
        !          2371: 
        !          2372: 
        !          2373: 
        !          2374: 
        !          2375: 
        !          2376: 
        !          2377: 
        !          2378: 
        !          2379: MATLAB, page 36
        !          2380: 
        !          2381: 
        !          2382: 
        !          2383:        p(z) = C(1)*z**n + ... + C(n)*z + C(n+1) .
        !          2384: 
        !          2385: The algorithm simply involves computing the  eigenvalues  of  the
        !          2386: companion matrix:
        !          2387: 
        !          2388:       A = 0*ONES(n,n)
        !          2389:       for j = 1:n, A(1,j) = -C(j+1)/C(1);
        !          2390:       for i = 2:n, A(i,i-1) = 1;
        !          2391:       EIG(A)
        !          2392: 
        !          2393: It is possible to prove that the results produced are  the  exact
        !          2394: eigenvalues  of  a  matrix within roundoff error of the companion
        !          2395: matrix A, but this does not mean that they are the exact roots of
        !          2396: a  polynomial with coefficients within roundoff error of those in
        !          2397: C .  There are more accurate, more efficient methods for  finding
        !          2398: polynomial  roots,  but  this  approach has the crucial advantage
        !          2399: that it does not require very much additional code.
        !          2400: 
        !          2401:      The elementary functions EXP, LOG, SQRT, SIN, COS  and  ATAN
        !          2402: are  applied  to  square  matrices  by  diagonalizing the matrix,
        !          2403: applying the functions to the  individual  eigenvalues  and  then
        !          2404: transforming back.  For example, EXP(A) is computed by
        !          2405: 
        !          2406:       <X,D> = EIG(A);
        !          2407:       for j = 1:n, D(j,j) = EXP(D(j,j));
        !          2408:       X*D/X
        !          2409: 
        !          2410: This is essentially method number 14  out  of  the  19  'dubious'
        !          2411: possibilities described in [8].  It is dubious because it doesn't
        !          2412: always work.  The matrix of eigenvectors  X  can  be  arbitrarily
        !          2413: badly  conditioned  and  all  accuracy lost in the computation of
        !          2414: X*D/X.  A warning message is printed if RCOND(X) is  very  small,
        !          2415: but  this  only  catches the extreme cases.  An example of a case
        !          2416: not detected is when A has a double eigenvalue, but theoretically
        !          2417: only  one  linearly  independent  eigenvector associated with it.
        !          2418: The computed eigenvalues will be separated by  something  on  the
        !          2419: order  of the square root of the roundoff level.  This separation
        !          2420: will be reflected in RCOND(X) which will probably  not  be  small
        !          2421: enough to trigger the error message.  The computed EXP(A) will be
        !          2422: accurate to only half precision.  Better methods  are  known  for
        !          2423: computing EXP(A), but they do not easily extend to the other five
        !          2424: functions and would require a considerable amount  of  additional
        !          2425: code.
        !          2426: 
        !          2427:      The expression A**p is evaluated by repeated  multiplication
        !          2428: if p is an integer greater than 1.  Otherwise it is evaluated by
        !          2429: 
        !          2430:       <X,D> = EIG(A);
        !          2431:       for j = 1:n, D(j,j) = EXP(p*LOG(D(j,j)))
        !          2432:       X*D/X
        !          2433: 
        !          2434: This suffers from the same potential loss of  accuracy  if  X  is
        !          2435: badly conditioned.  It was partly for this reason that the case p
        !          2436: 
        !          2437: 
        !          2438: 
        !          2439: 
        !          2440: 
        !          2441: 
        !          2442: 
        !          2443: 
        !          2444: 
        !          2445: MATLAB, page 37
        !          2446: 
        !          2447: 
        !          2448: 
        !          2449: = 1 is included in the general case.  Comparison of A**1  with  A
        !          2450: gives some idea of the loss of accuracy for other values of p and
        !          2451: for the elementary functions.
        !          2452: 
        !          2453:      RREF, the reduced row echelon form, is of some  interest  in
        !          2454: theoretical  linear algebra, although it has little computational
        !          2455: value.  It is included in MATLAB for  pedagogical  reasons.   The
        !          2456: algorithm  is essentially Gauss-Jordan elimination with detection
        !          2457: of negligible columns applied to rectangular matrices.
        !          2458: 
        !          2459:      There are three separate places in MATLAB where the rank  of
        !          2460: a  matrix  is  implicitly  computed:  in RREF(A), in A\B for non-
        !          2461: square A, and in  the  pseudoinverse  PINV(A).   Three  different
        !          2462: algorithms  with  three  different criteria for negligibility are
        !          2463: used and so it is possible that three different values  could  be
        !          2464: produced for the same matrix.  With RREF(A), the rank of A is the
        !          2465: number of nonzero rows.  The elimination algorithm used for  RREF
        !          2466: is  the  fastest of the three rank-determining algorithms, but it
        !          2467: is the least sophisticated numerically and  the  least  reliable.
        !          2468: With  A\B,  the  algorithm  is  essentially  that used by example
        !          2469: subroutine SQRST  in  chapter  9  of  the  LINPACK  guide.   With
        !          2470: PINV(A),   the   algorithm   is   based  on  the  singular  value
        !          2471: decomposition and is described  in  chapter  11  of  the  LINPACK
        !          2472: guide.   The  SVD  algorithm  is the most time-consuming, but the
        !          2473: most reliable and is therefore also used for RANK(A).
        !          2474: 
        !          2475:      The  uniformly  distributed  random  numbers  in  RAND   are
        !          2476: obtained  from  the  machine-independent  random number generator
        !          2477: URAND described in [9].  It is possible  to  switch  to  normally
        !          2478: distributed   random   numbers,   which   are  obtained  using  a
        !          2479: transformation also described in [9].
        !          2480: 
        !          2481:      The computation of
        !          2482: 
        !          2483:                 2    2
        !          2484:           sqrt(a  + b )
        !          2485: 
        !          2486: is  required  in  many  matrix  algorithms,  particularly   those
        !          2487: involving  complex  arithmetic.   A  new approach to carrying out
        !          2488: this operation is described by Moler and Morrison [10].  It is  a
        !          2489: cubically  convergent  algorithm  which  starts with  a  and  b ,
        !          2490: rather than with their squares, and  thereby  avoids  destructive
        !          2491: arithmetic underflows and overflows.  In MATLAB, the algorithm is
        !          2492: used for complex modulus, Euclidean vector norm, plane rotations,
        !          2493: and  the  shift  calculation in the eigenvalue and singular value
        !          2494: iterations.
        !          2495: 
        !          2496: 
        !          2497: 12.  FLOP and CHOP
        !          2498: 
        !          2499:      Detailed information about the amount of  work  involved  in
        !          2500: matrix  calculations  and  the  resulting accuracy is provided by
        !          2501: FLOP and CHOP.  The basic unit of work is the "flop", or floating
        !          2502: 
        !          2503: 
        !          2504: 
        !          2505: 
        !          2506: 
        !          2507: 
        !          2508: 
        !          2509: 
        !          2510: 
        !          2511: MATLAB, page 38
        !          2512: 
        !          2513: 
        !          2514: 
        !          2515: point operation.  Roughly, one flop is one execution of a Fortran
        !          2516: statement like
        !          2517: 
        !          2518:       S = S + X(I)*Y(I)
        !          2519: 
        !          2520: or
        !          2521: 
        !          2522:       Y(I) = Y(I) + T*X(I)
        !          2523: 
        !          2524: In other words, it consists of one floating point multiplication,
        !          2525: together  with  one  floating  point  addition and the associated
        !          2526: indexing and storage reference operations.
        !          2527: 
        !          2528:      MATLAB will  print  the  number  of  flops  required  for  a
        !          2529: particular statement when the statement is terminated by an extra
        !          2530: comma.  For example, the line
        !          2531: 
        !          2532:       n = 20;  RAND(n)*RAND(n);,
        !          2533: 
        !          2534: ends with an extra comma.  Two  20  by  20  random  matrices  are
        !          2535: generated  and  multiplied  together.   The result is assigned to
        !          2536: ANS, but the semicolon suppresses its printing.  The only  output
        !          2537: is
        !          2538: 
        !          2539:         8800 flops
        !          2540: 
        !          2541: This is  n**3 + 2*n**2  flops,  n**2  for each random matrix  and
        !          2542: n**3 for the product.
        !          2543: 
        !          2544:      FLOP is a predefined vector with two components.  FLOP(1) is
        !          2545: the number of flops used by the most recently executed statement,
        !          2546: except that statements with zero flops are ignored.  For example,
        !          2547: after executing the previous statement,
        !          2548: 
        !          2549:       flop(1)/n**3
        !          2550: 
        !          2551: results in
        !          2552: 
        !          2553:       ANS   =
        !          2554: 
        !          2555:           1.1000
        !          2556: 
        !          2557: 
        !          2558:      FLOP(2) is the cumulative total of all the flops used  since
        !          2559: the beginning of the MATLAB session.  The statement
        !          2560: 
        !          2561:       FLOP = <0 0>
        !          2562: 
        !          2563: resets the total.
        !          2564: 
        !          2565:      There are several difficulties  associated  with  keeping  a
        !          2566: precise  count  of  floating  point  operations.  An  addition or
        !          2567: subtraction that is not paired with a multiplication  is  usually
        !          2568: 
        !          2569: 
        !          2570: 
        !          2571: 
        !          2572: 
        !          2573: 
        !          2574: 
        !          2575: 
        !          2576: 
        !          2577: MATLAB, page 39
        !          2578: 
        !          2579: 
        !          2580: 
        !          2581: counted as a flop. The same is true of an isolated multiplication
        !          2582: that is  not  paired  with  an  addition.   Each  floating  point
        !          2583: division counts as a flop.  But the number of operations required
        !          2584: by system dependent library functions such as square root  cannot
        !          2585: be  counted, so most elementary functions are arbitrarily counted
        !          2586: as using only one flop.
        !          2587: 
        !          2588:      The  biggest  difficulty  occurs  with  complex  arithmetic.
        !          2589: Almost  all operations on the real parts of matrices are counted.
        !          2590: However, the operations on the  complex  parts  of  matrices  are
        !          2591: counted only when they involve nonzero elements.  This means that
        !          2592: simple operations on nonreal matrices require only about twice as
        !          2593: many  flops as the same operations on real matrices.  This factor
        !          2594: of two is not necessarily an accurate  measure  of  the  relative
        !          2595: costs of real and complex arithmetic.
        !          2596: 
        !          2597:      The result of each floating  point  operation  may  also  be
        !          2598: "chopped" to simulate a computer with a shorter word length.  The
        !          2599: details of this chopping operation depend upon the format of  the
        !          2600: floating point word.  Usually, the fraction in the floating point
        !          2601: word  can  be  regarded  as  consisting  of  several   octal   or
        !          2602: hexadecimal digits.  The least significant of these digits can be
        !          2603: set to zero by a logical masking operation.  Thus the statement
        !          2604: 
        !          2605:       CHOP(p)
        !          2606: 
        !          2607: causes the  p  least significant octal or hexadecimal  digits  in
        !          2608: the  result  of  each floating point operation to be set to zero.
        !          2609: For example, if the computer being  used  has  an  IBM  360  long
        !          2610: floating  point  word with 14 hexadecimal digits in the fraction,
        !          2611: then CHOP(8) results in simulation of  a  computer  with  only  6
        !          2612: hexadecimal  digits  in the fraction, i.e. a short floating point
        !          2613: word. On a computer such as the CDC 6600 with  16  octal  digits,
        !          2614: CHOP(8)  results in about the same accuracy because the remaining
        !          2615: 8 octal digits represent the same number of bits as 6 hexadecimal
        !          2616: digits.
        !          2617: 
        !          2618:      Some idea of the effect of CHOP on any particular system can
        !          2619: be obtained by executing the following statements.
        !          2620: 
        !          2621:       long,   t = 1/10
        !          2622:       long z, t = 1/10
        !          2623:       chop(8)
        !          2624:       long,   t = 1/10
        !          2625:       long z, t = 1/10
        !          2626: 
        !          2627: 
        !          2628:      The following Fortran subprograms illustrate more details of
        !          2629: FLOP  and CHOP. The first subprogram is a simplified example of a
        !          2630: system-dependent function used within MATLAB itself.  The  common
        !          2631: variable  FLP  is essentially the first component of the variable
        !          2632: FLOP.  The common variable CHP is initially zero, but it  is  set
        !          2633: to  p  by the statement  CHOP(p).  To shorten the DATA statement,
        !          2634: 
        !          2635: 
        !          2636: 
        !          2637: 
        !          2638: 
        !          2639: 
        !          2640: 
        !          2641: 
        !          2642: 
        !          2643: MATLAB, page 40
        !          2644: 
        !          2645: 
        !          2646: 
        !          2647: we assume there are only 6 hexadecimal digits.  We also assume an
        !          2648: extension  of  Fortran  that  allows .AND. to be used as a binary
        !          2649: operation between two real variables.
        !          2650: 
        !          2651:       REAL FUNCTION FLOP(X)
        !          2652:       REAL X
        !          2653:       INTEGER FLP,CHP
        !          2654:       COMMON FLP,CHP
        !          2655:       REAL MASK(5)
        !          2656:       DATA MASK/ZFFFFFFF0,ZFFFFFF00,ZFFFFF000,ZFFFF0000,ZFFF00000/
        !          2657:       FLP = FLP + 1
        !          2658:       IF (CHP .EQ. 0) FLOP = X
        !          2659:       IF (CHP .GE. 1 .AND. CHP .LE. 5) FLOP = X .AND. MASK(CHP)
        !          2660:       IF (CHP .GE. 6) FLOP = 0.0
        !          2661:       RETURN
        !          2662:       END
        !          2663: 
        !          2664: 
        !          2665:      The following subroutine illustrates a typical  use  of  the
        !          2666: previous  function  within MATLAB.  It is a simplified version of
        !          2667: the Basic Linear Algebra Subprogram that adds a  scalar  multiple
        !          2668: of  one  vector  to another.  We assume here that the vectors are
        !          2669: stored with a memory increment of one.
        !          2670: 
        !          2671:       SUBROUTINE SAXPY(N,TR,TI,XR,XI,YR,YI)
        !          2672:       REAL TR,TI,XR(N),XI(N),YR(N),YI(N),FLOP
        !          2673:       IF (N .LE. 0) RETURN
        !          2674:       IF (TR .EQ. 0.0 .AND. TI .EQ. 0.0) RETURN
        !          2675:       DO 10 I = 1, N
        !          2676:          YR(I) = FLOP(YR(I) + TR*XR(I) - TI*XI(I))
        !          2677:          YI(I) = YI(I) + TR*XI(I) + TI*XR(I)
        !          2678:          IF (YI(I) .NE. 0.0D0) YI(I) = FLOP(YI(I))
        !          2679:    10 CONTINUE
        !          2680:       RETURN
        !          2681:       END
        !          2682: 
        !          2683: 
        !          2684:      The  saxpy  operation  is  perhaps  the   most   fundamental
        !          2685: operation  within  LINPACK.  It is used in the computation of the
        !          2686: LU, the QR and the  SVD  factorizations,  and  in  several  other
        !          2687: places.   We  see  that  adding  a  multiple of one vector with n
        !          2688: components to another uses n flops if the vectors  are  real  and
        !          2689: between  n  and  2*n  flops if the vectors have nonzero imaginary
        !          2690: components.
        !          2691: 
        !          2692:      The permanent MATLAB variable EPS is reset by the  statement
        !          2693: CHOP(p).   Its new value is usually the smallest inverse power of
        !          2694: two that satisfies the Fortran logical test
        !          2695: 
        !          2696:             FLOP(1.0+EPS) .GT. 1.0
        !          2697: 
        !          2698: However, if EPS had been directly reset to a  larger  value,  the
        !          2699: old value is retained.
        !          2700: 
        !          2701: 
        !          2702: 
        !          2703: 
        !          2704: 
        !          2705: 
        !          2706: 
        !          2707: 
        !          2708: 
        !          2709: MATLAB, page 41
        !          2710: 
        !          2711: 
        !          2712: 
        !          2713: 
        !          2714: 
        !          2715: 13.  Communicating with other programs
        !          2716: 
        !          2717:      There  are  four  different  ways  MATLAB  can  be  used  in
        !          2718: conjunction with other programs:
        !          2719:       -- USER,
        !          2720:       -- EXEC,
        !          2721:       -- SAVE and LOAD,
        !          2722:       -- MATZ, CALL and RETURN .
        !          2723: 
        !          2724:      Let us illustrate each of  these  by  the  following  simple
        !          2725: example.
        !          2726: 
        !          2727:       n = 6
        !          2728:       for i = 1:n, for j = 1:n, a(i,j) = abs(i-j);
        !          2729:       A
        !          2730:       X = inv(A)
        !          2731: 
        !          2732: 
        !          2733:      The example  A  could be introduced into MATLAB  by  writing
        !          2734: the following Fortran subroutine.
        !          2735: 
        !          2736:          SUBROUTINE USER(A,M,N,S,T)
        !          2737:          DOUBLE PRECISION A(1),S,T
        !          2738:          N = IDINT(A(1))
        !          2739:          M = N
        !          2740:          DO 10 J = 1, N
        !          2741:          DO 10 I = 1, N
        !          2742:             K = I + (J-1)*M
        !          2743:             A(K) = IABS(I-J)
        !          2744:       10 CONTINUE
        !          2745:          RETURN
        !          2746:          END
        !          2747: 
        !          2748: This subroutine should be compiled  and  linked  into  MATLAB  in
        !          2749: place   of  the  original  version  of  USER.   Then  the  MATLAB
        !          2750: statements
        !          2751: 
        !          2752:       n = 6
        !          2753:       A = user(n)
        !          2754:       X = inv(A)
        !          2755: 
        !          2756: do the job.
        !          2757: 
        !          2758:      The example A could be generated by  storing  the  following
        !          2759: text in a file named, say, EXAMPLE .
        !          2760: 
        !          2761:       for i = 1:n, for j = 1:n, a(i,j) = abs(i-j);
        !          2762: 
        !          2763: Then the MATLAB statements
        !          2764: 
        !          2765:       n = 6
        !          2766: 
        !          2767: 
        !          2768: 
        !          2769: 
        !          2770: 
        !          2771: 
        !          2772: 
        !          2773: 
        !          2774: 
        !          2775: MATLAB, page 42
        !          2776: 
        !          2777: 
        !          2778: 
        !          2779:       exec('EXAMPLE',0)
        !          2780:       X = inv(A)
        !          2781: 
        !          2782: have the desired effect.  The 0 as the optional second  parameter
        !          2783: of exec indicates that the text in the file should not be printed
        !          2784: on the terminal.
        !          2785: 
        !          2786:      The matrices A and X could also be  stored  in  files.   Two
        !          2787: separate main programs would be involved.  The first is:
        !          2788: 
        !          2789:          PROGRAM MAINA
        !          2790:          DOUBLE PRECISION A(10,10)
        !          2791:          N = 6
        !          2792:          DO 10 J = 1, N
        !          2793:          DO 10 I = 1, N
        !          2794:             A(I,J) = IABS(I-J)
        !          2795:       10 CONTINUE
        !          2796:          OPEN(UNIT=1,FILE='A')
        !          2797:          WRITE(1,101) N,N
        !          2798:      101 FORMAT('A   ',2I4)
        !          2799:          DO 20 J = 1, N
        !          2800:             WRITE(1,102) (A(I,J),I=1,N)
        !          2801:       20 CONTINUE
        !          2802:      102 FORMAT(4Z18)
        !          2803:          END
        !          2804: 
        !          2805: The OPEN statement may take different forms on different systems.
        !          2806: It  attaches  Fortran  logical unit number 1 to the file named A.
        !          2807: The FORMAT  number  102  may  also  be  system  dependent.   This
        !          2808: particular one is appropriate for hexadecimal computers with an 8
        !          2809: byte double precision floating point  word.   Check,  or  modify,
        !          2810: MATLAB subroutine SAVLOD.
        !          2811: 
        !          2812:      After this program is executed, enter MATLAB  and  give  the
        !          2813: following statements:
        !          2814: 
        !          2815:       load('A')
        !          2816:       X = inv(A)
        !          2817:       save('X',X)
        !          2818: 
        !          2819: If all goes according to plan, this will read the matrix  A  from
        !          2820: the  file A, invert it, store the inverse in X and then write the
        !          2821: matrix X on the file X .  The following program can then access X
        !          2822: .
        !          2823: 
        !          2824:          PROGRAM MAINX
        !          2825:          DOUBLE PRECISION X(10,10)
        !          2826:          OPEN(UNIT=1,FILE='X')
        !          2827:          REWIND 1
        !          2828:          READ (1,101) ID,M,N
        !          2829:      101 FORMAT(A4,2I4)
        !          2830:          DO 10 J = 1, N
        !          2831:             READ(1,102) (X(I,J),I=1,M)
        !          2832: 
        !          2833: 
        !          2834: 
        !          2835: 
        !          2836: 
        !          2837: 
        !          2838: 
        !          2839: 
        !          2840: 
        !          2841: MATLAB, page 43
        !          2842: 
        !          2843: 
        !          2844: 
        !          2845:       10 CONTINUE
        !          2846:      102 FORMAT(4Z18)
        !          2847:          ...
        !          2848:          ...
        !          2849: 
        !          2850: 
        !          2851:      The most elaborate mechanism  involves  using  MATLAB  as  a
        !          2852: subroutine within another program.  Communication with the MATLAB
        !          2853: stack is accomplished using subroutine MATZ which is  distributed
        !          2854: with  MATLAB,  but  which  is  not  used  by  MATLAB itself.  The
        !          2855: preample of MATZ is:
        !          2856: 
        !          2857:       SUBROUTINE MATZ(A,LDA,M,N,IDA,JOB,IERR)
        !          2858:       INTEGER LDA,M,N,IDA(1),JOB,IERR
        !          2859:       DOUBLE PRECISION A(LDA,N)
        !          2860: C
        !          2861: C     ACCESS MATLAB VARIABLE STACK
        !          2862: C     A IS AN M BY N MATRIX, STORED IN AN ARRAY WITH
        !          2863: C         LEADING DIMENSION LDA.
        !          2864: C     IDA IS THE NAME OF A.
        !          2865: C         IF IDA IS AN INTEGER K LESS THAN 10, THEN THE NAME IS 'A'K
        !          2866: C         OTHERWISE, IDA(1:4) IS FOUR CHARACTERS, FORMAT 4A1.
        !          2867: C     JOB =  0  GET REAL A FROM MATLAB,
        !          2868: C         =  1  PUT REAL A INTO MATLAB,
        !          2869: C         = 10  GET IMAG PART OF A FROM MATLAB,
        !          2870: C         = 11  PUT IMAG PART OF A INTO MATLAB.
        !          2871: C     RETURN WITH NONZERO IERR AFTER MATLAB ERROR MESSAGE.
        !          2872: C
        !          2873: C     USES MATLAB ROUTINES STACKG, STACKP AND ERROR
        !          2874: 
        !          2875: 
        !          2876:      The preample of subroutine MATLAB is:
        !          2877: 
        !          2878: 
        !          2879:       SUBROUTINE MATLAB(INIT)
        !          2880: C     INIT = 0 FOR FIRST ENTRY, NONZERO FOR SUBSEQUENT ENTRIES
        !          2881: 
        !          2882: 
        !          2883:      To do our example, write the following program:
        !          2884: 
        !          2885:          DOUBLE PRECISION A(10,10),X(10,10)
        !          2886:          INTEGER IDA(4),IDX(4)
        !          2887:          DATA LDA/10/
        !          2888:          DATA IDA/'A',' ',' ',' '/, IDX/'X',' ',' ',' '/
        !          2889:          CALL MATLAB(0)
        !          2890:          N = 6
        !          2891:          DO 10 J = 1, N
        !          2892:          DO 10 I = 1, N
        !          2893:             A(I,J) = IABS(I-J)
        !          2894:       10 CONTINUE
        !          2895:          CALL MATZ(A,LDA,N,N,IDA,1,IERR)
        !          2896:          IF (IERR .NE. 0) GO TO ...
        !          2897:          CALL MATLAB(1)
        !          2898: 
        !          2899: 
        !          2900: 
        !          2901: 
        !          2902: 
        !          2903: 
        !          2904: 
        !          2905: 
        !          2906: 
        !          2907: MATLAB, page 44
        !          2908: 
        !          2909: 
        !          2910: 
        !          2911:          CALL MATZ(X,LDA,N,N,IDX,0,IERR)
        !          2912:          IF (IERR .NE. 0) GO TO ...
        !          2913:          ...
        !          2914:          ...
        !          2915: 
        !          2916: When this program is executed, the call to MATLAB(0) produces the
        !          2917: MATLAB greeting, then waits for input.  The command
        !          2918: 
        !          2919:          return
        !          2920: 
        !          2921: sends control back to our  example  program.   The  matrix  A  is
        !          2922: generated  by the program and sent to the stack by the first call
        !          2923: to MATZ.  The call to MATLAB(1) produces the MATLAB prompt.  Then
        !          2924: the statements
        !          2925: 
        !          2926:          X = inv(A)
        !          2927:          return
        !          2928: 
        !          2929: will invert our matrix, put the result on the stack and  go  back
        !          2930: to our program.  The second call to MATZ will retrieve X .
        !          2931: 
        !          2932:      By the way, this matrix  X  is interesting. Take a  look  at
        !          2933: round(2*(n-1)*X).
        !          2934: 
        !          2935: 
        !          2936: 
        !          2937: 
        !          2938: Acknowledgement.
        !          2939: 
        !          2940: 
        !          2941:      Most of the work on MATLAB  has  been  carried  out  at  the
        !          2942: University  of  New  Mexico,  where  it is being supported by the
        !          2943: National Science Foundation. Additional work has been done during
        !          2944: visits  to  Stanford  Linear Accelerator Center, Argonne National
        !          2945: Laboratory and Los Alamos Scientific  Laboratory,  where  support
        !          2946: has been provided by NSF and the Department of Energy.
        !          2947: 
        !          2948: 
        !          2949: References
        !          2950: 
        !          2951: [1]  J. J. Dongarra, J. R. Bunch, C. B. Moler and G. W.  Stewart,
        !          2952:      LINPACK  Users'  Guide,  Society  for Industrial and Applied
        !          2953:      Mathematics, Philadelphia, 1979.
        !          2954: 
        !          2955: [2]  B. T. Smith, J. M. Boyle, J. J. Dongarra, B. S.  Garbow,  Y.
        !          2956:      Ikebe, V. C. Klema, C. B. Moler, Matrix Eigensystem Routines
        !          2957:      -- EISPACK Guide, Lecture Notes in Computer Science,  volume
        !          2958:      6, second edition, Springer-Verlag, 1976.
        !          2959: 
        !          2960: [3]  B. S. Garbow, J. M. Boyle, J.  J.  Dongarra,  C.  B.  Moler,
        !          2961:      Matrix  Eigensystem  Routines  --  EISPACK  Guide Extension,
        !          2962:      Lecture Notes in  Computer  Science,  volume  51,  Springer-
        !          2963:      Verlag, 1977.
        !          2964: 
        !          2965: 
        !          2966: 
        !          2967: 
        !          2968: 
        !          2969: 
        !          2970: 
        !          2971: 
        !          2972: 
        !          2973: MATLAB, page 45
        !          2974: 
        !          2975: 
        !          2976: 
        !          2977: [4]  S. Cohen and  S.  Piper,  SPEAKEASY  III  Reference  Manual,
        !          2978:      Speakeasy Computing Corp., Chicago, Ill., 1979.
        !          2979: 
        !          2980: [5]  J. H. Wilkinson  and  C.  Reinsch,  Handbook  for  Automatic
        !          2981:      Computation,  volume  II,  Linear  Algebra, Springer-Verlag,
        !          2982:      1971.
        !          2983: 
        !          2984: [6]  Niklaus Wirth, Algorithms  +  Data  Structures  =  Programs,
        !          2985:      Prentice-Hall, 1976.
        !          2986: 
        !          2987: [7]  H. B. Keller and D. Sachs, "Calculations of the Conductivity
        !          2988:      of  a  Medium Containing Cylindrical Inclusions", J. Applied
        !          2989:      Physics 35, 537-538, 1964.
        !          2990: 
        !          2991: [8]  C. B. Moler and C. F. Van Loan,  Nineteen  Dubious  Ways  to
        !          2992:      Compute  the  Exponential  of a Matrix, SIAM Review 20, 801-
        !          2993:      836, 1979.
        !          2994: 
        !          2995: [9]  G. E. Forsythe, M. A. Malcolm  and  C.  B.  Moler,  Computer
        !          2996:      Methods for Mathematical Computations, Prentice-Hall, 1977.
        !          2997: 
        !          2998: [10] C. B. Moler and D. R. Morrison, "Replacing square  roots  by
        !          2999:      Pythagorean   sums",  University  of  New  Mexico,  Computer
        !          3000:      Science  Department,   technical   report,   submitted   for
        !          3001:      publication, 1980.
        !          3002: 
        !          3003: 
        !          3004: 
        !          3005: 
        !          3006: 
        !          3007: 
        !          3008: 
        !          3009: 
        !          3010: 
        !          3011: 
        !          3012: 
        !          3013: 
        !          3014: 
        !          3015: 
        !          3016: 
        !          3017: 
        !          3018: 
        !          3019: 
        !          3020: 
        !          3021: 
        !          3022: 
        !          3023: 
        !          3024: 
        !          3025: 
        !          3026: 
        !          3027: 
        !          3028: 
        !          3029: 
        !          3030: 
        !          3031: 
        !          3032: 
        !          3033: 
        !          3034: 
        !          3035: 
        !          3036: 
        !          3037: 
        !          3038: 
        !          3039: MATLAB, page 46
        !          3040: 
        !          3041: 
        !          3042: 
        !          3043: Appendix.  The HELP document
        !          3044: 
        !          3045: NEWS  MATLAB NEWS dated May, 1981.
        !          3046:       This describes recent or local changes.
        !          3047:       The new features added since the November,  1980,  printing
        !          3048:       of the Users' Guide include DIARY, EDIT, KRON, MACRO, PLOT,
        !          3049:       RAT, TRIL, TRIU and six element-by-element operations:
        !          3050:             .*   ./   .\   .*.   ./.   .\.
        !          3051:       Some additional  capabilities  have  been  added  to  EXIT,
        !          3052:       RANDOM, RCOND, SIZE and SVD.
        !          3053: 
        !          3054: INTRO Welcome to MATLAB.
        !          3055: 
        !          3056:       Here are a few sample statements:
        !          3057: 
        !          3058:       A = <1 2; 3 4>
        !          3059:       b = <5 6>'
        !          3060:       x = A\b
        !          3061:       <V,D> = eig(A),  norm(A-V*D/V)
        !          3062:       help \ , help eig
        !          3063:       exec('demo',7)
        !          3064: 
        !          3065:       For more information, see the MATLAB Users' Guide which  is
        !          3066:       contained in file ...  or may be obtained from ... .
        !          3067: 
        !          3068: HELP  HELP gives assistance.
        !          3069:       HELP HELP obviously prints this message.
        !          3070:       To see all the HELP messages, list the file ... .
        !          3071: 
        !          3072: <     < > Brackets used in forming vectors and matrices.
        !          3073:       <6.9  9.64  SQRT(-1)>  is  a  vector  with  three  elements
        !          3074:       separated  by  blanks.   <6.9,  9.64, sqrt(-1)> is the same
        !          3075:       thing.  <1+I 2-I 3>  and  <1 +I 2 -I 3>  are not the  same.
        !          3076:       The first has three elements, the second has five.
        !          3077:       <11 12 13; 21 22 23>  is a 2 by 3 matrix .   The  semicolon
        !          3078:       ends the first row.
        !          3079: 
        !          3080:       Vectors and matrices can be used inside < > brackets.
        !          3081:       <A B; C>  is allowed if the number of rows  of   A   equals
        !          3082:       the  number  of rows of  B  and the number of columns of  A
        !          3083:       plus the number of columns of   B   equals  the  number  of
        !          3084:       columns  of   C  .   This  rule  generalizes in a hopefully
        !          3085:       obvious way to allow fairly complicated constructions.
        !          3086: 
        !          3087:       A = < >  stores an empty matrix in  A , thereby removing it
        !          3088:       from the list of current variables.
        !          3089: 
        !          3090:       For the use of < and > on the left of  the  =  in  multiple
        !          3091:       assignment statements, see LU, EIG, SVD and so on.
        !          3092: 
        !          3093:       In WHILE and IF clauses, <>  means  less  than  or  greater
        !          3094:       than,  i.e.  not  equal, < means less than, > means greater
        !          3095:       than, <= means less than or equal, >= means greater than or
        !          3096: 
        !          3097: 
        !          3098: 
        !          3099: 
        !          3100: 
        !          3101: 
        !          3102: 
        !          3103: 
        !          3104: 
        !          3105: MATLAB, page 47
        !          3106: 
        !          3107: 
        !          3108: 
        !          3109:       equal.
        !          3110: 
        !          3111:       For the use of > and < to delineate macros, see MACRO.
        !          3112: 
        !          3113: >     See < .  Also see MACRO.
        !          3114: 
        !          3115: (     ( ) Used to indicate precedence in  arithmetic  expressions
        !          3116:       in  the  usual way.  Used to enclose arguments of functions
        !          3117:       in the usual way.  Used to enclose  subscripts  of  vectors
        !          3118:       and  matrices  in  a  manner somewhat more general than the
        !          3119:       usual way.  If  X   and   V  are  vectors,  then   X(V)  is
        !          3120:       <X(V(1)),  X(V(2)),  ...,  X(V(N))> .  The components of  V
        !          3121:       are rounded to nearest integers and used as subscripts.  An
        !          3122:       error  occurs  if  any  such  subscript  is  less than 1 or
        !          3123:       greater than the dimension of  X .  Some examples:
        !          3124:       X(3)  is the third element of  X .
        !          3125:       X(<1 2 3>)  is the first three elements of  X .  So is
        !          3126:       X(<SQRT(2), SQRT(3), 4*ATAN(1)>)  .
        !          3127:       If  X  has  N  components,  X(N:-1:1) reverses them.
        !          3128:       The same indirect subscripting is used in matrices.  If   V
        !          3129:       has   M  components and  W  has  N  components, then A(V,W)
        !          3130:       is the  M by N  matrix formed from the elements of A  whose
        !          3131:       subscripts are the elements of  V  and  W .  For example...
        !          3132:       A(<1,5>,:) = A(<5,1>,:)  interchanges rows 1 and 5 of  A .
        !          3133: 
        !          3134: )     See  ( .
        !          3135: 
        !          3136: =     Used in assignment statements and to mean equality in WHILE
        !          3137:       and IF clauses.
        !          3138: 
        !          3139: .     Decimal point.  314/100, 3.14  and   .314E1   are  all  the
        !          3140:       same.
        !          3141: 
        !          3142:       Element-by-element multiplicative operations  are  obtained
        !          3143:       using  .*  ,  ./  , or .\ .  For example, C = A ./ B is the
        !          3144:       matrix with elements  c(i,j) = a(i,j)/b(i,j) .
        !          3145: 
        !          3146:       Kronecker tensor products and quotients are  obtained  with
        !          3147:       .*. , ./.  and .\. .  See KRON.
        !          3148: 
        !          3149:       Two or  more  points  at  the  end  of  the  line  indicate
        !          3150:       continuation.    The   total  line  length  limit  is  1024
        !          3151:       characters.
        !          3152: 
        !          3153: ,     Used to separate matrix subscripts and function  arguments.
        !          3154:       Used  at  the  end  of  FOR, WHILE and IF clauses.  Used to
        !          3155:       separate statements  in  multi-statement  lines.   In  this
        !          3156:       situation,  it  may  be  replaced  by semicolon to suppress
        !          3157:       printing.
        !          3158: 
        !          3159: ;     Used inside brackets to end rows.
        !          3160:       Used after an expression or statement to suppress printing.
        !          3161:       See SEMI.
        !          3162: 
        !          3163: 
        !          3164: 
        !          3165: 
        !          3166: 
        !          3167: 
        !          3168: 
        !          3169: 
        !          3170: 
        !          3171: MATLAB, page 48
        !          3172: 
        !          3173: 
        !          3174: 
        !          3175: \     Backslash or matrix left division.   A\B   is  roughly  the
        !          3176:       same  as   INV(A)*B  , except it is computed in a different
        !          3177:       way.  If  A  is an N by N matrix and  B  is a column vector
        !          3178:       with  N  components, or a matrix with several such columns,
        !          3179:       then X = A\B  is the solution to  the  equation   A*X  =  B
        !          3180:       computed  by  Gaussian  elimination.   A warning message is
        !          3181:       printed if  A is badly scaled or nearly singular.
        !          3182:       A\EYE produces the inverse of  A .
        !          3183: 
        !          3184:       If  A  is an  M by N  matrix with  M < or > N  and  B  is a
        !          3185:       column vector with  M  components, or a matrix with several
        !          3186:       such columns, then  X = A\B  is the solution in  the  least
        !          3187:       squares  sense  to  the under- or overdetermined system  of
        !          3188:       equations A*X = B .  The  effective  rank,  K,  of   A   is
        !          3189:       determined  from  the  QR  decomposition  with pivoting.  A
        !          3190:       solution  X  is  computed  which  has  at  most  K  nonzero
        !          3191:       components  per column.  If  K < N this will usually not be
        !          3192:       the same solution as PINV(A)*B .
        !          3193:       A\EYE produces a generalized inverse of  A .
        !          3194: 
        !          3195:       If A and B have the  same  dimensions,  then  A  .\  B  has
        !          3196:       elements a(i,j)\b(i,j) .
        !          3197: 
        !          3198:       Also, see EDIT.
        !          3199: 
        !          3200: /     Slash or matrix right division.  B/A  is roughly  the  same
        !          3201:       as  B*INV(A) .  More precisely,  B/A = (A'\B')' .  See \ .
        !          3202: 
        !          3203:       IF A and B have the  same  dimensions,  then  A  ./  B  has
        !          3204:       elements a(i,j)/b(i,j) .
        !          3205: 
        !          3206:       Two or more slashes together on a line indicate  a  logical
        !          3207:       end of line.  Any following text is ignored.
        !          3208: 
        !          3209: '     Transpose.  X'  is the complex conjugate transpose of  X  .
        !          3210:       Quote.   'ANY  TEXT'   is a vector whose components are the
        !          3211:       MATLAB internal codes for the characters.  A  quote  within
        !          3212:       the text is indicated by two quotes.  See DISP and FILE .
        !          3213: 
        !          3214: +     Addition.  X + Y .  X and Y must have the same dimensions.
        !          3215: 
        !          3216: -     Subtraction.  X  -  Y  .   X  and  Y  must  have  the  same
        !          3217:       dimensions.
        !          3218: 
        !          3219: *     Matrix multiplication, X*Y .  Any scalar (1  by  1  matrix)
        !          3220:       may multiply anything.  Otherwise, the number of columns of
        !          3221:       X must equal the number of rows of Y .
        !          3222: 
        !          3223:       Element-by-element multiplication is obtained with X .* Y .
        !          3224: 
        !          3225:       The Kronecker tensor product is denoted by X .*. Y .
        !          3226: 
        !          3227:       Powers.  X**p  is  X  to the   p   power.   p   must  be  a
        !          3228: 
        !          3229: 
        !          3230: 
        !          3231: 
        !          3232: 
        !          3233: 
        !          3234: 
        !          3235: 
        !          3236: 
        !          3237: MATLAB, page 49
        !          3238: 
        !          3239: 
        !          3240: 
        !          3241:       scalar.  If  X  is a matrix, see  FUN .
        !          3242: 
        !          3243: :     Colon.  Used in subscripts,  FOR  iterations  and  possibly
        !          3244:       elsewhere.
        !          3245:       J:K  is the same as  <J, J+1, ..., K>
        !          3246:       J:K  is empty if  J > K .
        !          3247:       J:I:K  is the same as  <J, J+I, J+2I, ..., K>
        !          3248:       J:I:K  is empty if  I > 0 and J > K or if I < 0 and J < K .
        !          3249:       The colon notation can be used to pick out  selected  rows,
        !          3250:       columns and elements of vectors and matrices.
        !          3251:       A(:)  is all the  elements  of  A,  regarded  as  a  single
        !          3252:       column.
        !          3253:       A(:,J)  is the  J-th  column of A
        !          3254:       A(J:K)  is  A(J),A(J+1),...,A(K)
        !          3255:       A(:,J:K)  is  A(:,J),A(:,J+1),...,A(:,K) and so on.
        !          3256:       For the use of the colon in the FOR statement, See FOR .
        !          3257: 
        !          3258: ABS   ABS(X)  is the absolute value, or complex modulus,  of  the
        !          3259:       elements of X .
        !          3260: 
        !          3261: ANS   Variable created automatically  when  expressions  are  not
        !          3262:       assigned to anything else.
        !          3263: 
        !          3264: ATAN  ATAN(X)  is the arctangent of  X .  See FUN .
        !          3265: 
        !          3266: BASE  BASE(X,B) is a vector containing the base B  representation
        !          3267:       of   X  .   This is often used in conjunction with DISPLAY.
        !          3268:       DISPLAY(X,B)  is  the  same  as  DISPLAY(BASE(X,B)).    For
        !          3269:       example,    DISP(4*ATAN(1),16)   prints   the   hexadecimal
        !          3270:       representation of pi.
        !          3271: 
        !          3272: CHAR  CHAR(K)  requests  an  input  line  containing   a   single
        !          3273:       character  to  replace  MATLAB  character  number  K in the
        !          3274:       following table.  For example, CHAR(45) replaces backslash.
        !          3275:       CHAR(-K) replaces the alternate character number K.
        !          3276: 
        !          3277:                 K  character alternate name
        !          3278:               0 - 9   0 - 9    0 - 9   digits
        !          3279:              10 - 35  A - Z    a - z   letters
        !          3280:                36                      blank
        !          3281:                37       (        (     lparen
        !          3282:                38       )        )     rparen
        !          3283:                39       ;        ;     semi
        !          3284:                40       :        |     colon
        !          3285:                41       +        +     plus
        !          3286:                42       -        -     minus
        !          3287:                43       *        *     star
        !          3288:                44       /        /     slash
        !          3289:                45       \        $     backslash
        !          3290:                46       =        =     equal
        !          3291:                47       .        .     dot
        !          3292:                48       ,        ,     comma
        !          3293:                49       '        "     quote
        !          3294: 
        !          3295: 
        !          3296: 
        !          3297: 
        !          3298: 
        !          3299: 
        !          3300: 
        !          3301: 
        !          3302: 
        !          3303: MATLAB, page 50
        !          3304: 
        !          3305: 
        !          3306: 
        !          3307:                50       <        [     less
        !          3308:                51       >        ]     great
        !          3309: 
        !          3310: CHOL  Cholesky factorization.  CHOL(X)  uses  only  the  diagonal
        !          3311:       and upper triangle of  X .  The lower triangular is assumed
        !          3312:       to be the (complex conjugate) transpose of the  upper.   If
        !          3313:       X   is  positive  definite,  then  R = CHOL(X)  produces an
        !          3314:       upper triangular  R  so that  R'*R = X .   If   X   is  not
        !          3315:       positive definite, an error message is printed.
        !          3316: 
        !          3317: CHOP  Truncate arithmetic.  CHOP(P) causes P places to be chopped
        !          3318:       off   after   each   arithmetic   operation  in  subsequent
        !          3319:       computations.  This means  P  hexadecimal  digits  on  some
        !          3320:       computers  and  P octal digits on others.  CHOP(0) restores
        !          3321:       full precision.
        !          3322: 
        !          3323: CLEAR Erases all variables, except EPS, FLOP, EYE and RAND.
        !          3324:       X = <>  erases only variable  X .  So does CLEAR X .
        !          3325: 
        !          3326: COND  Condition number in 2-norm.  COND(X) is the  ratio  of  the
        !          3327:       largest singular value of  X  to the smallest.
        !          3328: 
        !          3329: CONJG CONJG(X)  is the complex conjugate of  X .
        !          3330: 
        !          3331: COS   COS(X)  is the cosine of  X .  See FUN .
        !          3332: 
        !          3333: DET   DET(X)  is the determinant of the square matrix  X .
        !          3334: 
        !          3335: DIAG  If  V  is  a  row  or  column  vector  with  N  components,
        !          3336:       DIAG(V,K)   is a square matrix of order  N+ABS(K)  with the
        !          3337:       elements of  V  on the K-th diagonal.  K = 0  is  the  main
        !          3338:       diagonal,  K  >  0  is above the main diagonal and K < 0 is
        !          3339:       below the main diagonal.  DIAG(V)  simply puts  V   on  the
        !          3340:       main diagonal.
        !          3341:       eg. DIAG(-M:M) + DIAG(ONES(2*M,1),1) + DIAG(ONES(2*M,1),-1)
        !          3342:       produces a tridiagonal matrix of order 2*M+1 .
        !          3343:       IF  X  is a matrix,  DIAG(X,K)  is a column  vector  formed
        !          3344:       from the elements of the K-th diagonal of  X .
        !          3345:       DIAG(X)  is the main diagonal of  X .
        !          3346:       DIAG(DIAG(X))  is a diagonal matrix .
        !          3347: 
        !          3348: DIARY DIARY('file') causes a  copy  of  all  subsequent  terminal
        !          3349:       input and most of the resulting output to be written on the
        !          3350:       file. DIARY(0) turns it off.  See FILE.
        !          3351: 
        !          3352: DISP  DISPLAY(X) prints X  in  a  compact  format.   If  all  the
        !          3353:       elements  of  X  are  integers  between 0 and 51, then X is
        !          3354:       interpreted  as  MATLAB  text  and   printed   accordingly.
        !          3355:       Otherwise,  +  ,  -   and  blank  are printed for positive,
        !          3356:       negative and zero elements.  Imaginary parts are ignored.
        !          3357:       DISP(X,B) is the same as DISP(BASE(X,B)).
        !          3358: 
        !          3359: EDIT  There  are  no   editing   features   available   on   most
        !          3360: 
        !          3361: 
        !          3362: 
        !          3363: 
        !          3364: 
        !          3365: 
        !          3366: 
        !          3367: 
        !          3368: 
        !          3369: MATLAB, page 51
        !          3370: 
        !          3371: 
        !          3372: 
        !          3373:       installations and EDIT is not a command.  However, on a few
        !          3374:       systems a command line consisting of a single  backslash  \
        !          3375:       will  cause  the local file editor to be called with a copy
        !          3376:       of the  previous  input  line.   When  the  editor  returns
        !          3377:       control to MATLAB, it will execute the line again.
        !          3378: 
        !          3379: EIG   Eigenvalues and eigenvectors.
        !          3380:       EIG(X) is a vector containing the eigenvalues of  a  square
        !          3381:       matrix  X .
        !          3382:       <V,D>  =  EIG(X)   produces  a  diagonal  matrix    D    of
        !          3383:       eigenvalues  and  a  full  matrix  V  whose columns are the
        !          3384:       corresponding eigenvectors so that  X*V = V*D .
        !          3385: 
        !          3386: ELSE  Used with IF .
        !          3387: 
        !          3388: END   Terminates the scope  of  FOR,  WHILE  and  IF  statements.
        !          3389:       Without  END's,  FOR  and WHILE repeat all statements up to
        !          3390:       the end of the line.  Each END is paired with  the  closest
        !          3391:       previous  unpaired FOR or WHILE and serves to terminate its
        !          3392:       scope.  The line
        !          3393:       FOR I=1:N, FOR J=1:N, A(I,J)=1/(I+J-1); A
        !          3394:       would cause A to be printed  N**2  times, once for each new
        !          3395:       element.  On the other hand, the line
        !          3396:       FOR I=1:N, FOR J=1:N, A(I,J)=1/(I+J-1); END, END, A
        !          3397:       will lead to only the final printing of  A .
        !          3398:       Similar considerations apply to WHILE.
        !          3399:       EXIT terminates execution of loops or of MATLAB itself.
        !          3400: 
        !          3401: EPS   Floating point relative  accuracy.   A  permanent  variable
        !          3402:       whose  value is initially the distance from 1.0 to the next
        !          3403:       largest floating point number.  The  value  is  changed  by
        !          3404:       CHOP,  and  other values may be assigned.  EPS is used as a
        !          3405:       default tolerance by PINV and RANK.
        !          3406: 
        !          3407: EXEC  EXEC('file',k) obtains  subsequent  MATLAB  input  from  an
        !          3408:       external  file.  The printing of input is controlled by the
        !          3409:       optional parameter k .
        !          3410:       If k = 1 , the input is echoed.
        !          3411:       If k = 2 , the MATLAB prompt <> is printed.
        !          3412:       If k = 4 , MATLAB pauses before each prompt and waits for a
        !          3413:       null line to continue.
        !          3414:       If k = 0 , there is no echo, prompt or pause.  This is  the
        !          3415:       default if the exec command is followed by a semicolon.
        !          3416:       If k = 7 , there will be echos, prompts and pauses. This is
        !          3417:       useful for demonstrations on video terminals.
        !          3418:       If k = 3 , there will be echos and prompts, but no  pauses.
        !          3419:       This is the the default if the exec command is not followed
        !          3420:       by a semicolon.
        !          3421:       EXEC(0) causes subsequent input to  be  obtained  from  the
        !          3422:       terminal. An end-of-file has the same effect.
        !          3423:       EXEC's may be nested, i.e. the text in the file may contain
        !          3424:       EXEC of another file.  EXEC's may also be driven by FOR and
        !          3425:       WHILE loops.
        !          3426: 
        !          3427: 
        !          3428: 
        !          3429: 
        !          3430: 
        !          3431: 
        !          3432: 
        !          3433: 
        !          3434: 
        !          3435: MATLAB, page 52
        !          3436: 
        !          3437: 
        !          3438: 
        !          3439: EXIT  Causes termination of a FOR or WHILE loop.
        !          3440:       If not in a loop, terminates execution of MATLAB.
        !          3441: 
        !          3442: EXP   EXP(X)  is the exponential of  X ,  e  to the X .  See  FUN
        !          3443:       .
        !          3444: 
        !          3445: EYE   Identity matrix.  EYE(N) is the N  by  N  identity  matrix.
        !          3446:       EYE(M,N)   is an M by N matrix with 1's on the diagonal and
        !          3447:       zeros elsewhere.  EYE(A)  is the same size  as   A  .   EYE
        !          3448:       with  no  arguments is an identity matrix of whatever order
        !          3449:       is appropriate in the context.   For  example,  A  +  3*EYE
        !          3450:       adds  3  to each diagonal element of  A .
        !          3451: 
        !          3452: FILE  The EXEC, SAVE, LOAD,  PRINT  and  DIARY  functions  access
        !          3453:       files.   The  'file'  parameter  takes  different forms for
        !          3454:       different operating systems.  On most systems,  'file'  may
        !          3455:       be a string of up to 32 characters in quotes.  For example,
        !          3456:       SAVE('A') or EXEC('matlab/demo.exec') .  The string will be
        !          3457:       used as the name of a file in the local operating system.
        !          3458:       On all systems, 'file' may be a positive integer   k   less
        !          3459:       than  10  which  will  be  used  as  a FORTRAN logical unit
        !          3460:       number. Some systems then automatically access a file  with
        !          3461:       a  name  like  FORT.k  or FORk.DAT. Other systems require a
        !          3462:       file with a name like FT0kF001 to be assigned  to  unit   k
        !          3463:       before  MATLAB  is  executed. Check your local installation
        !          3464:       for details.
        !          3465: 
        !          3466: FLOPS Count of floating point operations.
        !          3467:       FLOPS  is  a  permanently  defined  row  vector  with   two
        !          3468:       elements.    FLOPS(1)  is  the  number  of  floating  point
        !          3469:       operations counted during the previous statement.  FLOPS(2)
        !          3470:       is  a  cumulative total.  FLOPS can be used in the same way
        !          3471:       as any other vector.  FLOPS(2) = 0  resets  the  cumulative
        !          3472:       total.   In  addition,  FLOPS(1) will be printed whenever a
        !          3473:       statement is terminated by an extra comma.  For example,
        !          3474:       X = INV(A);,
        !          3475:       or
        !          3476:       COND(A),   (as the last statement on the line).
        !          3477:       HELP FLPS gives more details.
        !          3478: 
        !          3479: FLPS  More detail on FLOPS.
        !          3480:       It is not feasible to count absolutely all  floating  point
        !          3481:       operations,  but  most  of  the important ones are counted.
        !          3482:       Each multiply and add in a real vector operation such as  a
        !          3483:       dot  product  or  a 'saxpy' counts one flop.  Each multiply
        !          3484:       and add in a complex vector  operation  counts  two  flops.
        !          3485:       Other additions, subtractions and multiplications count one
        !          3486:       flop each if the result is real and two flops if it is not.
        !          3487:       Real  divisions  count one and complex divisions count two.
        !          3488:       Elementary functions count one if real and two if  complex.
        !          3489:       Some examples.  If A and B are real N by N matrices, then
        !          3490:       A + B  counts N**2 flops,
        !          3491:       A*B    counts N**3 flops,
        !          3492: 
        !          3493: 
        !          3494: 
        !          3495: 
        !          3496: 
        !          3497: 
        !          3498: 
        !          3499: 
        !          3500: 
        !          3501: MATLAB, page 53
        !          3502: 
        !          3503: 
        !          3504: 
        !          3505:       A**100 counts 99*N**3 flops,
        !          3506:       LU(A)  counts roughly (1/3)*N**3 flops.
        !          3507: 
        !          3508: FOR   Repeat statements a specific number of times.
        !          3509:       FOR variable = expr, statement, ..., statement, END
        !          3510:       The END at the end of a line may  be  omitted.   The  comma
        !          3511:       before  the  END  may  also be omitted.  The columns of the
        !          3512:       expression are stored one at a time  in  the  variable  and
        !          3513:       then the following statements, up to the END, are executed.
        !          3514:       The expression is often of the form X:Y, in which case  its
        !          3515:       columns  are  simply  scalars.  Some examples (assume N has
        !          3516:       already been assigned a value).
        !          3517:       FOR I = 1:N, FOR J = 1:N, A(I,J) = 1/(I+J-1);
        !          3518:       FOR J = 2:N-1, A(J,J) = J; END; A
        !          3519:       FOR S = 1.0: -0.1: 0.0, ...  steps S with increments of -0.1 .
        !          3520:       FOR E = EYE(N), ...   sets  E  to the unit N-vectors.
        !          3521:       FOR V = A, ...   has the same effect as
        !          3522:       FOR J = 1:N, V = A(:,J); ...  except J is also set here.
        !          3523: 
        !          3524: FUN   For matrix arguments  X , the  functions  SIN,  COS,  ATAN,
        !          3525:       SQRT,  LOG,  EXP and X**p are computed using eigenvalues  D
        !          3526:       and eigenvectors  V .  If  <V,D> =  EIG(X)   then   f(X)  =
        !          3527:       V*f(D)/V  .   This method may give inaccurate results if  V
        !          3528:       is badly conditioned.  Some idea of  the  accuracy  can  be
        !          3529:       obtained by comparing  X**1  with  X .
        !          3530:       For vector arguments,  the  function  is  applied  to  each
        !          3531:       component.
        !          3532: 
        !          3533: HESS  Hessenberg form.  The Hessenberg form of a matrix  is  zero
        !          3534:       below the first subdiagonal.  If the matrix is symmetric or
        !          3535:       Hermitian,  the  form  is  tridiagonal.   <P,H>  =  HESS(A)
        !          3536:       produces  a  unitary  matrix P and a Hessenberg matrix H so
        !          3537:       that A = P*H*P'.  By itself, HESS(A) returns H.
        !          3538: 
        !          3539: HILB  Inverse Hilbert matrix.  HILB(N)  is the inverse of  the  N
        !          3540:       by  N   matrix  with elements  1/(i+j-1), which is a famous
        !          3541:       example of a badly conditioned matrix.  The result is exact
        !          3542:       for  N  less than about 15, depending upon the computer.
        !          3543: 
        !          3544: IF    Conditionally execute statements.  Simple form...
        !          3545:       IF expression rop expression, statements
        !          3546:       where rop is =, <, >, <=, >=, or  <>  (not  equal)  .   The
        !          3547:       statements  are  executed  once if the indicated comparison
        !          3548:       between the real parts of the first components of  the  two
        !          3549:       expressions  is true, otherwise the statements are skipped.
        !          3550:       Example.
        !          3551:       IF ABS(I-J) = 1, A(I,J) = -1;
        !          3552:       More complicated forms use END in the same way it  is  used
        !          3553:       with FOR and WHILE and use ELSE as an abbreviation for END,
        !          3554:       IF expression not rop expression .  Example
        !          3555:       FOR I = 1:N, FOR J = 1:N, ...
        !          3556:          IF I = J, A(I,J) = 2; ELSE IF ABS(I-J) = 1, A(I,J) = -1; ...
        !          3557:          ELSE A(I,J) = 0;
        !          3558: 
        !          3559: 
        !          3560: 
        !          3561: 
        !          3562: 
        !          3563: 
        !          3564: 
        !          3565: 
        !          3566: 
        !          3567: MATLAB, page 54
        !          3568: 
        !          3569: 
        !          3570: 
        !          3571:       An easier way to accomplish the same thing is
        !          3572:       A = 2*EYE(N);
        !          3573:       FOR I = 1:N-1, A(I,I+1) = -1; A(I+1,I) = -1;
        !          3574: 
        !          3575: IMAG  IMAG(X)  is the imaginary part of  X .
        !          3576: 
        !          3577: INV   INV(X)  is the inverse of the square matrix  X .  A warning
        !          3578:       message  is  printed  if   X   is  badly  scaled  or nearly
        !          3579:       singular.
        !          3580: 
        !          3581: KRON  KRON(X,Y) is the Kronecker tensor product of X and Y  .  It
        !          3582:       is  also  denoted by X .*. Y . The result is a large matrix
        !          3583:       formed by taking all possible products between the elements
        !          3584:       of  X  and  those  of Y . For example, if X is 2 by 3, then
        !          3585:       X .*. Y is
        !          3586: 
        !          3587:             < x(1,1)*Y  x(1,2)*Y  x(1,3)*Y
        !          3588:               x(2,1)*Y  x(2,2)*Y  x(2,3)*Y >
        !          3589: 
        !          3590:       The five-point discrete Laplacian for an n-by-n grid can be
        !          3591:       generated by
        !          3592: 
        !          3593:             T = diag(ones(n-1,1),1);  T = T + T';  I = EYE(T);
        !          3594:             A = T.*.I + I.*.T - 4*EYE;
        !          3595: 
        !          3596:       Just  in  case  they  might  be  useful,  MATLAB   includes
        !          3597:       constructions called Kronecker tensor quotients, denoted by
        !          3598:       X ./. Y and X .\. Y .  They are obtained by  replacing  the
        !          3599:       elementwise multiplications in X .*. Y with divisions.
        !          3600: 
        !          3601: LINES An internal count is kept of the number of lines of  output
        !          3602:       since  the  last  input.   Whenever this count approaches a
        !          3603:       limit, the  user  is  asked  whether  or  not  to  suppress
        !          3604:       printing  until the next input.  Initially the limit is 25.
        !          3605:       LINES(N) resets the limit to N .
        !          3606: 
        !          3607: LOAD  LOAD('file') retrieves all the variables from  the  file  .
        !          3608:       See  FILE  and  SAVE for more details.  To prepare your own
        !          3609:       file for LOADing, change the READs to WRITEs  in  the  code
        !          3610:       given under SAVE.
        !          3611: 
        !          3612: LOG   LOG(X)  is the  natural  logarithm  of   X  .   See  FUN  .
        !          3613:       Complex results are produced if  X  is not positive, or has
        !          3614:       nonpositive eigenvalues.
        !          3615: 
        !          3616: LONG  Determine output format.   All  computations  are  done  in
        !          3617:       complex arithmetic and double precision if it is available.
        !          3618:       SHORT and  LONG  merely  switch  between  different  output
        !          3619:       formats.
        !          3620:       SHORT    Scaled fixed point format with about 5 digits.
        !          3621:       LONG     Scaled fixed point format with about 15 digits.
        !          3622:       SHORT E  Floating point format with about 5 digits.
        !          3623:       LONG E   Floating point format with about 15 digits.
        !          3624: 
        !          3625: 
        !          3626: 
        !          3627: 
        !          3628: 
        !          3629: 
        !          3630: 
        !          3631: 
        !          3632: 
        !          3633: MATLAB, page 55
        !          3634: 
        !          3635: 
        !          3636: 
        !          3637:       LONG Z   System dependent format, often hexadecimal.
        !          3638: 
        !          3639: LU    Factors from Gaussian elimination.  <L,U> = LU(X)  stores a
        !          3640:       upper triangular matrix in  U  and a 'psychologically lower
        !          3641:       triangular matrix', i.e. a product of lower triangular  and
        !          3642:       permutation matrices, in L , so that  X = L*U .  By itself,
        !          3643:       LU(X) returns the output from CGEFA .
        !          3644: 
        !          3645: MACRO The macro facility involves text and inward pointing  angle
        !          3646:       brackets.  If  STRING  is  the  source  text for any MATLAB
        !          3647:       expression or statement, then
        !          3648:             t = 'STRING';
        !          3649:       encodes the text as a vector of integers  and  stores  that
        !          3650:       vector in  t .  DISP(t) will print the text and
        !          3651:             >t<
        !          3652:       causes the text to be interpreted, either as a statement or
        !          3653:       as a factor in an expression.  For example
        !          3654:             t = '1/(i+j-1)';
        !          3655:             disp(t)
        !          3656:             for i = 1:n, for j = 1:n, a(i,j) = >t<;
        !          3657:       generates the Hilbert matrix of order n.
        !          3658:       Another example showing indexed text,
        !          3659:             S = <'x = 3            '
        !          3660:                  'y = 4            '
        !          3661:                  'z = sqrt(x*x+y*y)'>
        !          3662:             for k = 1:3, >S(k,:)<
        !          3663:       It is necessary that the strings making up  the  "rows"  of
        !          3664:       the "matrix"  S  have the same lengths.
        !          3665: 
        !          3666: MAGIC Magic square.  MAGIC(N) is an N  by  N  matrix  constructed
        !          3667:       from  the integers 1 through N**2 with equal row and column
        !          3668:       sums.
        !          3669: 
        !          3670: NORM  For matrices..
        !          3671:       NORM(X)  is the largest singular value of  X .
        !          3672:       NORM(X,1)  is the 1-norm of  X .
        !          3673:       NORM(X,2)  is the same as NORM(X) .
        !          3674:       NORM(X,'INF')  is the infinity norm of  X .
        !          3675:       NORM(X,'FRO')  is the F-norm, i.e.  SQRT(SUM(DIAG(X'*X))) .
        !          3676:       For vectors..
        !          3677:       NORM(V,P) = (SUM(V(I)**P))**(1/P) .
        !          3678:       NORM(V) = NORM(V,2) .
        !          3679:       NORM(V,'INF') = MAX(ABS(V(I))) .
        !          3680: 
        !          3681: ONES  All ones.  ONES(N)  is an N by N matrix of ones.  ONES(M,N)
        !          3682:       is an M by N matrix of ones .  ONES(A)  is the same size as
        !          3683:       A  and all ones .
        !          3684: 
        !          3685: ORTH  Orthogonalization.   Q  =  ORTH(X)   is   a   matrix   with
        !          3686:       orthonormal  columns,  i.e. Q'*Q = EYE, which span the same
        !          3687:       space as the columns of  X .
        !          3688: 
        !          3689: PINV  Pseudoinverse.  X = PINV(A) produces a matrix   X   of  the
        !          3690: 
        !          3691: 
        !          3692: 
        !          3693: 
        !          3694: 
        !          3695: 
        !          3696: 
        !          3697: 
        !          3698: 
        !          3699: MATLAB, page 56
        !          3700: 
        !          3701: 
        !          3702: 
        !          3703:       same  dimensions as  A' so that  A*X*A = A , X*A*X = X  and
        !          3704:       AX  and  XA  are Hermitian .  The computation is  based  on
        !          3705:       SVD(A)  and  any  singular values less than a tolerance are
        !          3706:       treated   as    zero.     The    default    tolerance    is
        !          3707:       NORM(SIZE(A),'inf')*NORM(A)*EPS.   This  tolerance  may  be
        !          3708:       overridden with X = PINV(A,tol).  See RANK.
        !          3709: 
        !          3710: PLOT  PLOT(X,Y) produces a plot of  the  elements  of  Y  against
        !          3711:       those  of X . PLOT(Y) is the same as PLOT(1:n,Y) where n is
        !          3712:       the  number  of   elements   in   Y   .    PLOT(X,Y,P)   or
        !          3713:       PLOT(X,Y,p1,...,pk)  passes the optional parameter vector P
        !          3714:       or scalars p1 through pk to the plot routine.  The  default
        !          3715:       plot  routine  is a crude printer-plot. It is hoped that an
        !          3716:       interface to local graphics equipment can be provided.
        !          3717:       An interesting example is
        !          3718:             t = 0:50;
        !          3719:             PLOT( t.*cos(t), t.*sin(t) )
        !          3720: 
        !          3721: POLY  Characteristic polynomial.
        !          3722:       If  A  is an N by N matrix, POLY(A) is a column vector with
        !          3723:       N+1   elements   which   are   the   coefficients   of  the
        !          3724:       characteristic polynomial,  DET(lambda*EYE - A) .
        !          3725:       If V is a vector, POLY(V) is a vector  whose  elements  are
        !          3726:       the  coefficients  of  the  polynomial  whose roots are the
        !          3727:       elements of V .  For vectors, ROOTS and  POLY  are  inverse
        !          3728:       functions  of  each  other,  up  to  ordering, scaling, and
        !          3729:       roundoff error.
        !          3730:       ROOTS(POLY(1:20)) generates Wilkinson's famous example.
        !          3731: 
        !          3732: PRINT PRINT('file',X) prints X on  the  file  using  the  current
        !          3733:       format determined by SHORT, LONG Z, etc.  See FILE.
        !          3734: 
        !          3735: PROD  PROD(X)  is the product of all the elements of  X .
        !          3736: 
        !          3737: QR    Orthogonal-triangular decomposition.
        !          3738:       <Q,R> = QR(X)  produces an upper triangular  matrix   R  of
        !          3739:       the  same dimension as  X  and a unitary matrix  Q  so that
        !          3740:       X = Q*R .
        !          3741:       <Q,R,E> = QR(X)  produces a  permutation  matrix   E  ,  an
        !          3742:       upper  triangular  R  with decreasing diagonal elements and
        !          3743:       a unitary  Q  so that  X*E = Q*R .
        !          3744:       By itself, QR(X) returns the output of CQRDC .  TRIU(QR(X))
        !          3745:       is R .
        !          3746: 
        !          3747: RAND  Random numbers and matrices.  RAND(N)  is an N by N  matrix
        !          3748:       with  random  entries.  RAND(M,N)  is an M by N matrix with
        !          3749:       random entries.  RAND(A)  is the same size as   A  .   RAND
        !          3750:       with no arguments is a scalar whose value changes each time
        !          3751:       it is referenced.
        !          3752:       Ordinarily,  random numbers are  uniformly  distributed  in
        !          3753:       the  interval  (0.0,1.0)  .   RAND('NORMAL')  switches to a
        !          3754:       normal distribution  with  mean  0.0  and  variance  1.0  .
        !          3755:       RAND('UNIFORM')  switches back to the uniform distribution.
        !          3756: 
        !          3757: 
        !          3758: 
        !          3759: 
        !          3760: 
        !          3761: 
        !          3762: 
        !          3763: 
        !          3764: 
        !          3765: MATLAB, page 57
        !          3766: 
        !          3767: 
        !          3768: 
        !          3769:       RAND('SEED') returns the current value of the seed for  the
        !          3770:       generator.    RAND('SEED',n)   sets   the   seed   to  n  .
        !          3771:       RAND('SEED',0) resets the seed to 0, its value when  MATLAB
        !          3772:       is first entered.
        !          3773: 
        !          3774: RANK  Rank.  K = RANK(X) is the number of singular values  of   X
        !          3775:       that are larger than NORM(SIZE(X),'inf')*NORM(X)*EPS.
        !          3776:       K = RANK(X,tol) is the number of singular values of  X that
        !          3777:       are larger than tol .
        !          3778: 
        !          3779: RCOND RCOND(X)   is  an  estimate  for  the  reciprocal  of   the
        !          3780:       condition  of   X   in  the  1-norm obtained by the LINPACK
        !          3781:       condition estimator.  If  X  is well conditioned,  RCOND(X)
        !          3782:       is  near  1.0  .   If  X  is badly conditioned, RCOND(X) is
        !          3783:       near 0.0 .
        !          3784:       <R, Z> = RCOND(A) sets  R  to RCOND(A) and also produces  a
        !          3785:       vector  Z so that
        !          3786:                  NORM(A*Z,1) = R*NORM(A,1)*NORM(Z,1)
        !          3787:       So, if RCOND(A) is small, then  Z  is an  approximate  null
        !          3788:       vector.
        !          3789: 
        !          3790: RAT   An experimental  function  which  attempts  to  remove  the
        !          3791:       roundoff   error  from  results  that  should  be  "simple"
        !          3792:       rational numbers.
        !          3793:       RAT(X) approximates each  element  of   X  by  a  continued
        !          3794:       fraction of the form
        !          3795: 
        !          3796:                 a/b = d1 + 1/(d2 + 1/(d3 + ... + 1/dk))
        !          3797: 
        !          3798:       with k <= len, integer di and abs(di) <= max .  The default
        !          3799:       values of the parameters are len = 5 and max = 100.
        !          3800:       RAT(len,max) changes the default values.  Increasing either
        !          3801:       len or max increases the number of possible fractions.
        !          3802:       <A,B> = RAT(X) produces integer matrices A and B so that
        !          3803: 
        !          3804:                 A ./ B  =  RAT(X)
        !          3805: 
        !          3806:       Some examples:
        !          3807: 
        !          3808:             long
        !          3809:             T = hilb(6), X = inv(T)
        !          3810:             <A,B> = rat(X)
        !          3811:             H = A ./ B, S = inv(H)
        !          3812: 
        !          3813:             short e
        !          3814:             d = 1:8,  e = ones(d),  A = abs(d'*e - e'*d)
        !          3815:             X = inv(A)
        !          3816:             rat(X)
        !          3817:             display(ans)
        !          3818: 
        !          3819: 
        !          3820: REAL  REAL(X)  is the real part of  X .
        !          3821: 
        !          3822: 
        !          3823: 
        !          3824: 
        !          3825: 
        !          3826: 
        !          3827: 
        !          3828: 
        !          3829: 
        !          3830: 
        !          3831: MATLAB, page 58
        !          3832: 
        !          3833: 
        !          3834: 
        !          3835: RETURN  From the terminal, causes return to the operating  system
        !          3836:       or  other  program  which  invoked  MATLAB.  From inside an
        !          3837:       EXEC, causes  return  to  the  invoking  EXEC,  or  to  the
        !          3838:       terminal.
        !          3839: 
        !          3840: RREF  RREF(A) is the reduced row echelon form of the  rectangular
        !          3841:       matrix.  RREF(A,B) is the same as RREF(<A,B>) .
        !          3842: 
        !          3843: ROOTS Find polynomial roots.  ROOTS(C)  computes the roots of the
        !          3844:       polynomial  whose  coefficients  are  the  elements  of the
        !          3845:       vector  C .  If  C  has  N+1  components, the polynomial is
        !          3846:       C(1)*X**N + ... + C(N)*X + C(N+1) .  See POLY.
        !          3847: 
        !          3848: ROUND ROUND(X)  rounds  the  elements  of   X   to  the   nearest
        !          3849:       integers.
        !          3850: 
        !          3851: SAVE  SAVE('file') stores all the current variables in a file.
        !          3852:       SAVE('file',X) saves only X .  See FILE .
        !          3853:       The variables may be retrieved later by LOAD('file') or  by
        !          3854:       your  own program using the following code for each matrix.
        !          3855:       The lines involving XIMAG may be eliminated  if  everything
        !          3856:       is known to be real.
        !          3857: 
        !          3858:             attach lunit to 'file'
        !          3859:             REAL or DOUBLE PRECISION XREAL(MMAX,NMAX)
        !          3860:             REAL or DOUBLE PRECISION XIMAG(MMAX,NMAX)
        !          3861:             READ(lunit,101) ID,M,N,IMG
        !          3862:             DO 10 J = 1, N
        !          3863:                READ(lunit,102) (XREAL(I,J), I=1,M)
        !          3864:                IF (IMG .NE. 0) READ(lunit,102) (XIMAG(I,J),I=1,M)
        !          3865:          10 CONTINUE
        !          3866: 
        !          3867:       The formats used are system dependent.  The  following  are
        !          3868:       typical.     See    SUBROUTINE   SAVLOD   in   your   local
        !          3869:       implementation of MATLAB.
        !          3870: 
        !          3871:         101 FORMAT(4A1,3I4)
        !          3872:         102 FORMAT(4Z18)
        !          3873:         102 FORMAT(4O20)
        !          3874:         102 FORMAT(4D25.18)
        !          3875: 
        !          3876: SCHUR Schur decomposition.  <U,T> = SCHUR(X)  produces  an  upper
        !          3877:       triangular  matrix   T , with the eigenvalues of  X  on the
        !          3878:       diagonal, and a unitary matrix  U so that  X =  U*T*U'  and
        !          3879:       U'*U = EYE .  By itself, SCHUR(X) returns  T .
        !          3880: 
        !          3881: SHORT See LONG .
        !          3882: 
        !          3883: SEMI  Semicolons at the end of  lines  will  cause,  rather  than
        !          3884:       suppress,  printing.   A  second  SEMI restores the initial
        !          3885:       interpretation.
        !          3886: 
        !          3887: SIN   SIN(X)  is the sine of  X .  See FUN .
        !          3888: 
        !          3889: 
        !          3890: 
        !          3891: 
        !          3892: 
        !          3893: 
        !          3894: 
        !          3895: 
        !          3896: 
        !          3897: MATLAB, page 59
        !          3898: 
        !          3899: 
        !          3900: 
        !          3901: SIZE  If X is an M by N matrix, then SIZE(X) is <M, N> .
        !          3902:       Can also be used with a multiple assignment,
        !          3903:             <M, N> = SIZE(X) .
        !          3904: 
        !          3905: SQRT  SQRT(X)  is the square root of  X .   See  FUN  .   Complex
        !          3906:       results  are  produced  if   X   is  not  positive,  or has
        !          3907:       nonpositive eigenvalues.
        !          3908: 
        !          3909: STOP  Use EXIT instead.
        !          3910: 
        !          3911: SUM   SUM(X)   is  the  sum  of  all  the  elements   of    X   .
        !          3912:       SUM(DIAG(X))  is the trace of  X .
        !          3913: 
        !          3914: SVD   Singular value decomposition.  <U,S,V> = SVD(X)  produces a
        !          3915:       diagonal  matrix  S , of the same dimension as  X  and with
        !          3916:       nonnegative diagonal  elements  in  decreasing  order,  and
        !          3917:       unitary matrices  U  and  V  so that  X = U*S*V' .
        !          3918:       By itself, SVD(X) returns a vector containing the  singular
        !          3919:       values.
        !          3920:       <U,S,V>   =   SVD(X,0)   produces   the   "economy    size"
        !          3921:       decomposition.   If  X  is m by n with m > n, then only the
        !          3922:       first n columns of U are computed and S is n by n .
        !          3923: 
        !          3924: TRIL  Lower triangle.  TRIL(X) is the lower triangular part of X.
        !          3925:       TRIL(X,K) is the elements on and below the K-th diagonal of
        !          3926:       X.  K = 0 is the main diagonal, K > 0  is  above  the  main
        !          3927:       diagonal and K < 0 is below the main diagonal.
        !          3928: 
        !          3929: TRIU  Upper triangle.  TRIU(X) is the upper triangular part of X.
        !          3930:       TRIU(X,K) is the elements on and above the K-th diagonal of
        !          3931:       X.  K = 0 is the main diagonal, K > 0  is  above  the  main
        !          3932:       diagonal and K < 0 is below the main diagonal.
        !          3933: 
        !          3934: USER  Allows personal  Fortran  subroutines  to  be  linked  into
        !          3935:       MATLAB .  The subroutine should have the heading
        !          3936: 
        !          3937:                SUBROUTINE USER(A,M,N,S,T)
        !          3938:                REAL or DOUBLE PRECISION A(M,N),S,T
        !          3939: 
        !          3940:       The MATLAB statement  Y = USER(X,s,t)  results in a call to
        !          3941:       the  subroutine with a copy of the matrix  X  stored in the
        !          3942:       argument  A , its column and row dimensions in  M  and  N ,
        !          3943:       and  the scalar parameters  s  and  t  stored in  S  and  T
        !          3944:       . If  s and t  are omitted, they are set to  0.0  .   After
        !          3945:       the  return,   A  is stored in  Y .  The dimensions  M  and
        !          3946:       N  may be reset within the subroutine.  The statement  Y  =
        !          3947:       USER(K)  results in a call with M = 1, N = 1  and  A(1,1) =
        !          3948:       FLOAT(K) .  After the subroutine has been written, it  must
        !          3949:       be compiled and linked to the MATLAB object code within the
        !          3950:       local operating system.
        !          3951: 
        !          3952: WHAT  Lists commands and functions currently available.
        !          3953: 
        !          3954: 
        !          3955: 
        !          3956: 
        !          3957: 
        !          3958: 
        !          3959: 
        !          3960: 
        !          3961: 
        !          3962: 
        !          3963: MATLAB, page 60
        !          3964: 
        !          3965: 
        !          3966: 
        !          3967: WHILE Repeat statements an indefinite number of times.
        !          3968:       WHILE expr rop expr, statement, ..., statement, END
        !          3969:       where rop is =, <, >, <=, >=, or <> (not equal) .  The  END
        !          3970:       at  the end of a line may be omitted.  The comma before the
        !          3971:       END may also be omitted.  The commas  may  be  replaced  by
        !          3972:       semicolons   to   avoid   printing.    The  statements  are
        !          3973:       repeatedly executed as long  as  the  indicated  comparison
        !          3974:       between  the  real parts of the first components of the two
        !          3975:       expressions is true.   Example  (assume  a  matrix   A   is
        !          3976:       already defined).
        !          3977:       E = 0*A; F = E + EYE; N = 1;
        !          3978:       WHILE NORM(E+F-E,1) > 0, E = E + F; F = A*F/N; N = N + 1;
        !          3979:       E
        !          3980: 
        !          3981: WHO   Lists current variables.
        !          3982: 
        !          3983: WHY   Provides succinct answers to any questions.
        !          3984: 
        !          3985: //
        !          3986: 
        !          3987: 
        !          3988: 
        !          3989: 
        !          3990: 
        !          3991: 
        !          3992: 
        !          3993: 
        !          3994: 
        !          3995: 
        !          3996: 
        !          3997: 
        !          3998: 
        !          3999: 
        !          4000: 
        !          4001: 
        !          4002: 
        !          4003: 
        !          4004: 
        !          4005: 
        !          4006: 
        !          4007: 
        !          4008: 
        !          4009: 
        !          4010: 
        !          4011: 
        !          4012: 
        !          4013: 
        !          4014: 
        !          4015: 
        !          4016: 
        !          4017: 
        !          4018: 
        !          4019: 
        !          4020: 
        !          4021: 
        !          4022: 
        !          4023: 
        !          4024: 
        !          4025: 
        !          4026: 

This archive runs on limited infrastructure. Preserving old code on modern bandwidth. Automated agents are requested to crawl responsibly.