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