This is octave.info, produced by makeinfo version 5.2 from octave.texi. START-INFO-DIR-ENTRY * Octave: (octave). Interactive language for numerical computations. END-INFO-DIR-ENTRY Copyright (C) 1996-2013 John W. Eaton. Permission is granted to make and distribute verbatim copies of this manual provided the copyright notice and this permission notice are preserved on all copies. Permission is granted to copy and distribute modified versions of this manual under the conditions for verbatim copying, provided that the entire resulting derived work is distributed under the terms of a permission notice identical to this one. Permission is granted to copy and distribute translations of this manual into another language, under the above conditions for modified versions.  File: octave.info, Node: Surface Properties, Prev: Patch Properties, Up: Graphics Object Properties 15.3.3.8 Surface Properties ........................... The 'surface' properties are: '__modified__': "off" | {"on"} 'alphadata': scalar | matrix, def. '1' Transparency is not yet implemented for surface objects. 'alphadata' is unused. 'alphadatamapping': "direct" | "none" | {"scaled"} Transparency is not yet implemented for surface objects. 'alphadatamapping' is unused. 'ambientstrength': def. '0.30000' Light is not yet implemented for surface objects. 'ambientstrength' is unused. 'backfacelighting': "lit" | {"reverselit"} | "unlit" Light is not yet implemented for surface objects. 'backfacelighting' is unused. 'beingdeleted': {"off"} | "on" 'busyaction': "cancel" | {"queue"} 'buttondownfcn': string | function handle, def. '[](0x0)' 'cdata': matrix, def. 3-by-3 double 'cdatamapping': "direct" | {"scaled"} 'cdatasource': def. "" 'children' (read-only): vector of graphics handles, def. '[](0x1)' 'children' is unused. 'clipping': "off" | {"on"} If 'clipping' is "on", the surface is clipped in its parent axes limits. 'createfcn': string | function handle, def. '[](0x0)' Callback function executed immediately after surface has been created. Function is set by using default property on root object, e.g., 'set (0, "defaultsurfacecreatefcn", 'disp ("surface created!")')'. 'deletefcn': string | function handle, def. '[](0x0)' Callback function executed immediately before surface is deleted. 'diffusestrength': def. '0.60000' Light is not yet implemented for surface objects. 'diffusestrength' is unused. 'displayname': def. "" Text for the legend entry corresponding to this surface. 'edgealpha': scalar, def. '1' Transparency is not yet implemented for surface objects. 'edgealpha' is unused. 'edgecolor': def. '[0 0 0]' 'edgelighting': "flat" | "gouraud" | {"none"} | "phong" Light is not yet implemented for surface objects. 'edgelighting' is unused. 'erasemode': "background" | "none" | {"normal"} | "xor" 'erasemode' is unused. 'facealpha': scalar | matrix, def. '1' Transparency is not yet implemented for surface objects. 'facealpha' is unused. 'facecolor': {"flat"} | "interp" | "none" | "texturemap" 'facelighting': "flat" | "gouraud" | {"none"} | "phong" Light is not yet implemented for surface objects. 'facelighting' is unused. 'handlevisibility': "callback" | "off" | {"on"} If 'handlevisibility' is "off", the surface's handle is not visible in its parent's "children" property. 'hittest': "off" | {"on"} 'interpreter': "latex" | "none" | {"tex"} 'interruptible': "off" | {"on"} 'linestyle': {"-"} | "-" | "-." | ":" | "none" *Note Line Styles::. 'linewidth': def. '0.50000' *Note line linewidth property: XREFlinelinewidth. 'marker': "*" | "+" | "." | "<" | ">" | "^" | "d" | "diamond" | "h" | "hexagram" | {"none"} | "o" | "p" | "pentagram" | "s" | "square" | "v" | "x" *Note Marker Styles::. 'markeredgecolor': {"auto"} | "flat" | "none" *Note line markeredgecolor property: XREFlinemarkeredgecolor. 'markerfacecolor': "auto" | "flat" | {"none"} *Note line markerfacecolor property: XREFlinemarkerfacecolor. 'markersize': scalar, def. '6' *Note line markersize property: XREFlinemarkersize. 'meshstyle': {"both"} | "column" | "row" 'normalmode': {"auto"} | "manual" 'parent': graphics handle Handle of the parent graphics object. 'selected': {"off"} | "on" 'selectionhighlight': "off" | {"on"} 'specularcolorreflectance': def. '1' Light is not yet implemented for surface objects. 'specularcolorreflectance' is unused. 'specularexponent': def. '10' Light is not yet implemented for surface objects. 'specularexponent' is unused. 'specularstrength': def. '0.90000' Light is not yet implemented for surface objects. 'specularstrength' is unused. 'tag': string, def. "" A user-defined string to label the graphics object. 'type' (read-only): string Class name of the graphics object. 'type' is always "surface" 'uicontextmenu': graphics handle, def. '[](0x0)' Graphics handle of the uicontextmenu object that is currently associated to this surface object. 'userdata': Any Octave data, def. '[](0x0)' User-defined data to associate with the graphics object. 'vertexnormals': def. 3-by-3-by-3 double 'visible': "off" | {"on"} If 'visible' is "off", the surface is not rendered on screen. 'xdata': matrix, def. '[1 2 3]' 'xdatasource': def. "" 'ydata': matrix, def. '[1; 2; 3]' 'ydatasource': def. "" 'zdata': matrix, def. 3-by-3 double 'zdatasource': def. ""  File: octave.info, Node: Searching Properties, Next: Managing Default Properties, Prev: Graphics Object Properties, Up: Graphics Data Structures 15.3.4 Searching Properties --------------------------- -- Function File: H = findobj () -- Function File: H = findobj (PROP_NAME, PROP_VALUE, ...) -- Function File: H = findobj (PROP_NAME, PROP_VALUE, "-LOGICAL_OP", PROP_NAME, PROP_VALUE) -- Function File: H = findobj ("-property", PROP_NAME) -- Function File: H = findobj ("-regexp", PROP_NAME, PATTERN) -- Function File: H = findobj (HLIST, ...) -- Function File: H = findobj (HLIST, "flat", ...) -- Function File: H = findobj (HLIST, "-depth", D, ...) Find graphics object with specified property values. The simplest form is findobj (PROP_NAME, PROP_VALUE) which returns the handles of all objects which have a property named PROP_NAME that has the value PROP_VALUE. If multiple property/value pairs are specified then only objects meeting all of the conditions are returned. The search can be limited to a particular set of objects and their descendants, by passing a handle or set of handles HLIST as the first argument. The depth of the object hierarchy to search can be limited with the "-depth" argument. An example of searching only three generations of children is: findobj (HLIST, "-depth", 3, PROP_NAME, PROP_VALUE) Specifying a depth D of 0, limits the search to the set of objects passed in HLIST. A depth D of 0 is equivalent to the "flat" argument. A specified logical operator may be applied to the pairs of PROP_NAME and PROP_VALUE. The supported logical operators are: "-and", "-or", "-xor", "-not". Objects may also be matched by comparing a regular expression to the property values, where property values that match 'regexp (PROP_VALUE, PATTERN)' are returned. Finally, objects may be matched by property name only by using the "-property" option. Implementation Note: The search only includes objects with visible handles (HandleVisibility = "on"). *Note findall: XREFfindall, to search for all objects including hidden ones. See also: *note findall: XREFfindall, *note allchild: XREFallchild, *note get: XREFget, *note set: XREFset. -- Function File: H = findall () -- Function File: H = findall (PROP_NAME, PROP_VALUE, ...) -- Function File: H = findall (PROP_NAME, PROP_VALUE, "-LOGICAL_OP", PROP_NAME, PROP_VALUE) -- Function File: H = findall ("-property", PROP_NAME) -- Function File: H = findall ("-regexp", PROP_NAME, PATTERN) -- Function File: H = findall (HLIST, ...) -- Function File: H = findall (HLIST, "flat", ...) -- Function File: H = findall (HLIST, "-depth", D, ...) Find graphics object, including hidden ones, with specified property values. The return value H is a list of handles to the found graphic objects. 'findall' performs the same search as 'findobj', but it includes hidden objects (HandleVisibility = "off"). For full documentation, *note findobj: XREFfindobj. See also: *note findobj: XREFfindobj, *note allchild: XREFallchild, *note get: XREFget, *note set: XREFset.  File: octave.info, Node: Managing Default Properties, Prev: Searching Properties, Up: Graphics Data Structures 15.3.5 Managing Default Properties ---------------------------------- Object properties have two classes of default values, "factory defaults" (the initial values) and "user-defined defaults", which may override the factory defaults. Although default values may be set for any object, they are set in parent objects and apply to child objects, of the specified object type. For example, setting the default 'color' property of 'line' objects to "green", for the 'root' object, will result in all 'line' objects inheriting the 'color' "green" as the default value. set (0, "defaultlinecolor", "green"); sets the default line color for all objects. The rule for constructing the property name to set a default value is default + OBJECT-TYPE + PROPERTY-NAME This rule can lead to some strange looking names, for example 'defaultlinelinewidth"' specifies the default 'linewidth' property for 'line' objects. The example above used the root figure object, 0, so the default property value will apply to all line objects. However, default values are hierarchical, so defaults set in a figure objects override those set in the root figure object. Likewise, defaults set in axes objects override those set in figure or root figure objects. For example, subplot (2, 1, 1); set (0, "defaultlinecolor", "red"); set (1, "defaultlinecolor", "green"); set (gca (), "defaultlinecolor", "blue"); line (1:10, rand (1, 10)); subplot (2, 1, 2); line (1:10, rand (1, 10)); figure (2) line (1:10, rand (1, 10)); produces two figures. The line in first subplot window of the first figure is blue because it inherits its color from its parent axes object. The line in the second subplot window of the first figure is green because it inherits its color from its parent figure object. The line in the second figure window is red because it inherits its color from the global root figure parent object. To remove a user-defined default setting, set the default property to the value "remove". For example, set (gca (), "defaultlinecolor", "remove"); removes the user-defined default line color setting from the current axes object. To quickly remove all user-defined defaults use the 'reset' function. -- Built-in Function: reset (H) Reset the properties of the graphic object H to their default values. For figures, the properties "position", "units", "windowstyle", and "paperunits" are not affected. For axes, the properties "position" and "units" are not affected. The input H may also be a vector of graphic handles in which case each individual object will be reset. See also: *note cla: XREFcla, *note clf: XREFclf, *note newplot: XREFnewplot. Getting the "default" property of an object returns a list of user-defined defaults set for the object. For example, get (gca (), "default"); returns a list of user-defined default values for the current axes object. Factory default values are stored in the root figure object. The command get (0, "factory"); returns a list of factory defaults.  File: octave.info, Node: Advanced Plotting, Prev: Graphics Data Structures, Up: Plotting 15.4 Advanced Plotting ====================== * Menu: * Colors:: * Line Styles:: * Marker Styles:: * Callbacks:: * Application-defined Data:: * Object Groups:: * Graphics Toolkits::  File: octave.info, Node: Colors, Next: Line Styles, Up: Advanced Plotting 15.4.1 Colors ------------- Colors may be specified as RGB triplets with values ranging from zero to one, or by name. Recognized color names include "blue", "black", "cyan", "green", "magenta", "red", "white", and "yellow".  File: octave.info, Node: Line Styles, Next: Marker Styles, Prev: Colors, Up: Advanced Plotting 15.4.2 Line Styles ------------------ Line styles are specified by the following properties: 'linestyle' May be one of "-" Solid line. [default] "--" Dashed line. ":" Dotted line. "-." A dash-dot line. "none" No line. Points will still be marked using the current Marker Style. 'linewidth' A number specifying the width of the line. The default is 1. A value of 2 is twice as wide as the default, etc.  File: octave.info, Node: Marker Styles, Next: Callbacks, Prev: Line Styles, Up: Advanced Plotting 15.4.3 Marker Styles -------------------- Marker styles are specified by the following properties: 'marker' A character indicating a plot marker to be place at each data point, or "none", meaning no markers should be displayed. 'markeredgecolor' The color of the edge around the marker, or "auto", meaning that the edge color is the same as the face color. *Note Colors::. 'markerfacecolor' The color of the marker, or "none" to indicate that the marker should not be filled. *Note Colors::. 'markersize' A number specifying the size of the marker. The default is 1. A value of 2 is twice as large as the default, etc. The 'colstyle' function will parse a 'plot'-style specification and will return the color, line, and marker values that would result. -- Function File: [STYLE, COLOR, MARKER, MSG] = colstyle (LINESPEC) Parse LINESPEC and return the line style, color, and markers given. In the case of an error, the string MSG will return the text of the error.  File: octave.info, Node: Callbacks, Next: Application-defined Data, Prev: Marker Styles, Up: Advanced Plotting 15.4.4 Callbacks ---------------- Callback functions can be associated with graphics objects and triggered after certain events occur. The basic structure of all callback function is function mycallback (src, data) ... endfunction where 'src' gives a handle to the source of the callback, and 'code' gives some event specific data. This can then be associated with an object either at the objects creation or later with the 'set' function. For example, plot (x, "DeleteFcn", @(s, e) disp ("Window Deleted")) where at the moment that the plot is deleted, the message "Window Deleted" will be displayed. Additional user arguments can be passed to callback functions, and will be passed after the 2 default arguments. For example: plot (x, "DeleteFcn", {@mycallback, "1"}) ... function mycallback (src, data, a1) fprintf ("Closing plot %d\n", a1); endfunction The basic callback functions that are available for all graphics objects are * CreateFcn This is the callback that is called at the moment of the objects creation. It is not called if the object is altered in any way, and so it only makes sense to define this callback in the function call that defines the object. Callbacks that are added to 'CreateFcn' later with the 'set' function will never be executed. * DeleteFcn This is the callback that is called at the moment an object is deleted. * ButtonDownFcn This is the callback that is called if a mouse button is pressed while the pointer is over this object. Note, that the gnuplot interface does not respect this callback. The object and figure that the event occurred in that resulted in the callback being called can be found with the 'gcbo' and 'gcbf' functions. -- Function File: H = gcbo () -- Function File: [H, FIG] = gcbo () Return a handle to the object whose callback is currently executing. If no callback is executing, this function returns the empty matrix. This handle is obtained from the root object property "CallbackObject". When called with a second output argument, return the handle of the figure containing the object whose callback is currently executing. If no callback is executing the second output is also set to the empty matrix. See also: *note gcbf: XREFgcbf, *note gco: XREFgco, *note gca: XREFgca, *note gcf: XREFgcf, *note get: XREFget, *note set: XREFset. -- Function File: FIG = gcbf () Return a handle to the figure containing the object whose callback is currently executing. If no callback is executing, this function returns the empty matrix. The handle returned by this function is the same as the second output argument of 'gcbo'. See also: *note gcbo: XREFgcbo, *note gcf: XREFgcf, *note gco: XREFgco, *note gca: XREFgca, *note get: XREFget, *note set: XREFset. Callbacks can equally be added to properties with the 'addlistener' function described below.  File: octave.info, Node: Application-defined Data, Next: Object Groups, Prev: Callbacks, Up: Advanced Plotting 15.4.5 Application-defined Data ------------------------------- Octave has a provision for attaching application-defined data to a graphics handle. The data can be anything which is meaningful to the application, and will be completely ignored by Octave. -- Function File: setappdata (H, NAME, VALUE) -- Function File: setappdata (H, NAME1, VALUE1, NAME2, VALUE3, ...) Set the application data NAME to VALUE for the graphics object with handle H. H may also be a vector of graphics handles. If the application data with the specified NAME does not exist, it is created. Multiple NAME/VALUE pairs can be specified at a time. See also: *note getappdata: XREFgetappdata, *note isappdata: XREFisappdata, *note rmappdata: XREFrmappdata, *note guidata: XREFguidata, *note get: XREFget, *note set: XREFset, *note getpref: XREFgetpref, *note setpref: XREFsetpref. -- Function File: VALUE = getappdata (H, NAME) -- Function File: APPDATA = getappdata (H) Return the VALUE of the application data NAME for the graphics object with handle H. H may also be a vector of graphics handles. If no second argument NAME is given then 'getappdata' returns a structure, APPDATA, whose fields correspond to the appdata properties. See also: *note setappdata: XREFsetappdata, *note isappdata: XREFisappdata, *note rmappdata: XREFrmappdata, *note guidata: XREFguidata, *note get: XREFget, *note set: XREFset, *note getpref: XREFgetpref, *note setpref: XREFsetpref. -- Function File: rmappdata (H, NAME) -- Function File: rmappdata (H, NAME1, NAME2, ...) Delete the application data NAME from the graphics object with handle H. H may also be a vector of graphics handles. Multiple application data names may be supplied to delete several properties at once. See also: *note setappdata: XREFsetappdata, *note getappdata: XREFgetappdata, *note isappdata: XREFisappdata. -- Function File: VALID = isappdata (H, NAME) Return true if the named application data, NAME, exists for the graphics object with handle H. H may also be a vector of graphics handles. See also: *note getappdata: XREFgetappdata, *note setappdata: XREFsetappdata, *note rmappdata: XREFrmappdata, *note guidata: XREFguidata, *note get: XREFget, *note set: XREFset, *note getpref: XREFgetpref, *note setpref: XREFsetpref.  File: octave.info, Node: Object Groups, Next: Graphics Toolkits, Prev: Application-defined Data, Up: Advanced Plotting 15.4.6 Object Groups -------------------- A number of Octave high level plot functions return groups of other graphics objects or they return graphics objects that have their properties linked in such a way that changes to one of the properties results in changes in the others. A graphic object that groups other objects is an 'hggroup' -- Function File: hggroup () -- Function File: hggroup (HAX) -- Function File: hggroup (..., PROPERTY, VALUE, ...) -- Function File: H = hggroup (...) Create handle graphics group object with axes parent HAX. If no parent is specified, the group is created in the current axes. Multiple property/value pairs may be specified for the hggroup, but they must appear in pairs. The optional return value H is a graphics handle to the created hggroup object. Programming Note: An hggroup is a way to group base graphics objects such as line objects or patch objects into a single unit which can react appropriately. For example, the individual lines of a contour plot are collected into a single hggroup so that they can be made visible/invisible with a single command, 'set (hg_handle, "visible", "off")'. See also: *note addproperty: XREFaddproperty, *note addlistener: XREFaddlistener. For example a simple use of a 'hggroup' might be x = 0:0.1:10; hg = hggroup (); plot (x, sin (x), "color", [1, 0, 0], "parent", hg); hold on plot (x, cos (x), "color", [0, 1, 0], "parent", hg); set (hg, "visible", "off"); which groups the two plots into a single object and controls their visibility directly. The default properties of an 'hggroup' are the same as the set of common properties for the other graphics objects. Additional properties can be added with the 'addproperty' function. -- Built-in Function: addproperty (NAME, H, TYPE) -- Built-in Function: addproperty (NAME, H, TYPE, ARG, ...) Create a new property named NAME in graphics object H. TYPE determines the type of the property to create. ARGS usually contains the default value of the property, but additional arguments might be given, depending on the type of the property. The supported property types are: 'string' A string property. ARG contains the default string value. 'any' An un-typed property. This kind of property can hold any octave value. ARGS contains the default value. 'radio' A string property with a limited set of accepted values. The first argument must be a string with all accepted values separated by a vertical bar ('|'). The default value can be marked by enclosing it with a '{' '}' pair. The default value may also be given as an optional second string argument. 'boolean' A boolean property. This property type is equivalent to a radio property with "on|off" as accepted values. ARG contains the default property value. 'double' A scalar double property. ARG contains the default value. 'handle' A handle property. This kind of property holds the handle of a graphics object. ARG contains the default handle value. When no default value is given, the property is initialized to the empty matrix. 'data' A data (matrix) property. ARG contains the default data value. When no default value is given, the data is initialized to the empty matrix. 'color' A color property. ARG contains the default color value. When no default color is given, the property is set to black. An optional second string argument may be given to specify an additional set of accepted string values (like a radio property). TYPE may also be the concatenation of a core object type and a valid property name for that object type. The property created then has the same characteristics as the referenced property (type, possible values, hidden state...). This allows one to clone an existing property into the graphics object H. Examples: addproperty ("my_property", gcf, "string", "a string value"); addproperty ("my_radio", gcf, "radio", "val_1|val_2|{val_3}"); addproperty ("my_style", gcf, "linelinestyle", "--"); See also: *note addlistener: XREFaddlistener, *note hggroup: XREFhggroup. Once a property in added to an 'hggroup', it is not linked to any other property of either the children of the group, or any other graphics object. Add so to control the way in which this newly added property is used, the 'addlistener' function is used to define a callback function that is executed when the property is altered. -- Built-in Function: addlistener (H, PROP, FCN) Register FCN as listener for the property PROP of the graphics object H. Property listeners are executed (in order of registration) when the property is set. The new value is already available when the listeners are executed. PROP must be a string naming a valid property in H. FCN can be a function handle, a string or a cell array whose first element is a function handle. If FCN is a function handle, the corresponding function should accept at least 2 arguments, that will be set to the object handle and the empty matrix respectively. If FCN is a string, it must be any valid octave expression. If FCN is a cell array, the first element must be a function handle with the same signature as described above. The next elements of the cell array are passed as additional arguments to the function. Example: function my_listener (h, dummy, p1) fprintf ("my_listener called with p1=%s\n", p1); endfunction addlistener (gcf, "position", {@my_listener, "my string"}) See also: *note addproperty: XREFaddproperty, *note hggroup: XREFhggroup. -- Built-in Function: dellistener (H, PROP, FCN) Remove the registration of FCN as a listener for the property PROP of the graphics object H. The function FCN must be the same variable (not just the same value), as was passed to the original call to 'addlistener'. If FCN is not defined then all listener functions of PROP are removed. Example: function my_listener (h, dummy, p1) fprintf ("my_listener called with p1=%s\n", p1); endfunction c = {@my_listener, "my string"}; addlistener (gcf, "position", c); dellistener (gcf, "position", c); An example of the use of these two functions might be x = 0:0.1:10; hg = hggroup (); h = plot (x, sin (x), "color", [1, 0, 0], "parent", hg); addproperty ("linestyle", hg, "linelinestyle", get (h, "linestyle")); addlistener (hg, "linestyle", @update_props); hold on plot (x, cos (x), "color", [0, 1, 0], "parent", hg); function update_props (h, d) set (get (h, "children"), "linestyle", get (h, "linestyle")); endfunction that adds a 'linestyle' property to the 'hggroup' and propagating any changes its value to the children of the group. The 'linkprop' function can be used to simplify the above to be x = 0:0.1:10; hg = hggroup (); h1 = plot (x, sin (x), "color", [1, 0, 0], "parent", hg); addproperty ("linestyle", hg, "linelinestyle", get (h, "linestyle")); hold on h2 = plot (x, cos (x), "color", [0, 1, 0], "parent", hg); hlink = linkprop ([hg, h1, h2], "color"); -- Function File: HLINK = linkprop (H, "PROP") -- Function File: HLINK = linkprop (H, {"PROP1", "PROP2", ...}) Link graphic object properties, such that a change in one is propagated to the others. The input H is a vector of graphic handles to link. PROP may be a string when linking a single property, or a cell array of strings for multiple properties. During the linking process all properties in PROP will initially be set to the values that exist on the first object in the list H. The function returns HLINK which is a special object describing the link. As long as the reference HLINK exists the link between graphic objects will be active. This means that HLINK must be preserved in a workspace variable, a global variable, or otherwise stored using a function such as 'setappdata', 'guidata'. To unlink properties, execute 'clear HLINK'. An example of the use of 'linkprop' is x = 0:0.1:10; subplot (1,2,1); h1 = plot (x, sin (x)); subplot (1,2,2); h2 = plot (x, cos (x)); hlink = linkprop ([h1, h2], {"color","linestyle"}); set (h1, "color", "green"); set (h2, "linestyle", "--"); See also: *note linkaxes: XREFlinkaxes. -- Function File: linkaxes (HAX) -- Function File: linkaxes (HAX, OPTSTR) Link the axis limits of 2-D plots such that a change in one is propagated to the others. The axes handles to be linked are passed as the first argument HAX. The optional second argument is a string which defines which axis limits will be linked. The possible values for OPTSTR are: "x" Link x-axes "y" Link y-axes "xy" (default) Link both axes "off" Turn off linking If unspecified the default is to link both X and Y axes. When linking, the limits from the first axes in HAX are applied to the other axes in the list. Subsequent changes to any one of the axes will be propagated to the others. See also: *note linkprop: XREFlinkprop, *note addproperty: XREFaddproperty. These capabilities are used in a number of basic graphics objects. The 'hggroup' objects created by the functions of Octave contain one or more graphics object and are used to: * group together multiple graphics objects, * create linked properties between different graphics objects, and * to hide the nominal user data, from the actual data of the objects. For example the 'stem' function creates a stem series where each 'hggroup' of the stem series contains two line objects representing the body and head of the stem. The 'ydata' property of the 'hggroup' of the stem series represents the head of the stem, whereas the body of the stem is between the baseline and this value. For example h = stem (1:4) get (h, "xdata") => [ 1 2 3 4]' get (get (h, "children")(1), "xdata") => [ 1 1 NaN 2 2 NaN 3 3 NaN 4 4 NaN]' shows the difference between the 'xdata' of the 'hggroup' of a stem series object and the underlying line. The basic properties of such group objects is that they consist of one or more linked 'hggroup', and that changes in certain properties of these groups are propagated to other members of the group. Whereas, certain properties of the members of the group only apply to the current member. In addition the members of the group can also be linked to other graphics objects through callback functions. For example the baseline of the 'bar' or 'stem' functions is a line object, whose length and position are automatically adjusted, based on changes to the corresponding hggroup elements. * Menu: * Data Sources in Object Groups:: * Area Series:: * Bar Series:: * Contour Groups:: * Error Bar Series:: * Line Series:: * Quiver Group:: * Scatter Group:: * Stair Group:: * Stem Series:: * Surface Group::  File: octave.info, Node: Data Sources in Object Groups, Next: Area Series, Up: Object Groups 15.4.6.1 Data Sources in Object Groups ...................................... All of the group objects contain data source parameters. There are string parameters that contain an expression that is evaluated to update the relevant data property of the group when the 'refreshdata' function is called. -- Function File: refreshdata () -- Function File: refreshdata (H) -- Function File: refreshdata (H, WORKSPACE) Evaluate any 'datasource' properties of the current figure and update the plot if the corresponding data has changed. If the first argument H is a list of graphic handles, then operate on these objects rather than the current figure returned by 'gcf'. The optional second argument WORKSPACE can take the following values: "base" Evaluate the datasource properties in the base workspace. (default). "caller" Evaluate the datasource properties in the workspace of the function that called 'refreshdata'. An example of the use of 'refreshdata' is: x = 0:0.1:10; y = sin (x); plot (x, y, "ydatasource", "y"); for i = 1 : 100 pause (0.1); y = sin (x + 0.1*i); refreshdata (); endfor  File: octave.info, Node: Area Series, Next: Bar Series, Prev: Data Sources in Object Groups, Up: Object Groups 15.4.6.2 Area Series .................... Area series objects are created by the 'area' function. Each of the 'hggroup' elements contains a single patch object. The properties of the area series are 'basevalue' The value where the base of the area plot is drawn. 'linewidth' 'linestyle' The line width and style of the edge of the patch objects making up the areas. *Note Line Styles::. 'edgecolor' 'facecolor' The line and fill color of the patch objects making up the areas. *Note Colors::. 'xdata' 'ydata' The x and y coordinates of the original columns of the data passed to 'area' prior to the cumulative summation used in the 'area' function. 'xdatasource' 'ydatasource' Data source variables.  File: octave.info, Node: Bar Series, Next: Contour Groups, Prev: Area Series, Up: Object Groups 15.4.6.3 Bar Series ................... Bar series objects are created by the 'bar' or 'barh' functions. Each 'hggroup' element contains a single patch object. The properties of the bar series are 'showbaseline' 'baseline' 'basevalue' The property 'showbaseline' flags whether the baseline of the bar series is displayed (default is "on"). The handle of the graphics object representing the baseline is given by the 'baseline' property and the y-value of the baseline by the 'basevalue' property. Changes to any of these properties are propagated to the other members of the bar series and to the baseline itself. Equally, changes in the properties of the base line itself are propagated to the members of the corresponding bar series. 'barwidth' 'barlayout' 'horizontal' The property 'barwidth' is the width of the bar corresponding to the WIDTH variable passed to 'bar' or BARH. Whether the bar series is "grouped" or "stacked" is determined by the 'barlayout' property and whether the bars are horizontal or vertical by the 'horizontal' property. Changes to any of these property are propagated to the other members of the bar series. 'linewidth' 'linestyle' The line width and style of the edge of the patch objects making up the bars. *Note Line Styles::. 'edgecolor' 'facecolor' The line and fill color of the patch objects making up the bars. *Note Colors::. 'xdata' The nominal x positions of the bars. Changes in this property and propagated to the other members of the bar series. 'ydata' The y value of the bars in the 'hggroup'. 'xdatasource' 'ydatasource' Data source variables.  File: octave.info, Node: Contour Groups, Next: Error Bar Series, Prev: Bar Series, Up: Object Groups 15.4.6.4 Contour Groups ....................... Contour group objects are created by the 'contour', 'contourf' and 'contour3' functions. The are equally one of the handles returned by the 'surfc' and 'meshc' functions. The properties of the contour group are 'contourmatrix' A read only property that contains the data return by 'contourc' used to create the contours of the plot. 'fill' A radio property that can have the values "on" or "off" that flags whether the contours to plot are to be filled. 'zlevelmode' 'zlevel' The radio property 'zlevelmode' can have the values "none", "auto", or "manual". When its value is "none" there is no z component to the plotted contours. When its value is "auto" the z value of the plotted contours is at the same value as the contour itself. If the value is "manual", then the z value at which to plot the contour is determined by the 'zlevel' property. 'levellistmode' 'levellist' 'levelstepmode' 'levelstep' If 'levellistmode' is "manual", then the levels at which to plot the contours is determined by 'levellist'. If 'levellistmode' is set to "auto", then the distance between contours is determined by 'levelstep'. If both 'levellistmode' and 'levelstepmode' are set to "auto", then there are assumed to be 10 equal spaced contours. 'textlistmode' 'textlist' 'textstepmode' 'textstep' If 'textlistmode' is "manual", then the labeled contours is determined by 'textlist'. If 'textlistmode' is set to "auto", then the distance between labeled contours is determined by 'textstep'. If both 'textlistmode' and 'textstepmode' are set to "auto", then there are assumed to be 10 equal spaced labeled contours. 'showtext' Flag whether the contour labels are shown or not. 'labelspacing' The distance between labels on a single contour in points. 'linewidth' 'linestyle' 'linecolor' The properties of the contour lines. The properties 'linewidth' and 'linestyle' are similar to the corresponding properties for lines. The property 'linecolor' is a color property (*note Colors::), that can also have the values of "none" or "auto". If 'linecolor' is "none", then no contour line is drawn. If 'linecolor' is "auto" then the line color is determined by the colormap. 'xdata' 'ydata' 'zdata' The original x, y, and z data of the contour lines. 'xdatasource' 'ydatasource' 'zdatasource' Data source variables.  File: octave.info, Node: Error Bar Series, Next: Line Series, Prev: Contour Groups, Up: Object Groups 15.4.6.5 Error Bar Series ......................... Error bar series are created by the 'errorbar' function. Each 'hggroup' element contains two line objects representing the data and the errorbars separately. The properties of the error bar series are 'color' The RGB color or color name of the line objects of the error bars. *Note Colors::. 'linewidth' 'linestyle' The line width and style of the line objects of the error bars. *Note Line Styles::. 'marker' 'markeredgecolor' 'markerfacecolor' 'markersize' The line and fill color of the markers on the error bars. *Note Colors::. 'xdata' 'ydata' 'ldata' 'udata' 'xldata' 'xudata' The original x, y, l, u, xl, xu data of the error bars. 'xdatasource' 'ydatasource' 'ldatasource' 'udatasource' 'xldatasource' 'xudatasource' Data source variables.  File: octave.info, Node: Line Series, Next: Quiver Group, Prev: Error Bar Series, Up: Object Groups 15.4.6.6 Line Series .................... Line series objects are created by the 'plot' and 'plot3' functions and are of the type 'line'. The properties of the line series with the ability to add data sources. 'color' The RGB color or color name of the line objects. *Note Colors::. 'linewidth' 'linestyle' The line width and style of the line objects. *Note Line Styles::. 'marker' 'markeredgecolor' 'markerfacecolor' 'markersize' The line and fill color of the markers. *Note Colors::. 'xdata' 'ydata' 'zdata' The original x, y and z data. 'xdatasource' 'ydatasource' 'zdatasource' Data source variables.  File: octave.info, Node: Quiver Group, Next: Scatter Group, Prev: Line Series, Up: Object Groups 15.4.6.7 Quiver Group ..................... Quiver series objects are created by the 'quiver' or 'quiver3' functions. Each 'hggroup' element of the series contains three line objects as children representing the body and head of the arrow, together with a marker as the point of origin of the arrows. The properties of the quiver series are 'autoscale' 'autoscalefactor' Flag whether the length of the arrows is scaled or defined directly from the U, V and W data. If the arrow length is flagged as being scaled by the 'autoscale' property, then the length of the autoscaled arrow is controlled by the 'autoscalefactor'. 'maxheadsize' This property controls the size of the head of the arrows in the quiver series. The default value is 0.2. 'showarrowhead' Flag whether the arrow heads are displayed in the quiver plot. 'color' The RGB color or color name of the line objects of the quiver. *Note Colors::. 'linewidth' 'linestyle' The line width and style of the line objects of the quiver. *Note Line Styles::. 'marker' 'markerfacecolor' 'markersize' The line and fill color of the marker objects at the original of the arrows. *Note Colors::. 'xdata' 'ydata' 'zdata' The origins of the values of the vector field. 'udata' 'vdata' 'wdata' The values of the vector field to plot. 'xdatasource' 'ydatasource' 'zdatasource' 'udatasource' 'vdatasource' 'wdatasource' Data source variables.  File: octave.info, Node: Scatter Group, Next: Stair Group, Prev: Quiver Group, Up: Object Groups 15.4.6.8 Scatter Group ...................... Scatter series objects are created by the 'scatter' or 'scatter3' functions. A single hggroup element contains as many children as there are points in the scatter plot, with each child representing one of the points. The properties of the stem series are 'linewidth' The line width of the line objects of the points. *Note Line Styles::. 'marker' 'markeredgecolor' 'markerfacecolor' The line and fill color of the markers of the points. *Note Colors::. 'xdata' 'ydata' 'zdata' The original x, y and z data of the stems. 'cdata' The color data for the points of the plot. Each point can have a separate color, or a unique color can be specified. 'sizedata' The size data for the points of the plot. Each point can its own size or a unique size can be specified. 'xdatasource' 'ydatasource' 'zdatasource' 'cdatasource' 'sizedatasource' Data source variables.  File: octave.info, Node: Stair Group, Next: Stem Series, Prev: Scatter Group, Up: Object Groups 15.4.6.9 Stair Group .................... Stair series objects are created by the 'stair' function. Each 'hggroup' element of the series contains a single line object as a child representing the stair. The properties of the stair series are 'color' The RGB color or color name of the line objects of the stairs. *Note Colors::. 'linewidth' 'linestyle' The line width and style of the line objects of the stairs. *Note Line Styles::. 'marker' 'markeredgecolor' 'markerfacecolor' 'markersize' The line and fill color of the markers on the stairs. *Note Colors::. 'xdata' 'ydata' The original x and y data of the stairs. 'xdatasource' 'ydatasource' Data source variables.  File: octave.info, Node: Stem Series, Next: Surface Group, Prev: Stair Group, Up: Object Groups 15.4.6.10 Stem Series ..................... Stem series objects are created by the 'stem' or 'stem3' functions. Each 'hggroup' element contains a single line object as a child representing the stems. The properties of the stem series are 'showbaseline' 'baseline' 'basevalue' The property 'showbaseline' flags whether the baseline of the stem series is displayed (default is "on"). The handle of the graphics object representing the baseline is given by the 'baseline' property and the y-value (or z-value for 'stem3') of the baseline by the 'basevalue' property. Changes to any of these property are propagated to the other members of the stem series and to the baseline itself. Equally changes in the properties of the base line itself are propagated to the members of the corresponding stem series. 'color' The RGB color or color name of the line objects of the stems. *Note Colors::. 'linewidth' 'linestyle' The line width and style of the line objects of the stems. *Note Line Styles::. 'marker' 'markeredgecolor' 'markerfacecolor' 'markersize' The line and fill color of the markers on the stems. *Note Colors::. 'xdata' 'ydata' 'zdata' The original x, y and z data of the stems. 'xdatasource' 'ydatasource' 'zdatasource' Data source variables.  File: octave.info, Node: Surface Group, Prev: Stem Series, Up: Object Groups 15.4.6.11 Surface Group ....................... Surface group objects are created by the 'surf' or 'mesh' functions, but are equally one of the handles returned by the 'surfc' or 'meshc' functions. The surface group is of the type 'surface'. The properties of the surface group are 'edgecolor' 'facecolor' The RGB color or color name of the edges or faces of the surface. *Note Colors::. 'linewidth' 'linestyle' The line width and style of the lines on the surface. *Note Line Styles::. 'marker' 'markeredgecolor' 'markerfacecolor' 'markersize' The line and fill color of the markers on the surface. *Note Colors::. 'xdata' 'ydata' 'zdata' 'cdata' The original x, y, z and c data. 'xdatasource' 'ydatasource' 'zdatasource' 'cdatasource' Data source variables.  File: octave.info, Node: Graphics Toolkits, Prev: Object Groups, Up: Advanced Plotting 15.4.7 Graphics Toolkits ------------------------ -- Function File: NAME = graphics_toolkit () -- Function File: NAME = graphics_toolkit (HLIST) -- Function File: graphics_toolkit (NAME) -- Function File: graphics_toolkit (HLIST, NAME) Query or set the default graphics toolkit which is assigned to new figures. With no inputs, return the current default graphics toolkit. If the input is a list of figure graphic handles, HLIST, then return the name of the graphics toolkit in use for each figure. When called with a single input NAME set the default graphics toolkit to NAME. If the toolkit is not already loaded, it is initialized by calling the function '__init_NAME__'. If the first input is a list of figure handles, HLIST, then the graphics toolkit is set to NAME for these figures only. See also: *note available_graphics_toolkits: XREFavailable_graphics_toolkits. -- Built-in Function: available_graphics_toolkits () Return a cell array of registered graphics toolkits. See also: *note graphics_toolkit: XREFgraphics_toolkit, *note register_graphics_toolkit: XREFregister_graphics_toolkit. -- Built-in Function: loaded_graphics_toolkits () Return a cell array of the currently loaded graphics toolkits. See also: *note available_graphics_toolkits: XREFavailable_graphics_toolkits. -- Built-in Function: register_graphics_toolkit (TOOLKIT) List TOOLKIT as an available graphics toolkit. See also: *note available_graphics_toolkits: XREFavailable_graphics_toolkits. * Menu: * Customizing Toolkit Behavior::  File: octave.info, Node: Customizing Toolkit Behavior, Up: Graphics Toolkits 15.4.7.1 Customizing Toolkit Behavior ..................................... The specific behavior of the backend toolkit may be modified using the following utility functions. Note: Not all functions apply to every graphics toolkit. -- Loadable Function: [PROG, ARGS] = gnuplot_binary () -- Loadable Function: [OLD_PROG, OLD_ARGS] = gnuplot_binary (NEW_PROG, ARG1, ...) Query or set the name of the program invoked by the plot command when the graphics toolkit is set to "gnuplot". Additional arguments to pass to the external plotting program may also be given. The default value is "gnuplot" with no additional arguments. *Note Installation::. See also: *note graphics_toolkit: XREFgraphics_toolkit.  File: octave.info, Node: Matrix Manipulation, Next: Arithmetic, Prev: Plotting, Up: Top 16 Matrix Manipulation ********************** There are a number of functions available for checking to see if the elements of a matrix meet some condition, and for rearranging the elements of a matrix. For example, Octave can easily tell you if all the elements of a matrix are finite, or are less than some specified value. Octave can also rotate the elements, extract the upper- or lower-triangular parts, or sort the columns of a matrix. * Menu: * Finding Elements and Checking Conditions:: * Rearranging Matrices:: * Special Utility Matrices:: * Famous Matrices::  File: octave.info, Node: Finding Elements and Checking Conditions, Next: Rearranging Matrices, Up: Matrix Manipulation 16.1 Finding Elements and Checking Conditions ============================================= The functions 'any' and 'all' are useful for determining whether any or all of the elements of a matrix satisfy some condition. The 'find' function is also useful in determining which elements of a matrix meet a specified condition. -- Built-in Function: any (X) -- Built-in Function: any (X, DIM) For a vector argument, return true (logical 1) if any element of the vector is nonzero. For a matrix argument, return a row vector of logical ones and zeros with each element indicating whether any of the elements of the corresponding column of the matrix are nonzero. For example: any (eye (2, 4)) => [ 1, 1, 0, 0 ] If the optional argument DIM is supplied, work along dimension DIM. For example: any (eye (2, 4), 2) => [ 1; 1 ] See also: *note all: XREFall. -- Built-in Function: all (X) -- Built-in Function: all (X, DIM) For a vector argument, return true (logical 1) if all elements of the vector are nonzero. For a matrix argument, return a row vector of logical ones and zeros with each element indicating whether all of the elements of the corresponding column of the matrix are nonzero. For example: all ([2, 3; 1, 0]) => [ 1, 0 ] If the optional argument DIM is supplied, work along dimension DIM. See also: *note any: XREFany. Since the comparison operators (*note Comparison Ops::) return matrices of ones and zeros, it is easy to test a matrix for many things, not just whether the elements are nonzero. For example, all (all (rand (5) < 0.9)) => 0 tests a random 5 by 5 matrix to see if all of its elements are less than 0.9. Note that in conditional contexts (like the test clause of 'if' and 'while' statements) Octave treats the test as if you had typed 'all (all (condition))'. -- Function File: Z = xor (X, Y) -- Function File: Z = xor (X1, X2, ...) Return the "exclusive or" of X and Y. For boolean expressions X and Y, 'xor (X, Y)' is true if and only if one of X or Y is true. Otherwise, if X and Y are both true or both false, 'xor' returns false. The truth table for the xor operation is X Y Z - - - 0 0 0 1 0 1 0 1 1 1 1 0 If more than two arguments are given the xor operation is applied cumulatively from left to right: (...((x1 XOR x2) XOR x3) XOR ...) See also: *note and: XREFand, *note or: XREFor, *note not: XREFnot. -- Built-in Function: diff (X) -- Built-in Function: diff (X, K) -- Built-in Function: diff (X, K, DIM) If X is a vector of length n, 'diff (X)' is the vector of first differences X(2) - X(1), ..., X(n) - X(n-1). If X is a matrix, 'diff (X)' is the matrix of column differences along the first non-singleton dimension. The second argument is optional. If supplied, 'diff (X, K)', where K is a non-negative integer, returns the K-th differences. It is possible that K is larger than the first non-singleton dimension of the matrix. In this case, 'diff' continues to take the differences along the next non-singleton dimension. The dimension along which to take the difference can be explicitly stated with the optional variable DIM. In this case the K-th order differences are calculated along this dimension. In the case where K exceeds 'size (X, DIM)' an empty matrix is returned. See also: *note sort: XREFsort, *note merge: XREFmerge. -- Mapping Function: isinf (X) Return a logical array which is true where the elements of X are infinite and false where they are not. For example: isinf ([13, Inf, NA, NaN]) => [ 0, 1, 0, 0 ] See also: *note isfinite: XREFisfinite, *note isnan: XREFisnan, *note isna: XREFisna. -- Mapping Function: isnan (X) Return a logical array which is true where the elements of X are NaN values and false where they are not. NA values are also considered NaN values. For example: isnan ([13, Inf, NA, NaN]) => [ 0, 0, 1, 1 ] See also: *note isna: XREFisna, *note isinf: XREFisinf, *note isfinite: XREFisfinite. -- Mapping Function: isfinite (X) Return a logical array which is true where the elements of X are finite values and false where they are not. For example: isfinite ([13, Inf, NA, NaN]) => [ 1, 0, 0, 0 ] See also: *note isinf: XREFisinf, *note isnan: XREFisnan, *note isna: XREFisna. -- Function File: [ERR, Y1, ...] = common_size (X1, ...) Determine if all input arguments are either scalar or of common size. If true, ERR is zero, and YI is a matrix of the common size with all entries equal to XI if this is a scalar or XI otherwise. If the inputs cannot be brought to a common size, ERR is 1, and YI is XI. For example: [errorcode, a, b] = common_size ([1 2; 3 4], 5) => errorcode = 0 => a = [ 1, 2; 3, 4 ] => b = [ 5, 5; 5, 5 ] This is useful for implementing functions where arguments can either be scalars or of common size. -- Built-in Function: IDX = find (X) -- Built-in Function: IDX = find (X, N) -- Built-in Function: IDX = find (X, N, DIRECTION) -- Built-in Function: [i, j] = find (...) -- Built-in Function: [i, j, v] = find (...) Return a vector of indices of nonzero elements of a matrix, as a row if X is a row vector or as a column otherwise. To obtain a single index for each matrix element, Octave pretends that the columns of a matrix form one long vector (like Fortran arrays are stored). For example: find (eye (2)) => [ 1; 4 ] If two inputs are given, N indicates the maximum number of elements to find from the beginning of the matrix or vector. If three inputs are given, DIRECTION should be one of "first" or "last", requesting only the first or last N indices, respectively. However, the indices are always returned in ascending order. If two outputs are requested, 'find' returns the row and column indices of nonzero elements of a matrix. For example: [i, j] = find (2 * eye (2)) => i = [ 1; 2 ] => j = [ 1; 2 ] If three outputs are requested, 'find' also returns a vector containing the nonzero values. For example: [i, j, v] = find (3 * eye (2)) => i = [ 1; 2 ] => j = [ 1; 2 ] => v = [ 3; 3 ] Note that this function is particularly useful for sparse matrices, as it extracts the nonzero elements as vectors, which can then be used to create the original matrix. For example: sz = size (a); [i, j, v] = find (a); b = sparse (i, j, v, sz(1), sz(2)); See also: *note nonzeros: XREFnonzeros. -- Built-in Function: IDX = lookup (TABLE, Y) -- Built-in Function: IDX = lookup (TABLE, Y, OPT) Lookup values in a sorted table. This function is usually used as a prelude to interpolation. If table is increasing and 'idx = lookup (table, y)', then 'table(idx(i)) <= y(i) < table(idx(i+1))' for all 'y(i)' within the table. If 'y(i) < table(1)' then 'idx(i)' is 0. If 'y(i) >= table(end)' or 'isnan (y(i))' then 'idx(i)' is 'n'. If the table is decreasing, then the tests are reversed. For non-strictly monotonic tables, empty intervals are always skipped. The result is undefined if TABLE is not monotonic, or if TABLE contains a NaN. The complexity of the lookup is O(M*log(N)) where N is the size of TABLE and M is the size of Y. In the special case when Y is also sorted, the complexity is O(min(M*log(N),M+N)). TABLE and Y can also be cell arrays of strings (or Y can be a single string). In this case, string lookup is performed using lexicographical comparison. If OPTS is specified, it must be a string with letters indicating additional options. 'm' 'table(idx(i)) == val(i)' if 'val(i)' occurs in table; otherwise, 'idx(i)' is zero. 'b' 'idx(i)' is a logical 1 or 0, indicating whether 'val(i)' is contained in table or not. 'l' For numeric lookups the leftmost subinterval shall be extended to infinity (i.e., all indices at least 1) 'r' For numeric lookups the rightmost subinterval shall be extended to infinity (i.e., all indices at most n-1). If you wish to check if a variable exists at all, instead of properties its elements may have, consult *note Status of Variables::.  File: octave.info, Node: Rearranging Matrices, Next: Special Utility Matrices, Prev: Finding Elements and Checking Conditions, Up: Matrix Manipulation 16.2 Rearranging Matrices ========================= -- Function File: fliplr (X) Flip array left to right. Return a copy of X with the order of the columns reversed. In other words, X is flipped left-to-right about a vertical axis. For example: fliplr ([1, 2; 3, 4]) => 2 1 4 3 See also: *note flipud: XREFflipud, *note flip: XREFflip, *note rot90: XREFrot90, *note rotdim: XREFrotdim. -- Function File: flipud (X) Flip array upside down. Return a copy of X with the order of the rows reversed. In other words, X is flipped upside-down about a horizontal axis. For example: flipud ([1, 2; 3, 4]) => 3 4 1 2 See also: *note fliplr: XREFfliplr, *note flip: XREFflip, *note rot90: XREFrot90, *note rotdim: XREFrotdim. -- Function File: flip (X) -- Function File: flip (X, DIM) Flip array across dimension DIM. Return a copy of X flipped about the dimension DIM. DIM defaults to the first non-singleton dimension. For example: flip ([1 2 3 4]) => 4 3 2 1 flip ([1; 2; 3; 4]) => 4 3 2 1 flip ([1 2; 3 4]) => 3 4 1 2 flip ([1 2; 3 4], 2) => 2 1 4 3 See also: *note fliplr: XREFfliplr, *note flipud: XREFflipud, *note rot90: XREFrot90, *note rotdim: XREFrotdim, *note permute: XREFpermute, *note transpose: XREFtranspose. -- Function File: rot90 (A) -- Function File: rot90 (A, K) Rotate array by 90 degree increments. Return a copy of A with the elements rotated counterclockwise in 90-degree increments. The second argument is optional, and specifies how many 90-degree rotations are to be applied (the default value is 1). Negative values of K rotate the matrix in a clockwise direction. For example, rot90 ([1, 2; 3, 4], -1) => 3 1 4 2 rotates the given matrix clockwise by 90 degrees. The following are all equivalent statements: rot90 ([1, 2; 3, 4], -1) rot90 ([1, 2; 3, 4], 3) rot90 ([1, 2; 3, 4], 7) The rotation is always performed on the plane of the first two dimensions, i.e., rows and columns. To perform a rotation on any other plane, use 'rotdim'. See also: *note rotdim: XREFrotdim, *note fliplr: XREFfliplr, *note flipud: XREFflipud, *note flip: XREFflip. -- Function File: rotdim (X) -- Function File: rotdim (X, N) -- Function File: rotdim (X, N, PLANE) Return a copy of X with the elements rotated counterclockwise in 90-degree increments. The second argument N is optional, and specifies how many 90-degree rotations are to be applied (the default value is 1). Negative values of N rotate the matrix in a clockwise direction. The third argument is also optional and defines the plane of the rotation. If present, PLANE is a two element vector containing two different valid dimensions of the matrix. When PLANE is not given the first two non-singleton dimensions are used. For example, rotdim ([1, 2; 3, 4], -1, [1, 2]) => 3 1 4 2 rotates the given matrix clockwise by 90 degrees. The following are all equivalent statements: rotdim ([1, 2; 3, 4], -1, [1, 2]) rotdim ([1, 2; 3, 4], 3, [1, 2]) rotdim ([1, 2; 3, 4], 7, [1, 2]) See also: *note rot90: XREFrot90, *note fliplr: XREFfliplr, *note flipud: XREFflipud, *note flip: XREFflip. -- Built-in Function: cat (DIM, ARRAY1, ARRAY2, ..., ARRAYN) Return the concatenation of N-D array objects, ARRAY1, ARRAY2, ..., ARRAYN along dimension DIM. A = ones (2, 2); B = zeros (2, 2); cat (2, A, B) => 1 1 0 0 1 1 0 0 Alternatively, we can concatenate A and B along the second dimension in the following way: [A, B] DIM can be larger than the dimensions of the N-D array objects and the result will thus have DIM dimensions as the following example shows: cat (4, ones (2, 2), zeros (2, 2)) => ans(:,:,1,1) = 1 1 1 1 ans(:,:,1,2) = 0 0 0 0 See also: *note horzcat: XREFhorzcat, *note vertcat: XREFvertcat. -- Built-in Function: horzcat (ARRAY1, ARRAY2, ..., ARRAYN) Return the horizontal concatenation of N-D array objects, ARRAY1, ARRAY2, ..., ARRAYN along dimension 2. Arrays may also be concatenated horizontally using the syntax for creating new matrices. For example: HCAT = [ ARRAY1, ARRAY2, ... ] See also: *note cat: XREFcat, *note vertcat: XREFvertcat. -- Built-in Function: vertcat (ARRAY1, ARRAY2, ..., ARRAYN) Return the vertical concatenation of N-D array objects, ARRAY1, ARRAY2, ..., ARRAYN along dimension 1. Arrays may also be concatenated vertically using the syntax for creating new matrices. For example: VCAT = [ ARRAY1; ARRAY2; ... ] See also: *note cat: XREFcat, *note horzcat: XREFhorzcat. -- Built-in Function: permute (A, PERM) Return the generalized transpose for an N-D array object A. The permutation vector PERM must contain the elements '1:ndims (A)' (in any order, but each element must appear only once). The Nth dimension of A gets remapped to dimension 'PERM(N)'. For example: X = zeros ([2, 3, 5, 7]); size (X) => 2 3 5 7 size (permute (X, [2, 1, 3, 4])) => 3 2 5 7 size (permute (X, [1, 3, 4, 2])) => 2 5 7 3 ## The identity permutation size (permute (X, [1, 2, 3, 4])) => 2 3 5 7 See also: *note ipermute: XREFipermute. -- Built-in Function: ipermute (A, IPERM) The inverse of the 'permute' function. The expression ipermute (permute (A, perm), perm) returns the original array A. See also: *note permute: XREFpermute. -- Built-in Function: reshape (A, M, N, ...) -- Built-in Function: reshape (A, [M N ...]) -- Built-in Function: reshape (A, ..., [], ...) -- Built-in Function: reshape (A, SIZE) Return a matrix with the specified dimensions (M, N, ...) whose elements are taken from the matrix A. The elements of the matrix are accessed in column-major order (like Fortran arrays are stored). The following code demonstrates reshaping a 1x4 row vector into a 2x2 square matrix. reshape ([1, 2, 3, 4], 2, 2) => 1 3 2 4 Note that the total number of elements in the original matrix ('prod (size (A))') must match the total number of elements in the new matrix ('prod ([M N ...])'). A single dimension of the return matrix may be left unspecified and Octave will determine its size automatically. An empty matrix ([]) is used to flag the unspecified dimension. See also: *note resize: XREFresize, *note vec: XREFvec, *note postpad: XREFpostpad, *note cat: XREFcat, *note squeeze: XREFsqueeze. -- Built-in Function: resize (X, M) -- Built-in Function: resize (X, M, N, ...) -- Built-in Function: resize (X, [M N ...]) Resize X cutting off elements as necessary. In the result, element with certain indices is equal to the corresponding element of X if the indices are within the bounds of X; otherwise, the element is set to zero. In other words, the statement y = resize (x, dv) is equivalent to the following code: y = zeros (dv, class (x)); sz = min (dv, size (x)); for i = 1:length (sz) idx{i} = 1:sz(i); endfor y(idx{:}) = x(idx{:}); but is performed more efficiently. If only M is supplied, and it is a scalar, the dimension of the result is M-by-M. If M, N, ... are all scalars, then the dimensions of the result are M-by-N-by-.... If given a vector as input, then the dimensions of the result are given by the elements of that vector. An object can be resized to more dimensions than it has; in such case the missing dimensions are assumed to be 1. Resizing an object to fewer dimensions is not possible. See also: *note reshape: XREFreshape, *note postpad: XREFpostpad, *note prepad: XREFprepad, *note cat: XREFcat. -- Function File: Y = circshift (X, N) Circularly shift the values of the array X. N must be a vector of integers no longer than the number of dimensions in X. The values of N can be either positive or negative, which determines the direction in which the values or X are shifted. If an element of N is zero, then the corresponding dimension of X will not be shifted. For example: x = [1, 2, 3; 4, 5, 6; 7, 8, 9]; circshift (x, 1) => 7, 8, 9 1, 2, 3 4, 5, 6 circshift (x, -2) => 7, 8, 9 1, 2, 3 4, 5, 6 circshift (x, [0,1]) => 3, 1, 2 6, 4, 5 9, 7, 8 See also: *note permute: XREFpermute, *note ipermute: XREFipermute, *note shiftdim: XREFshiftdim. -- Function File: shift (X, B) -- Function File: shift (X, B, DIM) If X is a vector, perform a circular shift of length B of the elements of X. If X is a matrix, do the same for each column of X. If the optional DIM argument is given, operate along this dimension. -- Function File: Y = shiftdim (X, N) -- Function File: [Y, NS] = shiftdim (X) Shift the dimensions of X by N, where N must be an integer scalar. When N is positive, the dimensions of X are shifted to the left, with the leading dimensions circulated to the end. If N is negative, then the dimensions of X are shifted to the right, with N leading singleton dimensions added. Called with a single argument, 'shiftdim', removes the leading singleton dimensions, returning the number of dimensions removed in the second output argument NS. For example: x = ones (1, 2, 3); size (shiftdim (x, -1)) => [1, 1, 2, 3] size (shiftdim (x, 1)) => [2, 3] [b, ns] = shiftdim (x) => b = [1, 1, 1; 1, 1, 1] => ns = 1 See also: *note reshape: XREFreshape, *note permute: XREFpermute, *note ipermute: XREFipermute, *note circshift: XREFcircshift, *note squeeze: XREFsqueeze. -- Built-in Function: [S, I] = sort (X) -- Built-in Function: [S, I] = sort (X, DIM) -- Built-in Function: [S, I] = sort (X, MODE) -- Built-in Function: [S, I] = sort (X, DIM, MODE) Return a copy of X with the elements arranged in increasing order. For matrices, 'sort' orders the elements within columns For example: sort ([1, 2; 2, 3; 3, 1]) => 1 1 2 2 3 3 If the optional argument DIM is given, then the matrix is sorted along the dimension defined by DIM. The optional argument 'mode' defines the order in which the values will be sorted. Valid values of 'mode' are "ascend" or "descend". The 'sort' function may also be used to produce a matrix containing the original row indices of the elements in the sorted matrix. For example: [s, i] = sort ([1, 2; 2, 3; 3, 1]) => s = 1 1 2 2 3 3 => i = 1 3 2 1 3 2 For equal elements, the indices are such that equal elements are listed in the order in which they appeared in the original list. Sorting of complex entries is done first by magnitude ('abs (Z)') and for any ties by phase angle ('angle (z)'). For example: sort ([1+i; 1; 1-i]) => 1 + 0i 1 - 1i 1 + 1i NaN values are treated as being greater than any other value and are sorted to the end of the list. The 'sort' function may also be used to sort strings and cell arrays of strings, in which case ASCII dictionary order (uppercase 'A' precedes lowercase 'a') of the strings is used. The algorithm used in 'sort' is optimized for the sorting of partially ordered lists. See also: *note sortrows: XREFsortrows, *note issorted: XREFissorted. -- Function File: [S, I] = sortrows (A) -- Function File: [S, I] = sortrows (A, C) Sort the rows of the matrix A according to the order of the columns specified in C. If C is omitted, a lexicographical sort is used. By default ascending order is used however if elements of C are negative then the corresponding column is sorted in descending order. See also: *note sort: XREFsort. -- Built-in Function: issorted (A) -- Built-in Function: issorted (A, MODE) -- Built-in Function: issorted (A, "rows", MODE) Return true if the array is sorted according to MODE, which may be either "ascending", "descending", or "either". By default, MODE is "ascending". NaNs are treated in the same manner as 'sort'. If the optional argument "rows" is supplied, check whether the array is sorted by rows as output by the function 'sortrows' (with no options). This function does not support sparse matrices. See also: *note sort: XREFsort, *note sortrows: XREFsortrows. -- Built-in Function: nth_element (X, N) -- Built-in Function: nth_element (X, N, DIM) Select the n-th smallest element of a vector, using the ordering defined by 'sort'. The result is equivalent to 'sort(X)(N)'. N can also be a contiguous range, either ascending 'l:u' or descending 'u:-1:l', in which case a range of elements is returned. If X is an array, 'nth_element' operates along the dimension defined by DIM, or the first non-singleton dimension if DIM is not given. Programming Note: nth_element encapsulates the C++ standard library algorithms nth_element and partial_sort. On average, the complexity of the operation is O(M*log(K)), where 'M = size (X, DIM)' and 'K = length (N)'. This function is intended for cases where the ratio K/M is small; otherwise, it may be better to use 'sort'. See also: *note sort: XREFsort, *note min: XREFmin, *note max: XREFmax. -- Function File: tril (A) -- Function File: tril (A, K) -- Function File: tril (A, K, PACK) -- Function File: triu (A) -- Function File: triu (A, K) -- Function File: triu (A, K, PACK) Return a new matrix formed by extracting the lower ('tril') or upper ('triu') triangular part of the matrix A, and setting all other elements to zero. The second argument is optional, and specifies how many diagonals above or below the main diagonal should also be set to zero. The default value of K is zero, so that 'triu' and 'tril' normally include the main diagonal as part of the result. If the value of K is nonzero integer, the selection of elements starts at an offset of K diagonals above or below the main diagonal; above for positive K and below for negative K. The absolute value of K must not be greater than the number of subdiagonals or superdiagonals. For example: tril (ones (3), -1) => 0 0 0 1 0 0 1 1 0 and tril (ones (3), 1) => 1 1 0 1 1 1 1 1 1 If the option "pack" is given as third argument, the extracted elements are not inserted into a matrix, but rather stacked column-wise one above other. See also: *note diag: XREFdiag. -- Built-in Function: V = vec (X) -- Built-in Function: V = vec (X, DIM) Return the vector obtained by stacking the columns of the matrix X one above the other. Without DIM this is equivalent to 'X(:)'. If DIM is supplied, the dimensions of V are set to DIM with all elements along the last dimension. This is equivalent to 'shiftdim (X(:), 1-DIM)'. See also: *note vech: XREFvech, *note resize: XREFresize, *note cat: XREFcat. -- Function File: vech (X) Return the vector obtained by eliminating all superdiagonal elements of the square matrix X and stacking the result one column above the other. This has uses in matrix calculus where the underlying matrix is symmetric and it would be pointless to keep values above the main diagonal. See also: *note vec: XREFvec. -- Function File: prepad (X, L) -- Function File: prepad (X, L, C) -- Function File: prepad (X, L, C, DIM) Prepend the scalar value C to the vector X until it is of length L. If C is not given, a value of 0 is used. If 'length (X) > L', elements from the beginning of X are removed until a vector of length L is obtained. If X is a matrix, elements are prepended or removed from each row. If the optional argument DIM is given, operate along this dimension. If DIM is larger than the dimensions of X, the result will have DIM dimensions. See also: *note postpad: XREFpostpad, *note cat: XREFcat, *note resize: XREFresize. -- Function File: postpad (X, L) -- Function File: postpad (X, L, C) -- Function File: postpad (X, L, C, DIM) Append the scalar value C to the vector X until it is of length L. If C is not given, a value of 0 is used. If 'length (X) > L', elements from the end of X are removed until a vector of length L is obtained. If X is a matrix, elements are appended or removed from each row. If the optional argument DIM is given, operate along this dimension. If DIM is larger than the dimensions of X, the result will have DIM dimensions. See also: *note prepad: XREFprepad, *note cat: XREFcat, *note resize: XREFresize. -- Built-in Function: M = diag (V) -- Built-in Function: M = diag (V, K) -- Built-in Function: M = diag (V, M, N) -- Built-in Function: V = diag (M) -- Built-in Function: V = diag (M, K) Return a diagonal matrix with vector V on diagonal K. The second argument is optional. If it is positive, the vector is placed on the K-th superdiagonal. If it is negative, it is placed on the -K-th subdiagonal. The default value of K is 0, and the vector is placed on the main diagonal. For example: diag ([1, 2, 3], 1) => 0 1 0 0 0 0 2 0 0 0 0 3 0 0 0 0 The 3-input form returns a diagonal matrix with vector V on the main diagonal and the resulting matrix being of size M rows x N columns. Given a matrix argument, instead of a vector, 'diag' extracts the K-th diagonal of the matrix. -- Function File: blkdiag (A, B, C, ...) Build a block diagonal matrix from A, B, C, ... All arguments must be numeric and either two-dimensional matrices or scalars. If any argument is of type sparse, the output will also be sparse. See also: *note diag: XREFdiag, *note horzcat: XREFhorzcat, *note vertcat: XREFvertcat, *note sparse: XREFsparse.  File: octave.info, Node: Special Utility Matrices, Next: Famous Matrices, Prev: Rearranging Matrices, Up: Matrix Manipulation 16.3 Special Utility Matrices ============================= -- Built-in Function: eye (N) -- Built-in Function: eye (M, N) -- Built-in Function: eye ([M N]) -- Built-in Function: eye (..., CLASS) Return an identity matrix. If invoked with a single scalar argument N, return a square NxN identity matrix. If supplied two scalar arguments (M, N), 'eye' takes them to be the number of rows and columns. If given a vector with two elements, 'eye' uses the values of the elements as the number of rows and columns, respectively. For example: eye (3) => 1 0 0 0 1 0 0 0 1 The following expressions all produce the same result: eye (2) == eye (2, 2) == eye (size ([1, 2; 3, 4])) The optional argument CLASS, allows 'eye' to return an array of the specified type, like val = zeros (n,m, "uint8") Calling 'eye' with no arguments is equivalent to calling it with an argument of 1. Any negative dimensions are treated as zero. These odd definitions are for compatibility with MATLAB. See also: *note speye: XREFspeye, *note ones: XREFones, *note zeros: XREFzeros. -- Built-in Function: ones (N) -- Built-in Function: ones (M, N) -- Built-in Function: ones (M, N, K, ...) -- Built-in Function: ones ([M N ...]) -- Built-in Function: ones (..., CLASS) Return a matrix or N-dimensional array whose elements are all 1. If invoked with a single scalar integer argument N, return a square NxN matrix. If invoked with two or more scalar integer arguments, or a vector of integer values, return an array with the given dimensions. To create a constant matrix whose values are all the same use an expression such as val_matrix = val * ones (m, n) The optional argument CLASS specifies the class of the return array and defaults to double. For example: val = ones (m,n, "uint8") See also: *note zeros: XREFzeros. -- Built-in Function: zeros (N) -- Built-in Function: zeros (M, N) -- Built-in Function: zeros (M, N, K, ...) -- Built-in Function: zeros ([M N ...]) -- Built-in Function: zeros (..., CLASS) Return a matrix or N-dimensional array whose elements are all 0. If invoked with a single scalar integer argument, return a square NxN matrix. If invoked with two or more scalar integer arguments, or a vector of integer values, return an array with the given dimensions. The optional argument CLASS specifies the class of the return array and defaults to double. For example: val = zeros (m,n, "uint8") See also: *note ones: XREFones. -- Function File: repmat (A, M) -- Function File: repmat (A, M, N) -- Function File: repmat (A, M, N, P ...) -- Function File: repmat (A, [M N]) -- Function File: repmat (A, [M N P ...]) Form a block matrix of size M by N, with a copy of matrix A as each element. If N is not specified, form an M by M block matrix. For copying along more than two dimensions, specify the number of times to copy across each dimension M, N, P, ..., in a vector in the second argument. See also: *note repelems: XREFrepelems. -- Built-in Function: repelems (X, R) Construct a vector of repeated elements from X. R is a 2xN integer matrix specifying which elements to repeat and how often to repeat each element. Entries in the first row, R(1,j), select an element to repeat. The corresponding entry in the second row, R(2,j), specifies the repeat count. If X is a matrix then the columns of X are imagined to be stacked on top of each other for purposes of the selection index. A row vector is always returned. Conceptually the result is calculated as follows: y = []; for i = 1:columns (R) y = [y, X(R(1,i)*ones(1, R(2,i)))]; endfor See also: *note repmat: XREFrepmat, *note cat: XREFcat. The functions 'linspace' and 'logspace' make it very easy to create vectors with evenly or logarithmically spaced elements. *Note Ranges::. -- Built-in Function: linspace (BASE, LIMIT) -- Built-in Function: linspace (BASE, LIMIT, N) Return a row vector with N linearly spaced elements between BASE and LIMIT. If the number of elements is greater than one, then the endpoints BASE and LIMIT are always included in the range. If BASE is greater than LIMIT, the elements are stored in decreasing order. If the number of points is not specified, a value of 100 is used. The 'linspace' function always returns a row vector if both BASE and LIMIT are scalars. If one, or both, of them are column vectors, 'linspace' returns a matrix. For compatibility with MATLAB, return the second argument (LIMIT) if fewer than two values are requested. See also: *note logspace: XREFlogspace. -- Function File: logspace (A, B) -- Function File: logspace (A, B, N) -- Function File: logspace (A, pi, N) Return a row vector with N elements logarithmically spaced from 10^A to 10^B. If N is unspecified it defaults to 50. If B is equal to pi, the points are between 10^A and pi, _not_ 10^A and 10^pi, in order to be compatible with the corresponding MATLAB function. Also for compatibility with MATLAB, return the second argument B if fewer than two values are requested. See also: *note linspace: XREFlinspace. -- Built-in Function: rand (N) -- Built-in Function: rand (M, N, ...) -- Built-in Function: rand ([M N ...]) -- Built-in Function: V = rand ("state") -- Built-in Function: rand ("state", V) -- Built-in Function: rand ("state", "reset") -- Built-in Function: V = rand ("seed") -- Built-in Function: rand ("seed", V) -- Built-in Function: rand ("seed", "reset") -- Built-in Function: rand (..., "single") -- Built-in Function: rand (..., "double") Return a matrix with random elements uniformly distributed on the interval (0, 1). The arguments are handled the same as the arguments for 'eye'. You can query the state of the random number generator using the form v = rand ("state") This returns a column vector V of length 625. Later, you can restore the random number generator to the state V using the form rand ("state", v) You may also initialize the state vector from an arbitrary vector of length <= 625 for V. This new state will be a hash based on the value of V, not V itself. By default, the generator is initialized from '/dev/urandom' if it is available, otherwise from CPU time, wall clock time, and the current fraction of a second. Note that this differs from MATLAB, which always initializes the state to the same state at startup. To obtain behavior comparable to MATLAB, initialize with a deterministic state vector in Octave's startup files (*note Startup Files::). To compute the pseudo-random sequence, 'rand' uses the Mersenne Twister with a period of 2^{19937}-1 (See M. Matsumoto and T. Nishimura, 'Mersenne Twister: A 623-dimensionally equidistributed uniform pseudorandom number generator', ACM Trans. on Modeling and Computer Simulation Vol. 8, No. 1, pp. 3-30, January 1998, ). Do *not* use for cryptography without securely hashing several returned values together, otherwise the generator state can be learned after reading 624 consecutive values. Older versions of Octave used a different random number generator. The new generator is used by default as it is significantly faster than the old generator, and produces random numbers with a significantly longer cycle time. However, in some circumstances it might be desirable to obtain the same random sequences as produced by the old generators. To do this the keyword "seed" is used to specify that the old generators should be used, as in rand ("seed", val) which sets the seed of the generator to VAL. The seed of the generator can be queried with s = rand ("seed") However, it should be noted that querying the seed will not cause 'rand' to use the old generators, only setting the seed will. To cause 'rand' to once again use the new generators, the keyword "state" should be used to reset the state of the 'rand'. The state or seed of the generator can be reset to a new random value using the "reset" keyword. The class of the value returned can be controlled by a trailing "double" or "single" argument. These are the only valid classes. See also: *note randn: XREFrandn, *note rande: XREFrande, *note randg: XREFrandg, *note randp: XREFrandp. -- Function File: randi (IMAX) -- Function File: randi (IMAX, N) -- Function File: randi (IMAX, M, N, ...) -- Function File: randi ([IMIN IMAX], ...) -- Function File: randi (..., "CLASS") Return random integers in the range 1:IMAX. Additional arguments determine the shape of the return matrix. When no arguments are specified a single random integer is returned. If one argument N is specified then a square matrix (N x N) is returned. Two or more arguments will return a multi-dimensional matrix (M x N x ...). The integer range may optionally be described by a two element matrix with a lower and upper bound in which case the returned integers will be on the interval [IMIN, IMAX]. The optional argument CLASS will return a matrix of the requested type. The default is "double". The following example returns 150 integers in the range 1-10. ri = randi (10, 150, 1) Implementation Note: 'randi' relies internally on 'rand' which uses class "double" to represent numbers. This limits the maximum integer (IMAX) and range (IMAX - IMIN) to the value returned by the 'bitmax' function. For IEEE floating point numbers this value is 2^{53} - 1. See also: *note rand: XREFrand. -- Built-in Function: randn (N) -- Built-in Function: randn (M, N, ...) -- Built-in Function: randn ([M N ...]) -- Built-in Function: V = randn ("state") -- Built-in Function: randn ("state", V) -- Built-in Function: randn ("state", "reset") -- Built-in Function: V = randn ("seed") -- Built-in Function: randn ("seed", V) -- Built-in Function: randn ("seed", "reset") -- Built-in Function: randn (..., "single") -- Built-in Function: randn (..., "double") Return a matrix with normally distributed random elements having zero mean and variance one. The arguments are handled the same as the arguments for 'rand'. By default, 'randn' uses the Marsaglia and Tsang "Ziggurat technique" to transform from a uniform to a normal distribution. The class of the value returned can be controlled by a trailing "double" or "single" argument. These are the only valid classes. Reference: G. Marsaglia and W.W. Tsang, 'Ziggurat Method for Generating Random Variables', J. Statistical Software, vol 5, 2000, See also: *note rand: XREFrand, *note rande: XREFrande, *note randg: XREFrandg, *note randp: XREFrandp. -- Built-in Function: rande (N) -- Built-in Function: rande (M, N, ...) -- Built-in Function: rande ([M N ...]) -- Built-in Function: V = rande ("state") -- Built-in Function: rande ("state", V) -- Built-in Function: rande ("state", "reset") -- Built-in Function: V = rande ("seed") -- Built-in Function: rande ("seed", V) -- Built-in Function: rande ("seed", "reset") -- Built-in Function: rande (..., "single") -- Built-in Function: rande (..., "double") Return a matrix with exponentially distributed random elements. The arguments are handled the same as the arguments for 'rand'. By default, 'randn' uses the Marsaglia and Tsang "Ziggurat technique" to transform from a uniform to a normal distribution. The class of the value returned can be controlled by a trailing "double" or "single" argument. These are the only valid classes. Reference: G. Marsaglia and W.W. Tsang, 'Ziggurat Method for Generating Random Variables', J. Statistical Software, vol 5, 2000, See also: *note rand: XREFrand, *note randn: XREFrandn, *note randg: XREFrandg, *note randp: XREFrandp. -- Built-in Function: randp (L, N) -- Built-in Function: randp (L, M, N, ...) -- Built-in Function: randp (L, [M N ...]) -- Built-in Function: V = randp ("state") -- Built-in Function: randp ("state", V) -- Built-in Function: randp ("state", "reset") -- Built-in Function: V = randp ("seed") -- Built-in Function: randp ("seed", V) -- Built-in Function: randp ("seed", "reset") -- Built-in Function: randp (..., "single") -- Built-in Function: randp (..., "double") Return a matrix with Poisson distributed random elements with mean value parameter given by the first argument, L. The arguments are handled the same as the arguments for 'rand', except for the argument L. Five different algorithms are used depending on the range of L and whether or not L is a scalar or a matrix. For scalar L <= 12, use direct method. W.H. Press, et al., 'Numerical Recipes in C', Cambridge University Press, 1992. For scalar L > 12, use rejection method.[1] W.H. Press, et al., 'Numerical Recipes in C', Cambridge University Press, 1992. For matrix L <= 10, use inversion method.[2] E. Stadlober, et al., WinRand source code, available via FTP. For matrix L > 10, use patchwork rejection method. E. Stadlober, et al., WinRand source code, available via FTP, or H. Zechner, 'Efficient sampling from continuous and discrete unimodal distributions', Doctoral Dissertation, 156pp., Technical University Graz, Austria, 1994. For L > 1e8, use normal approximation. L. Montanet, et al., 'Review of Particle Properties', Physical Review D 50 p1284, 1994. The class of the value returned can be controlled by a trailing "double" or "single" argument. These are the only valid classes. See also: *note rand: XREFrand, *note randn: XREFrandn, *note rande: XREFrande, *note randg: XREFrandg. -- Built-in Function: randg (N) -- Built-in Function: randg (M, N, ...) -- Built-in Function: randg ([M N ...]) -- Built-in Function: V = randg ("state") -- Built-in Function: randg ("state", V) -- Built-in Function: randg ("state", "reset") -- Built-in Function: V = randg ("seed") -- Built-in Function: randg ("seed", V) -- Built-in Function: randg ("seed", "reset") -- Built-in Function: randg (..., "single") -- Built-in Function: randg (..., "double") Return a matrix with 'gamma (A,1)' distributed random elements. The arguments are handled the same as the arguments for 'rand', except for the argument A. This can be used to generate many distributions: 'gamma (a, b)' for 'a > -1', 'b > 0' r = b * randg (a) 'beta (a, b)' for 'a > -1', 'b > -1' r1 = randg (a, 1) r = r1 / (r1 + randg (b, 1)) 'Erlang (a, n)' r = a * randg (n) 'chisq (df)' for 'df > 0' r = 2 * randg (df / 2) 't (df)' for '0 < df < inf' (use randn if df is infinite) r = randn () / sqrt (2 * randg (df / 2) / df) 'F (n1, n2)' for '0 < n1', '0 < n2' ## r1 equals 1 if n1 is infinite r1 = 2 * randg (n1 / 2) / n1 ## r2 equals 1 if n2 is infinite r2 = 2 * randg (n2 / 2) / n2 r = r1 / r2 negative 'binomial (n, p)' for 'n > 0', '0 < p <= 1' r = randp ((1 - p) / p * randg (n)) non-central 'chisq (df, L)', for 'df >= 0' and 'L > 0' (use chisq if 'L = 0') r = randp (L / 2) r(r > 0) = 2 * randg (r(r > 0)) r(df > 0) += 2 * randg (df(df > 0)/2) 'Dirichlet (a1, ... ak)' r = (randg (a1), ..., randg (ak)) r = r / sum (r) The class of the value returned can be controlled by a trailing "double" or "single" argument. These are the only valid classes. See also: *note rand: XREFrand, *note randn: XREFrandn, *note rande: XREFrande, *note randp: XREFrandp. The generators operate in the new or old style together, it is not possible to mix the two. Initializing any generator with "state" or "seed" causes the others to switch to the same style for future calls. The state of each generator is independent and calls to different generators can be interleaved without affecting the final result. For example, rand ("state", [11, 22, 33]); randn ("state", [44, 55, 66]); u = rand (100, 1); n = randn (100, 1); and rand ("state", [11, 22, 33]); randn ("state", [44, 55, 66]); u = zeros (100, 1); n = zeros (100, 1); for i = 1:100 u(i) = rand (); n(i) = randn (); end produce equivalent results. When the generators are initialized in the old style with "seed" only 'rand' and 'randn' are independent, because the old 'rande', 'randg' and 'randp' generators make calls to 'rand' and 'randn'. The generators are initialized with random states at start-up, so that the sequences of random numbers are not the same each time you run Octave.(1) If you really do need to reproduce a sequence of numbers exactly, you can set the state or seed to a specific value. If invoked without arguments, 'rand' and 'randn' return a single element of a random sequence. The original 'rand' and 'randn' functions use Fortran code from RANLIB, a library of Fortran routines for random number generation, compiled by Barry W. Brown and James Lovato of the Department of Biomathematics at The University of Texas, M.D. Anderson Cancer Center, Houston, TX 77030. -- Built-in Function: randperm (N) -- Built-in Function: randperm (N, M) Return a row vector containing a random permutation of '1:N'. If M is supplied, return M unique entries, sampled without replacement from '1:N'. The complexity is O(N) in memory and O(M) in time, unless M < N/5, in which case O(M) memory is used as well. The randomization is performed using rand(). All permutations are equally likely. See also: *note perms: XREFperms. ---------- Footnotes ---------- (1) The old versions of 'rand' and 'randn' obtain their initial seeds from the system clock.  File: octave.info, Node: Famous Matrices, Prev: Special Utility Matrices, Up: Matrix Manipulation 16.4 Famous Matrices ==================== The following functions return famous matrix forms. -- Function File: gallery (NAME) -- Function File: gallery (NAME, ARGS) Create interesting matrices for testing. -- Function File: C = gallery ("cauchy", X) -- Function File: C = gallery ("cauchy", X, Y) Create a Cauchy matrix. -- Function File: C = gallery ("chebspec", N) -- Function File: C = gallery ("chebspec", N, K) Create a Chebyshev spectral differentiation matrix. -- Function File: C = gallery ("chebvand", P) -- Function File: C = gallery ("chebvand", M, P) Create a Vandermonde-like matrix for the Chebyshev polynomials. -- Function File: A = gallery ("chow", N) -- Function File: A = gallery ("chow", N, ALPHA) -- Function File: A = gallery ("chow", N, ALPHA, DELTA) Create a Chow matrix - a singular Toeplitz lower Hessenberg matrix. -- Function File: C = gallery ("circul", V) Create a circulant matrix. -- Function File: A = gallery ("clement", N) -- Function File: A = gallery ("clement", N, K) Create a tridiagonal matrix with zero diagonal entries. -- Function File: C = gallery ("compar", A) -- Function File: C = gallery ("compar", A, K) Create a comparison matrix. -- Function File: A = gallery ("condex", N) -- Function File: A = gallery ("condex", N, K) -- Function File: A = gallery ("condex", N, K, THETA) Create a 'counterexample' matrix to a condition estimator. -- Function File: A = gallery ("cycol", [M N]) -- Function File: A = gallery ("cycol", N) -- Function File: A = gallery (..., K) Create a matrix whose columns repeat cyclically. -- Function File: [C, D, E] = gallery ("dorr", N) -- Function File: [C, D, E] = gallery ("dorr", N, THETA) -- Function File: A = gallery ("dorr", ...) Create a diagonally dominant, ill-conditioned, tridiagonal matrix. -- Function File: A = gallery ("dramadah", N) -- Function File: A = gallery ("dramadah", N, K) Create a (0, 1) matrix whose inverse has large integer entries. -- Function File: A = gallery ("fiedler", C) Create a symmetric Fiedler matrix. -- Function File: A = gallery ("forsythe", N) -- Function File: A = gallery ("forsythe", N, ALPHA) -- Function File: A = gallery ("forsythe", N, ALPHA, LAMBDA) Create a Forsythe matrix (a perturbed Jordan block). -- Function File: F = gallery ("frank", N) -- Function File: F = gallery ("frank", N, K) Create a Frank matrix (ill-conditioned eigenvalues). -- Function File: C = gallery ("gcdmat", N) Create a greatest common divisor matrix. C is an N-by-N matrix whose values correspond to the greatest common divisor of its coordinate values, i.e., C(i,j) correspond 'gcd (i, j)'. -- Function File: A = gallery ("gearmat", N) -- Function File: A = gallery ("gearmat", N, I) -- Function File: A = gallery ("gearmat", N, I, J) Create a Gear matrix. -- Function File: G = gallery ("grcar", N) -- Function File: G = gallery ("grcar", N, K) Create a Toeplitz matrix with sensitive eigenvalues. -- Function File: A = gallery ("hanowa", N) -- Function File: A = gallery ("hanowa", N, D) Create a matrix whose eigenvalues lie on a vertical line in the complex plane. -- Function File: V = gallery ("house", X) -- Function File: [V, BETA] = gallery ("house", X) Create a householder matrix. -- Function File: A = gallery ("integerdata", IMAX, [M N ...], J) -- Function File: A = gallery ("integerdata", IMAX, M, N, ..., J) -- Function File: A = gallery ("integerdata", [IMIN, IMAX], [M N ...], J) -- Function File: A = gallery ("integerdata", [IMIN, IMAX], M, N, ..., J) -- Function File: A = gallery ("integerdata", ..., "CLASS") Create a matrix with random integers in the range [1, IMAX]. If IMIN is given then the integers are in the range [IMIN, IMAX]. The second input is a matrix of dimensions describing the size of the output. The dimensions can also be input as comma-separated arguments. The input J is an integer index in the range [0, 2^32-1]. The values of the output matrix are always exactly the same (reproducibility) for a given size input and J index. The final optional argument determines the class of the resulting matrix. Possible values for CLASS: "uint8", "uint16", "uint32", "int8", "int16", int32", "single", "double". The default is "double". -- Function File: A = gallery ("invhess", X) -- Function File: A = gallery ("invhess", X, Y) Create the inverse of an upper Hessenberg matrix. -- Function File: A = gallery ("invol", N) Create an involutory matrix. -- Function File: A = gallery ("ipjfact", N) -- Function File: A = gallery ("ipjfact", N, K) Create a Hankel matrix with factorial elements. -- Function File: A = gallery ("jordbloc", N) -- Function File: A = gallery ("jordbloc", N, LAMBDA) Create a Jordan block. -- Function File: U = gallery ("kahan", N) -- Function File: U = gallery ("kahan", N, THETA) -- Function File: U = gallery ("kahan", N, THETA, PERT) Create a Kahan matrix (upper trapezoidal). -- Function File: A = gallery ("kms", N) -- Function File: A = gallery ("kms", N, RHO) Create a Kac-Murdock-Szego Toeplitz matrix. -- Function File: B = gallery ("krylov", A) -- Function File: B = gallery ("krylov", A, X) -- Function File: B = gallery ("krylov", A, X, J) Create a Krylov matrix. -- Function File: A = gallery ("lauchli", N) -- Function File: A = gallery ("lauchli", N, MU) Create a Lauchli matrix (rectangular). -- Function File: A = gallery ("lehmer", N) Create a Lehmer matrix (symmetric positive definite). -- Function File: T = gallery ("lesp", N) Create a tridiagonal matrix with real, sensitive eigenvalues. -- Function File: A = gallery ("lotkin", N) Create a Lotkin matrix. -- Function File: A = gallery ("minij", N) Create a symmetric positive definite matrix MIN(i,j). -- Function File: A = gallery ("moler", N) -- Function File: A = gallery ("moler", N, ALPHA) Create a Moler matrix (symmetric positive definite). -- Function File: [A, T] = gallery ("neumann", N) Create a singular matrix from the discrete Neumann problem (sparse). -- Function File: A = gallery ("normaldata", [M N ...], J) -- Function File: A = gallery ("normaldata", M, N, ..., J) -- Function File: A = gallery ("normaldata", ..., "CLASS") Create a matrix with random samples from the standard normal distribution (mean = 0, std = 1). The first input is a matrix of dimensions describing the size of the output. The dimensions can also be input as comma-separated arguments. The input J is an integer index in the range [0, 2^32-1]. The values of the output matrix are always exactly the same (reproducibility) for a given size input and J index. The final optional argument determines the class of the resulting matrix. Possible values for CLASS: "single", "double". The default is "double". -- Function File: Q = gallery ("orthog", N) -- Function File: Q = gallery ("orthog", N, K) Create orthogonal and nearly orthogonal matrices. -- Function File: A = gallery ("parter", N) Create a Parter matrix (a Toeplitz matrix with singular values near pi). -- Function File: P = gallery ("pei", N) -- Function File: P = gallery ("pei", N, ALPHA) Create a Pei matrix. -- Function File: A = gallery ("Poisson", N) Create a block tridiagonal matrix from Poisson's equation (sparse). -- Function File: A = gallery ("prolate", N) -- Function File: A = gallery ("prolate", N, W) Create a prolate matrix (symmetric, ill-conditioned Toeplitz matrix). -- Function File: H = gallery ("randhess", X) Create a random, orthogonal upper Hessenberg matrix. -- Function File: A = gallery ("rando", N) -- Function File: A = gallery ("rando", N, K) Create a random matrix with elements -1, 0 or 1. -- Function File: A = gallery ("randsvd", N) -- Function File: A = gallery ("randsvd", N, KAPPA) -- Function File: A = gallery ("randsvd", N, KAPPA, MODE) -- Function File: A = gallery ("randsvd", N, KAPPA, MODE, KL) -- Function File: A = gallery ("randsvd", N, KAPPA, MODE, KL, KU) Create a random matrix with pre-assigned singular values. -- Function File: A = gallery ("redheff", N) Create a zero and ones matrix of Redheffer associated with the Riemann hypothesis. -- Function File: A = gallery ("riemann", N) Create a matrix associated with the Riemann hypothesis. -- Function File: A = gallery ("ris", N) Create a symmetric Hankel matrix. -- Function File: A = gallery ("smoke", N) -- Function File: A = gallery ("smoke", N, K) Create a complex matrix, with a 'smoke ring' pseudospectrum. -- Function File: T = gallery ("toeppd", N) -- Function File: T = gallery ("toeppd", N, M) -- Function File: T = gallery ("toeppd", N, M, W) -- Function File: T = gallery ("toeppd", N, M, W, THETA) Create a symmetric positive definite Toeplitz matrix. -- Function File: P = gallery ("toeppen", N) -- Function File: P = gallery ("toeppen", N, A) -- Function File: P = gallery ("toeppen", N, A, B) -- Function File: P = gallery ("toeppen", N, A, B, C) -- Function File: P = gallery ("toeppen", N, A, B, C, D) -- Function File: P = gallery ("toeppen", N, A, B, C, D, E) Create a pentadiagonal Toeplitz matrix (sparse). -- Function File: A = gallery ("tridiag", X, Y, Z) -- Function File: A = gallery ("tridiag", N) -- Function File: A = gallery ("tridiag", N, C, D, E) Create a tridiagonal matrix (sparse). -- Function File: T = gallery ("triw", N) -- Function File: T = gallery ("triw", N, ALPHA) -- Function File: T = gallery ("triw", N, ALPHA, K) Create an upper triangular matrix discussed by Kahan, Golub, and Wilkinson. -- Function File: A = gallery ("uniformdata", [M N ...], J) -- Function File: A = gallery ("uniformdata", M, N, ..., J) -- Function File: A = gallery ("uniformdata", ..., "CLASS") Create a matrix with random samples from the standard uniform distribution (range [0,1]). The first input is a matrix of dimensions describing the size of the output. The dimensions can also be input as comma-separated arguments. The input J is an integer index in the range [0, 2^32-1]. The values of the output matrix are always exactly the same (reproducibility) for a given size input and J index. The final optional argument determines the class of the resulting matrix. Possible values for CLASS: "single", "double". The default is "double". -- Function File: A = gallery ("wathen", NX, NY) -- Function File: A = gallery ("wathen", NX, NY, K) Create the Wathen matrix. -- Function File: [A, B] = gallery ("wilk", N) Create various specific matrices devised/discussed by Wilkinson. -- Function File: hadamard (N) Construct a Hadamard matrix (Hn) of size N-by-N. The size N must be of the form 2^k * p in which p is one of 1, 12, 20 or 28. The returned matrix is normalized, meaning 'Hn(:,1) == 1' and 'Hn(1,:) == 1'. Some of the properties of Hadamard matrices are: * 'kron (Hm, Hn)' is a Hadamard matrix of size M-by-N. * 'Hn * Hn' = N * eye (N)'. * The rows of Hn are orthogonal. * 'det (A) <= abs (det (Hn))' for all A with 'abs (A(i, j)) <= 1'. * Multiplying any row or column by -1 and the matrix will remain a Hadamard matrix. See also: *note compan: XREFcompan, *note hankel: XREFhankel, *note toeplitz: XREFtoeplitz. -- Function File: hankel (C) -- Function File: hankel (C, R) Return the Hankel matrix constructed from the first column C, and (optionally) the last row R. If the last element of C is not the same as the first element of R, the last element of C is used. If the second argument is omitted, it is assumed to be a vector of zeros with the same size as C. A Hankel matrix formed from an m-vector C, and an n-vector R, has the elements H(i,j) = c(i+j-1), i+j-1 <= m; H(i,j) = r(i+j-m), otherwise See also: *note hadamard: XREFhadamard, *note toeplitz: XREFtoeplitz. -- Function File: hilb (N) Return the Hilbert matrix of order N. The i,j element of a Hilbert matrix is defined as H(i, j) = 1 / (i + j - 1) Hilbert matrices are close to being singular which make them difficult to invert with numerical routines. Comparing the condition number of a random matrix 5x5 matrix with that of a Hilbert matrix of order 5 reveals just how difficult the problem is. cond (rand (5)) => 14.392 cond (hilb (5)) => 4.7661e+05 See also: *note invhilb: XREFinvhilb. -- Function File: invhilb (N) Return the inverse of the Hilbert matrix of order N. This can be computed exactly using (i+j) /n+i-1\ /n+j-1\ /i+j-2\ 2 A(i,j) = -1 (i+j-1)( )( ) ( ) \ n-j / \ n-i / \ i-2 / = p(i) p(j) / (i+j-1) where k /k+n-1\ /n\ p(k) = -1 ( ) ( ) \ k-1 / \k/ The validity of this formula can easily be checked by expanding the binomial coefficients in both formulas as factorials. It can be derived more directly via the theory of Cauchy matrices. See J. W. Demmel, 'Applied Numerical Linear Algebra', p. 92. Compare this with the numerical calculation of 'inverse (hilb (n))', which suffers from the ill-conditioning of the Hilbert matrix, and the finite precision of your computer's floating point arithmetic. See also: *note hilb: XREFhilb. -- Function File: magic (N) Create an N-by-N magic square. A magic square is an arrangement of the integers '1:n^2' such that the row sums, column sums, and diagonal sums are all equal to the same value. Note: N must be greater than 2 for the magic square to exist. -- Function File: pascal (N) -- Function File: pascal (N, T) Return the Pascal matrix of order N if 'T = 0'. The default value of T is 0. When 'T = 1', return the pseudo-lower triangular Cholesky factor of the Pascal matrix (The sign of some columns may be negative). This matrix is its own inverse, that is 'pascal (N, 1) ^ 2 == eye (N)'. If 'T = -1', return the true Cholesky factor with strictly positive values on the diagonal. If 'T = 2', return a transposed and permuted version of 'pascal (N, 1)', which is the cube root of the identity matrix. That is, 'pascal (N, 2) ^ 3 == eye (N)'. See also: *note chol: XREFchol. -- Function File: rosser () Return the Rosser matrix. This is a difficult test case used to evaluate eigenvalue algorithms. See also: *note wilkinson: XREFwilkinson, *note eig: XREFeig. -- Function File: toeplitz (C) -- Function File: toeplitz (C, R) Return the Toeplitz matrix constructed from the first column C, and (optionally) the first row R. If the first element of R is not the same as the first element of C, the first element of C is used. If the second argument is omitted, the first row is taken to be the same as the first column. A square Toeplitz matrix has the form: c(0) r(1) r(2) ... r(n) c(1) c(0) r(1) ... r(n-1) c(2) c(1) c(0) ... r(n-2) . . . . . . . . . . . . . . . c(n) c(n-1) c(n-2) ... c(0) See also: *note hankel: XREFhankel. -- Function File: vander (C) -- Function File: vander (C, N) Return the Vandermonde matrix whose next to last column is C. If N is specified, it determines the number of columns; otherwise, N is taken to be equal to the length of C. A Vandermonde matrix has the form: c(1)^(n-1) ... c(1)^2 c(1) 1 c(2)^(n-1) ... c(2)^2 c(2) 1 . . . . . . . . . . . . . . . c(n)^(n-1) ... c(n)^2 c(n) 1 See also: *note polyfit: XREFpolyfit. -- Function File: wilkinson (N) Return the Wilkinson matrix of order N. Wilkinson matrices are symmetric and tridiagonal with pairs of nearly, but not exactly, equal eigenvalues. They are useful in testing the behavior and performance of eigenvalue solvers. See also: *note rosser: XREFrosser, *note eig: XREFeig.  File: octave.info, Node: Arithmetic, Next: Linear Algebra, Prev: Matrix Manipulation, Up: Top 17 Arithmetic ************* Unless otherwise noted, all of the functions described in this chapter will work for real and complex scalar, vector, or matrix arguments. Functions described as "mapping functions" apply the given operation individually to each element when given a matrix argument. For example: sin ([1, 2; 3, 4]) => 0.84147 0.90930 0.14112 -0.75680 * Menu: * Exponents and Logarithms:: * Complex Arithmetic:: * Trigonometry:: * Sums and Products:: * Utility Functions:: * Special Functions:: * Rational Approximations:: * Coordinate Transformations:: * Mathematical Constants::  File: octave.info, Node: Exponents and Logarithms, Next: Complex Arithmetic, Up: Arithmetic 17.1 Exponents and Logarithms ============================= -- Mapping Function: exp (X) Compute 'e^x' for each element of X. To compute the matrix exponential, see *note Linear Algebra::. See also: *note log: XREFlog. -- Mapping Function: expm1 (X) Compute 'exp (X) - 1' accurately in the neighborhood of zero. See also: *note exp: XREFexp. -- Mapping Function: log (X) Compute the natural logarithm, 'ln (X)', for each element of X. To compute the matrix logarithm, see *note Linear Algebra::. See also: *note exp: XREFexp, *note log1p: XREFlog1p, *note log2: XREFlog2, *note log10: XREFlog10, *note logspace: XREFlogspace. -- Function File: reallog (X) Return the real-valued natural logarithm of each element of X. If any element results in a complex return value 'reallog' aborts and issues an error. See also: *note log: XREFlog, *note realpow: XREFrealpow, *note realsqrt: XREFrealsqrt. -- Mapping Function: log1p (X) Compute 'log (1 + X)' accurately in the neighborhood of zero. See also: *note log: XREFlog, *note exp: XREFexp, *note expm1: XREFexpm1. -- Mapping Function: log10 (X) Compute the base-10 logarithm of each element of X. See also: *note log: XREFlog, *note log2: XREFlog2, *note logspace: XREFlogspace, *note exp: XREFexp. -- Mapping Function: log2 (X) -- Mapping Function: [F, E] = log2 (X) Compute the base-2 logarithm of each element of X. If called with two output arguments, split X into binary mantissa and exponent so that '1/2 <= abs(f) < 1' and E is an integer. If 'x = 0', 'f = e = 0'. See also: *note pow2: XREFpow2, *note log: XREFlog, *note log10: XREFlog10, *note exp: XREFexp. -- Function File: pow2 (X) -- Function File: pow2 (F, E) With one input argument, compute 2 .^ x for each element of X. With two input arguments, return f .* (2 .^ e). See also: *note log2: XREFlog2, *note nextpow2: XREFnextpow2, *note power: XREFpower. -- Function File: nextpow2 (X) Compute the exponent for the smallest power of two larger than the input. For each element in the input array X, return the first integer N such that 2^n >= abs (x). See also: *note pow2: XREFpow2, *note log2: XREFlog2. -- Function File: realpow (X, Y) Compute the real-valued, element-by-element power operator. This is equivalent to 'X .^ Y', except that 'realpow' reports an error if any return value is complex. See also: *note power: XREFpower, *note reallog: XREFreallog, *note realsqrt: XREFrealsqrt. -- Mapping Function: sqrt (X) Compute the square root of each element of X. If X is negative, a complex result is returned. To compute the matrix square root, see *note Linear Algebra::. See also: *note realsqrt: XREFrealsqrt, *note nthroot: XREFnthroot. -- Function File: realsqrt (X) Return the real-valued square root of each element of X. If any element results in a complex return value 'realsqrt' aborts and issues an error. See also: *note sqrt: XREFsqrt, *note realpow: XREFrealpow, *note reallog: XREFreallog. -- Mapping Function: cbrt (X) Compute the real cube root of each element of X. Unlike 'X^(1/3)', the result will be negative if X is negative. See also: *note nthroot: XREFnthroot. -- Function File: nthroot (X, N) Compute the real (non-complex) N-th root of X. X must have all real entries and N must be a scalar. If N is an even integer and X has negative entries then 'nthroot' aborts and issues an error. Example: nthroot (-1, 3) => -1 (-1) ^ (1 / 3) => 0.50000 - 0.86603i See also: *note realsqrt: XREFrealsqrt, *note sqrt: XREFsqrt, *note cbrt: XREFcbrt.  File: octave.info, Node: Complex Arithmetic, Next: Trigonometry, Prev: Exponents and Logarithms, Up: Arithmetic 17.2 Complex Arithmetic ======================= In the descriptions of the following functions, Z is the complex number X + IY, where I is defined as 'sqrt (-1)'. -- Mapping Function: abs (Z) Compute the magnitude of Z. The magnitude is defined as |Z| = 'sqrt (x^2 + y^2)'. For example: abs (3 + 4i) => 5 See also: *note arg: XREFarg. -- Mapping Function: arg (Z) -- Mapping Function: angle (Z) Compute the argument, i.e., angle of Z. This is defined as, THETA = 'atan2 (Y, X)', in radians. For example: arg (3 + 4i) => 0.92730 See also: *note abs: XREFabs. -- Mapping Function: conj (Z) Return the complex conjugate of Z. The complex conjugate is defined as 'conj (Z)' = X - IY. See also: *note real: XREFreal, *note imag: XREFimag. -- Function File: cplxpair (Z) -- Function File: cplxpair (Z, TOL) -- Function File: cplxpair (Z, TOL, DIM) Sort the numbers Z into complex conjugate pairs ordered by increasing real part. The negative imaginary complex numbers are placed first within each pair. All real numbers (those with 'abs (imag (Z) / Z) < TOL') are placed after the complex pairs. If TOL is unspecified the default value is 100*'eps'. By default the complex pairs are sorted along the first non-singleton dimension of Z. If DIM is specified, then the complex pairs are sorted along this dimension. Signal an error if some complex numbers could not be paired. Signal an error if all complex numbers are not exact conjugates (to within TOL). Note that there is no defined order for pairs with identical real parts but differing imaginary parts. cplxpair (exp(2i*pi*[0:4]'/5)) == exp(2i*pi*[3; 2; 4; 1; 0]/5) -- Mapping Function: imag (Z) Return the imaginary part of Z as a real number. See also: *note real: XREFreal, *note conj: XREFconj. -- Mapping Function: real (Z) Return the real part of Z. See also: *note imag: XREFimag, *note conj: XREFconj.  File: octave.info, Node: Trigonometry, Next: Sums and Products, Prev: Complex Arithmetic, Up: Arithmetic 17.3 Trigonometry ================= Octave provides the following trigonometric functions where angles are specified in radians. To convert from degrees to radians multiply by 'pi/180' (e.g., 'sin (30 * pi/180)' returns the sine of 30 degrees). As an alternative, Octave provides a number of trigonometric functions which work directly on an argument specified in degrees. These functions are named after the base trigonometric function with a 'd' suffix. For example, 'sin' expects an angle in radians while 'sind' expects an angle in degrees. Octave uses the C library trigonometric functions. It is expected that these functions are defined by the ISO/IEC 9899 Standard. This Standard is available at: . Section F.9.1 deals with the trigonometric functions. The behavior of most of the functions is relatively straightforward. However, there are some exceptions to the standard behavior. Many of the exceptions involve the behavior for -0. The most complex case is atan2. Octave exactly implements the behavior given in the Standard. Including 'atan2(+- 0, 0)' returns '+- pi'. It should be noted that MATLAB uses different definitions which apparently do not distinguish -0. -- Mapping Function: sin (X) Compute the sine for each element of X in radians. See also: *note asin: XREFasin, *note sind: XREFsind, *note sinh: XREFsinh. -- Mapping Function: cos (X) Compute the cosine for each element of X in radians. See also: *note acos: XREFacos, *note cosd: XREFcosd, *note cosh: XREFcosh. -- Mapping Function: tan (Z) Compute the tangent for each element of X in radians. See also: *note atan: XREFatan, *note tand: XREFtand, *note tanh: XREFtanh. -- Mapping Function: sec (X) Compute the secant for each element of X in radians. See also: *note asec: XREFasec, *note secd: XREFsecd, *note sech: XREFsech. -- Mapping Function: csc (X) Compute the cosecant for each element of X in radians. See also: *note acsc: XREFacsc, *note cscd: XREFcscd, *note csch: XREFcsch. -- Mapping Function: cot (X) Compute the cotangent for each element of X in radians. See also: *note acot: XREFacot, *note cotd: XREFcotd, *note coth: XREFcoth. -- Mapping Function: asin (X) Compute the inverse sine in radians for each element of X. See also: *note sin: XREFsin, *note asind: XREFasind. -- Mapping Function: acos (X) Compute the inverse cosine in radians for each element of X. See also: *note cos: XREFcos, *note acosd: XREFacosd. -- Mapping Function: atan (X) Compute the inverse tangent in radians for each element of X. See also: *note tan: XREFtan, *note atand: XREFatand. -- Mapping Function: asec (X) Compute the inverse secant in radians for each element of X. See also: *note sec: XREFsec, *note asecd: XREFasecd. -- Mapping Function: acsc (X) Compute the inverse cosecant in radians for each element of X. See also: *note csc: XREFcsc, *note acscd: XREFacscd. -- Mapping Function: acot (X) Compute the inverse cotangent in radians for each element of X. See also: *note cot: XREFcot, *note acotd: XREFacotd. -- Mapping Function: sinh (X) Compute the hyperbolic sine for each element of X. See also: *note asinh: XREFasinh, *note cosh: XREFcosh, *note tanh: XREFtanh. -- Mapping Function: cosh (X) Compute the hyperbolic cosine for each element of X. See also: *note acosh: XREFacosh, *note sinh: XREFsinh, *note tanh: XREFtanh. -- Mapping Function: tanh (X) Compute hyperbolic tangent for each element of X. See also: *note atanh: XREFatanh, *note sinh: XREFsinh, *note cosh: XREFcosh. -- Mapping Function: sech (X) Compute the hyperbolic secant of each element of X. See also: *note asech: XREFasech. -- Mapping Function: csch (X) Compute the hyperbolic cosecant of each element of X. See also: *note acsch: XREFacsch. -- Mapping Function: coth (X) Compute the hyperbolic cotangent of each element of X. See also: *note acoth: XREFacoth. -- Mapping Function: asinh (X) Compute the inverse hyperbolic sine for each element of X. See also: *note sinh: XREFsinh. -- Mapping Function: acosh (X) Compute the inverse hyperbolic cosine for each element of X. See also: *note cosh: XREFcosh. -- Mapping Function: atanh (X) Compute the inverse hyperbolic tangent for each element of X. See also: *note tanh: XREFtanh. -- Mapping Function: asech (X) Compute the inverse hyperbolic secant of each element of X. See also: *note sech: XREFsech. -- Mapping Function: acsch (X) Compute the inverse hyperbolic cosecant of each element of X. See also: *note csch: XREFcsch. -- Mapping Function: acoth (X) Compute the inverse hyperbolic cotangent of each element of X. See also: *note coth: XREFcoth. -- Mapping Function: atan2 (Y, X) Compute atan (Y / X) for corresponding elements of Y and X. Y and X must match in size and orientation. See also: *note tan: XREFtan, *note tand: XREFtand, *note tanh: XREFtanh, *note atanh: XREFatanh. Octave provides the following trigonometric functions where angles are specified in degrees. These functions produce true zeros at the appropriate intervals rather than the small round-off error that occurs when using radians. For example: cosd (90) => 0 cos (pi/2) => 6.1230e-17 -- Function File: sind (X) Compute the sine for each element of X in degrees. Returns zero for elements where 'X/180' is an integer. See also: *note asind: XREFasind, *note sin: XREFsin. -- Function File: cosd (X) Compute the cosine for each element of X in degrees. Returns zero for elements where '(X-90)/180' is an integer. See also: *note acosd: XREFacosd, *note cos: XREFcos. -- Function File: tand (X) Compute the tangent for each element of X in degrees. Returns zero for elements where 'X/180' is an integer and 'Inf' for elements where '(X-90)/180' is an integer. See also: *note atand: XREFatand, *note tan: XREFtan. -- Function File: secd (X) Compute the secant for each element of X in degrees. See also: *note asecd: XREFasecd, *note sec: XREFsec. -- Function File: cscd (X) Compute the cosecant for each element of X in degrees. See also: *note acscd: XREFacscd, *note csc: XREFcsc. -- Function File: cotd (X) Compute the cotangent for each element of X in degrees. See also: *note acotd: XREFacotd, *note cot: XREFcot. -- Function File: asind (X) Compute the inverse sine in degrees for each element of X. See also: *note sind: XREFsind, *note asin: XREFasin. -- Function File: acosd (X) Compute the inverse cosine in degrees for each element of X. See also: *note cosd: XREFcosd, *note acos: XREFacos. -- Function File: atand (X) Compute the inverse tangent in degrees for each element of X. See also: *note tand: XREFtand, *note atan: XREFatan. -- Function File: atan2d (Y, X) Compute atan2 (Y / X) in degrees for corresponding elements from Y and X. See also: *note tand: XREFtand, *note atan2: XREFatan2. -- Function File: asecd (X) Compute the inverse secant in degrees for each element of X. See also: *note secd: XREFsecd, *note asec: XREFasec. -- Function File: acscd (X) Compute the inverse cosecant in degrees for each element of X. See also: *note cscd: XREFcscd, *note acsc: XREFacsc. -- Function File: acotd (X) Compute the inverse cotangent in degrees for each element of X. See also: *note cotd: XREFcotd, *note acot: XREFacot.  File: octave.info, Node: Sums and Products, Next: Utility Functions, Prev: Trigonometry, Up: Arithmetic 17.4 Sums and Products ====================== -- Built-in Function: sum (X) -- Built-in Function: sum (X, DIM) -- Built-in Function: sum (..., "native") -- Built-in Function: sum (..., "double") -- Built-in Function: sum (..., "extra") Sum of elements along dimension DIM. If DIM is omitted, it defaults to the first non-singleton dimension. The optional "type" input determines the class of the variable used for calculations. If the argument "native" is given, then the operation is performed in the same type as the original argument, rather than the default double type. For example: sum ([true, true]) => 2 sum ([true, true], "native") => true On the contrary, if "double" is given, the sum is performed in double precision even for single precision inputs. For double precision inputs, the "extra" option will use a more accurate algorithm than straightforward summation. For single precision inputs, "extra" is the same as "double". Otherwise, "extra" has no effect. See also: *note cumsum: XREFcumsum, *note sumsq: XREFsumsq, *note prod: XREFprod. -- Built-in Function: prod (X) -- Built-in Function: prod (X, DIM) -- Built-in Function: prod (..., "native") -- Built-in Function: prod (..., "double") Product of elements along dimension DIM. If DIM is omitted, it defaults to the first non-singleton dimension. The optional "type" input determines the class of the variable used for calculations. If the argument "native" is given, then the operation is performed in the same type as the original argument, rather than the default double type. For example: prod ([true, true]) => 1 prod ([true, true], "native") => true On the contrary, if "double" is given, the operation is performed in double precision even for single precision inputs. See also: *note cumprod: XREFcumprod, *note sum: XREFsum. -- Built-in Function: cumsum (X) -- Built-in Function: cumsum (X, DIM) -- Built-in Function: cumsum (..., "native") -- Built-in Function: cumsum (..., "double") -- Built-in Function: cumsum (..., "extra") Cumulative sum of elements along dimension DIM. If DIM is omitted, it defaults to the first non-singleton dimension. See 'sum' for an explanation of the optional parameters "native", "double", and "extra". See also: *note sum: XREFsum, *note cumprod: XREFcumprod. -- Built-in Function: cumprod (X) -- Built-in Function: cumprod (X, DIM) Cumulative product of elements along dimension DIM. If DIM is omitted, it defaults to the first non-singleton dimension. See also: *note prod: XREFprod, *note cumsum: XREFcumsum. -- Built-in Function: sumsq (X) -- Built-in Function: sumsq (X, DIM) Sum of squares of elements along dimension DIM. If DIM is omitted, it defaults to the first non-singleton dimension. This function is conceptually equivalent to computing sum (x .* conj (x), dim) but it uses less memory and avoids calling 'conj' if X is real. See also: *note sum: XREFsum, *note prod: XREFprod.  File: octave.info, Node: Utility Functions, Next: Special Functions, Prev: Sums and Products, Up: Arithmetic 17.5 Utility Functions ====================== -- Mapping Function: ceil (X) Return the smallest integer not less than X. This is equivalent to rounding towards positive infinity. If X is complex, return 'ceil (real (X)) + ceil (imag (X)) * I'. ceil ([-2.7, 2.7]) => -2 3 See also: *note floor: XREFfloor, *note round: XREFround, *note fix: XREFfix. -- Mapping Function: fix (X) Truncate fractional portion of X and return the integer portion. This is equivalent to rounding towards zero. If X is complex, return 'fix (real (X)) + fix (imag (X)) * I'. fix ([-2.7, 2.7]) => -2 2 See also: *note ceil: XREFceil, *note floor: XREFfloor, *note round: XREFround. -- Mapping Function: floor (X) Return the largest integer not greater than X. This is equivalent to rounding towards negative infinity. If X is complex, return 'floor (real (X)) + floor (imag (X)) * I'. floor ([-2.7, 2.7]) => -3 2 See also: *note ceil: XREFceil, *note round: XREFround, *note fix: XREFfix. -- Mapping Function: round (X) Return the integer nearest to X. If X is complex, return 'round (real (X)) + round (imag (X)) * I'. If there are two nearest integers, return the one further away from zero. round ([-2.7, 2.7]) => -3 3 See also: *note ceil: XREFceil, *note floor: XREFfloor, *note fix: XREFfix, *note roundb: XREFroundb. -- Mapping Function: roundb (X) Return the integer nearest to X. If there are two nearest integers, return the even one (banker's rounding). If X is complex, return 'roundb (real (X)) + roundb (imag (X)) * I'. See also: *note round: XREFround. -- Built-in Function: max (X) -- Built-in Function: max (X, [], DIM) -- Built-in Function: [W, IW] = max (X) -- Built-in Function: max (X, Y) Find maximum values in the array X. For a vector argument, return the maximum value. For a matrix argument, return a row vector with the maximum value of each column. For a multi-dimensional array, 'max' operates along the first non-singleton dimension. If the optional third argument DIM is present then operate along this dimension. In this case the second argument is ignored and should be set to the empty matrix. For two matrices (or a matrix and a scalar), return the pairwise maximum. Thus, max (max (X)) returns the largest element of the 2-D matrix X, and max (2:5, pi) => 3.1416 3.1416 4.0000 5.0000 compares each element of the range '2:5' with 'pi', and returns a row vector of the maximum values. For complex arguments, the magnitude of the elements are used for comparison. If the magnitudes are identical, then the results are ordered by phase angle in the range (-pi, pi]. Hence, max ([-1 i 1 -i]) => -1 because all entries have magnitude 1, but -1 has the largest phase angle with value pi. If called with one input and two output arguments, 'max' also returns the first index of the maximum value(s). Thus, [x, ix] = max ([1, 3, 5, 2, 5]) => x = 5 ix = 3 See also: *note min: XREFmin, *note cummax: XREFcummax, *note cummin: XREFcummin. -- Built-in Function: min (X) -- Built-in Function: min (X, [], DIM) -- Built-in Function: [W, IW] = min (X) -- Built-in Function: min (X, Y) Find minimum values in the array X. For a vector argument, return the minimum value. For a matrix argument, return a row vector with the minimum value of each column. For a multi-dimensional array, 'min' operates along the first non-singleton dimension. If the optional third argument DIM is present then operate along this dimension. In this case the second argument is ignored and should be set to the empty matrix. For two matrices (or a matrix and a scalar), return the pairwise minimum. Thus, min (min (X)) returns the smallest element of the 2-D matrix X, and min (2:5, pi) => 2.0000 3.0000 3.1416 3.1416 compares each element of the range '2:5' with 'pi', and returns a row vector of the minimum values. For complex arguments, the magnitude of the elements are used for comparison. If the magnitudes are identical, then the results are ordered by phase angle in the range (-pi, pi]. Hence, min ([-1 i 1 -i]) => -i because all entries have magnitude 1, but -i has the smallest phase angle with value -pi/2. If called with one input and two output arguments, 'min' also returns the first index of the minimum value(s). Thus, [x, ix] = min ([1, 3, 0, 2, 0]) => x = 0 ix = 3 See also: *note max: XREFmax, *note cummin: XREFcummin, *note cummax: XREFcummax. -- Built-in Function: cummax (X) -- Built-in Function: cummax (X, DIM) -- Built-in Function: [W, IW] = cummax (...) Return the cumulative maximum values along dimension DIM. If DIM is unspecified it defaults to column-wise operation. For example: cummax ([1 3 2 6 4 5]) => 1 3 3 6 6 6 If called with two output arguments the index of the maximum value is also returned. [w, iw] = cummax ([1 3 2 6 4 5]) => w = 1 3 3 6 6 6 iw = 1 2 2 4 4 4 See also: *note cummin: XREFcummin, *note max: XREFmax, *note min: XREFmin. -- Built-in Function: cummin (X) -- Built-in Function: cummin (X, DIM) -- Built-in Function: [W, IW] = cummin (X) Return the cumulative minimum values along dimension DIM. If DIM is unspecified it defaults to column-wise operation. For example: cummin ([5 4 6 2 3 1]) => 5 4 4 2 2 1 If called with two output arguments the index of the minimum value is also returned. [w, iw] = cummin ([5 4 6 2 3 1]) => w = 5 4 4 2 2 1 iw = 1 2 2 4 4 6 See also: *note cummax: XREFcummax, *note min: XREFmin, *note max: XREFmax. -- Built-in Function: hypot (X, Y) -- Built-in Function: hypot (X, Y, Z, ...) Compute the element-by-element square root of the sum of the squares of X and Y. This is equivalent to 'sqrt (X.^2 + Y.^2)', but is calculated in a manner that avoids overflows for large values of X or Y. 'hypot' can also be called with more than 2 arguments; in this case, the arguments are accumulated from left to right: hypot (hypot (X, Y), Z) hypot (hypot (hypot (X, Y), Z), W), etc. -- Function File: DX = gradient (M) -- Function File: [DX, DY, DZ, ...] = gradient (M) -- Function File: [...] = gradient (M, S) -- Function File: [...] = gradient (M, X, Y, Z, ...) -- Function File: [...] = gradient (F, X0) -- Function File: [...] = gradient (F, X0, S) -- Function File: [...] = gradient (F, X0, X, Y, ...) Calculate the gradient of sampled data or a function. If M is a vector, calculate the one-dimensional gradient of M. If M is a matrix the gradient is calculated for each dimension. '[DX, DY] = gradient (M)' calculates the one-dimensional gradient for X and Y direction if M is a matrix. Additional return arguments can be use for multi-dimensional matrices. A constant spacing between two points can be provided by the S parameter. If S is a scalar, it is assumed to be the spacing for all dimensions. Otherwise, separate values of the spacing can be supplied by the X, ... arguments. Scalar values specify an equidistant spacing. Vector values for the X, ... arguments specify the coordinate for that dimension. The length must match their respective dimension of M. At boundary points a linear extrapolation is applied. Interior points are calculated with the first approximation of the numerical gradient y'(i) = 1/(x(i+1)-x(i-1)) * (y(i-1)-y(i+1)). If the first argument F is a function handle, the gradient of the function at the points in X0 is approximated using central difference. For example, 'gradient (@cos, 0)' approximates the gradient of the cosine function in the point x0 = 0. As with sampled data, the spacing values between the points from which the gradient is estimated can be set via the S or DX, DY, ... arguments. By default a spacing of 1 is used. See also: *note diff: XREFdiff, *note del2: XREFdel2. -- Built-in Function: dot (X, Y, DIM) Compute the dot product of two vectors. If X and Y are matrices, calculate the dot products along the first non-singleton dimension. If the optional argument DIM is given, calculate the dot products along this dimension. This is equivalent to 'sum (conj (X) .* Y, DIM)', but avoids forming a temporary array and is faster. When X and Y are column vectors, the result is equivalent to 'X' * Y'. See also: *note cross: XREFcross, *note divergence: XREFdivergence. -- Function File: cross (X, Y) -- Function File: cross (X, Y, DIM) Compute the vector cross product of two 3-dimensional vectors X and Y. If X and Y are matrices, the cross product is applied along the first dimension with three elements. The optional argument DIM forces the cross product to be calculated along the specified dimension. Example Code: cross ([1,1,0], [0,1,1]) => [ 1; -1; 1 ] See also: *note dot: XREFdot, *note curl: XREFcurl, *note divergence: XREFdivergence. -- Function File: DIV = divergence (X, Y, Z, FX, FY, FZ) -- Function File: DIV = divergence (FX, FY, FZ) -- Function File: DIV = divergence (X, Y, FX, FY) -- Function File: DIV = divergence (FX, FY) Calculate divergence of a vector field given by the arrays FX, FY, and FZ or FX, FY respectively. d d d div F(x,y,z) = -- F(x,y,z) + -- F(x,y,z) + -- F(x,y,z) dx dy dz The coordinates of the vector field can be given by the arguments X, Y, Z or X, Y respectively. See also: *note curl: XREFcurl, *note gradient: XREFgradient, *note del2: XREFdel2, *note dot: XREFdot. -- Function File: [CX, CY, CZ, V] = curl (X, Y, Z, FX, FY, FZ) -- Function File: [CZ, V] = curl (X, Y, FX, FY) -- Function File: [...] = curl (FX, FY, FZ) -- Function File: [...] = curl (FX, FY) -- Function File: V = curl (...) Calculate curl of vector field given by the arrays FX, FY, and FZ or FX, FY respectively. / d d d d d d \ curl F(x,y,z) = | -- Fz - -- Fy, -- Fx - -- Fz, -- Fy - -- Fx | \ dy dz dz dx dx dy / The coordinates of the vector field can be given by the arguments X, Y, Z or X, Y respectively. V calculates the scalar component of the angular velocity vector in direction of the z-axis for two-dimensional input. For three-dimensional input the scalar rotation is calculated at each grid point in direction of the vector field at that point. See also: *note divergence: XREFdivergence, *note gradient: XREFgradient, *note del2: XREFdel2, *note cross: XREFcross. -- Function File: D = del2 (M) -- Function File: D = del2 (M, H) -- Function File: D = del2 (M, DX, DY, ...) Calculate the discrete Laplace operator. For a 2-dimensional matrix M this is defined as 1 / d^2 d^2 \ D = --- * | --- M(x,y) + --- M(x,y) | 4 \ dx^2 dy^2 / For N-dimensional arrays the sum in parentheses is expanded to include second derivatives over the additional higher dimensions. The spacing between evaluation points may be defined by H, which is a scalar defining the equidistant spacing in all dimensions. Alternatively, the spacing in each dimension may be defined separately by DX, DY, etc. A scalar spacing argument defines equidistant spacing, whereas a vector argument can be used to specify variable spacing. The length of the spacing vectors must match the respective dimension of M. The default spacing value is 1. At least 3 data points are needed for each dimension. Boundary points are calculated from the linear extrapolation of interior points. See also: *note gradient: XREFgradient, *note diff: XREFdiff. -- Function File: factorial (N) Return the factorial of N where N is a real non-negative integer. If N is a scalar, this is equivalent to 'prod (1:N)'. For vector or matrix arguments, return the factorial of each element in the array. For non-integers see the generalized factorial function 'gamma'. Note that the factorial function grows large quite quickly, and even with double precision values overflow will occur if N > 171. For such cases consider 'gammaln'. See also: *note prod: XREFprod, *note gamma: XREFgamma, *note gammaln: XREFgammaln. -- Function File: PF = factor (Q) -- Function File: [PF, N] = factor (Q) Return the prime factorization of Q. The prime factorization is defined as 'prod (PF) == Q' where every element of PF is a prime number. If 'Q == 1', return 1. With two output arguments, return the unique prime factors PF and their multiplicities. That is, 'prod (PF .^ N) == Q'. Implementation Note: The input Q must be less than 'bitmax' (9.0072e+15) in order to factor correctly. See also: *note gcd: XREFgcd, *note lcm: XREFlcm, *note isprime: XREFisprime, *note primes: XREFprimes. -- Built-in Function: G = gcd (A1, A2, ...) -- Built-in Function: [G, V1, ...] = gcd (A1, A2, ...) Compute the greatest common divisor of A1, A2, .... If more than one argument is given then all arguments must be the same size or scalar. In this case the greatest common divisor is calculated for each element individually. All elements must be ordinary or Gaussian (complex) integers. Note that for Gaussian integers, the gcd is only unique up to a phase factor (multiplication by 1, -1, i, or -i), so an arbitrary greatest common divisor among the four possible is returned. Optional return arguments V1, ..., contain integer vectors such that, G = V1 .* A1 + V2 .* A2 + ... Example code: gcd ([15, 9], [20, 18]) => 5 9 See also: *note lcm: XREFlcm, *note factor: XREFfactor, *note isprime: XREFisprime. -- Mapping Function: lcm (X, Y) -- Mapping Function: lcm (X, Y, ...) Compute the least common multiple of X and Y, or of the list of all arguments. All elements must be numeric and of the same size or scalar. See also: *note factor: XREFfactor, *note gcd: XREFgcd, *note isprime: XREFisprime. -- Function File: chop (X, NDIGITS, BASE) Truncate elements of X to a length of NDIGITS such that the resulting numbers are exactly divisible by BASE. If BASE is not specified it defaults to 10. chop (-pi, 5, 10) => -3.14200000000000 chop (-pi, 5, 5) => -3.14150000000000 -- Mapping Function: rem (X, Y) Return the remainder of the division 'X / Y'. The remainder is computed using the expression x - y .* fix (x ./ y) An error message is printed if the dimensions of the arguments do not agree, or if either of the arguments is complex. See also: *note mod: XREFmod. -- Mapping Function: mod (X, Y) Compute the modulo of X and Y. Conceptually this is given by x - y .* floor (x ./ y) and is written such that the correct modulus is returned for integer types. This function handles negative values correctly. That is, 'mod (-1, 3)' is 2, not -1, as 'rem (-1, 3)' returns. 'mod (X, 0)' returns X. An error results if the dimensions of the arguments do not agree, or if either of the arguments is complex. See also: *note rem: XREFrem. -- Function File: primes (N) Return all primes up to N. The output data class (double, single, uint32, etc.) is the same as the input class of N. The algorithm used is the Sieve of Eratosthenes. Notes: If you need a specific number of primes you can use the fact that the distance from one prime to the next is, on average, proportional to the logarithm of the prime. Integrating, one finds that there are about k primes less than k*log (5*k). See also 'list_primes' if you need a specific number N of primes. See also: *note list_primes: XREFlist_primes, *note isprime: XREFisprime. -- Function File: list_primes () -- Function File: list_primes (N) List the first N primes. If N is unspecified, the first 25 primes are listed. See also: *note primes: XREFprimes, *note isprime: XREFisprime. -- Mapping Function: sign (X) Compute the "signum" function. This is defined as -1, x < 0; sign (x) = 0, x = 0; 1, x > 0. For complex arguments, 'sign' returns 'x ./ abs (X)'. Note that 'sign (-0.0)' is 0. Although IEEE 754 floating point allows zero to be signed, 0.0 and -0.0 compare equal. If you must test whether zero is signed, use the 'signbit' function. See also: *note signbit: XREFsignbit. -- Mapping Function: signbit (X) Return logical true if the value of X has its sign bit set and false otherwise. This behavior is consistent with the other logical functions. See *note Logical Values::. The behavior differs from the C language function which returns nonzero if the sign bit is set. This is not the same as 'x < 0.0', because IEEE 754 floating point allows zero to be signed. The comparison '-0.0 < 0.0' is false, but 'signbit (-0.0)' will return a nonzero value. See also: *note sign: XREFsign.  File: octave.info, Node: Special Functions, Next: Rational Approximations, Prev: Utility Functions, Up: Arithmetic 17.6 Special Functions ====================== -- Built-in Function: [A, IERR] = airy (K, Z, OPT) Compute Airy functions of the first and second kind, and their derivatives. K Function Scale factor (if "opt" is supplied) --- -------- --------------------------------------- 0 Ai (Z) exp ((2/3) * Z * sqrt (Z)) 1 dAi(Z)/dZ exp ((2/3) * Z * sqrt (Z)) 2 Bi (Z) exp (-abs (real ((2/3) * Z * sqrt (Z)))) 3 dBi(Z)/dZ exp (-abs (real ((2/3) * Z * sqrt (Z)))) The function call 'airy (Z)' is equivalent to 'airy (0, Z)'. The result is the same size as Z. If requested, IERR contains the following status information and is the same size as the result. 0. Normal return. 1. Input error, return 'NaN'. 2. Overflow, return 'Inf'. 3. Loss of significance by argument reduction results in less than half of machine accuracy. 4. Complete loss of significance by argument reduction, return 'NaN'. 5. Error--no computation, algorithm termination condition not met, return 'NaN'. -- Built-in Function: [J, IERR] = besselj (ALPHA, X, OPT) -- Built-in Function: [Y, IERR] = bessely (ALPHA, X, OPT) -- Built-in Function: [I, IERR] = besseli (ALPHA, X, OPT) -- Built-in Function: [K, IERR] = besselk (ALPHA, X, OPT) -- Built-in Function: [H, IERR] = besselh (ALPHA, K, X, OPT) Compute Bessel or Hankel functions of various kinds: 'besselj' Bessel functions of the first kind. If the argument OPT is 1 or true, the result is multiplied by 'exp (-abs (imag (X)))'. 'bessely' Bessel functions of the second kind. If the argument OPT is 1 or true, the result is multiplied by 'exp (-abs (imag (X)))'. 'besseli' Modified Bessel functions of the first kind. If the argument OPT is 1 or true, the result is multiplied by 'exp (-abs (real (X)))'. 'besselk' Modified Bessel functions of the second kind. If the argument OPT is 1 or true, the result is multiplied by 'exp (X)'. 'besselh' Compute Hankel functions of the first (K = 1) or second (K = 2) kind. If the argument OPT is 1 or true, the result is multiplied by 'exp (-I*X)' for K = 1 or 'exp (I*X)' for K = 2. If ALPHA is a scalar, the result is the same size as X. If X is a scalar, the result is the same size as ALPHA. If ALPHA is a row vector and X is a column vector, the result is a matrix with 'length (X)' rows and 'length (ALPHA)' columns. Otherwise, ALPHA and X must conform and the result will be the same size. The value of ALPHA must be real. The value of X may be complex. If requested, IERR contains the following status information and is the same size as the result. 0. Normal return. 1. Input error, return 'NaN'. 2. Overflow, return 'Inf'. 3. Loss of significance by argument reduction results in less than half of machine accuracy. 4. Complete loss of significance by argument reduction, return 'NaN'. 5. Error--no computation, algorithm termination condition not met, return 'NaN'. -- Mapping Function: beta (A, B) Compute the Beta function for real inputs A and B. The Beta function definition is beta (a, b) = gamma (a) * gamma (b) / gamma (a + b). The Beta function can grow quite large and it is often more useful to work with the logarithm of the output rather than the function directly. *Note betaln: XREFbetaln, for computing the logarithm of the Beta function in an efficient manner. See also: *note betaln: XREFbetaln, *note betainc: XREFbetainc, *note betaincinv: XREFbetaincinv. -- Mapping Function: betainc (X, A, B) Compute the regularized incomplete Beta function. The regularized incomplete Beta function is defined by x 1 / betainc (x, a, b) = ----------- | t^(a-1) (1-t)^(b-1) dt. beta (a, b) / t=0 If X has more than one component, both A and B must be scalars. If X is a scalar, A and B must be of compatible dimensions. See also: *note betaincinv: XREFbetaincinv, *note beta: XREFbeta, *note betaln: XREFbetaln. -- Mapping Function: betaincinv (Y, A, B) Compute the inverse of the incomplete Beta function. The inverse is the value X such that Y == betainc (X, A, B) See also: *note betainc: XREFbetainc, *note beta: XREFbeta, *note betaln: XREFbetaln. -- Mapping Function: betaln (A, B) Compute the natural logarithm of the Beta function for real inputs A and B. 'betaln' is defined as betaln (a, b) = log (beta (a, b)) and is calculated in a way to reduce the occurrence of underflow. The Beta function can grow quite large and it is often more useful to work with the logarithm of the output rather than the function directly. See also: *note beta: XREFbeta, *note betainc: XREFbetainc, *note betaincinv: XREFbetaincinv, *note gammaln: XREFgammaln. -- Mapping Function: bincoeff (N, K) Return the binomial coefficient of N and K, defined as / \ | n | n (n-1) (n-2) ... (n-k+1) | | = ------------------------- | k | k! \ / For example: bincoeff (5, 2) => 10 In most cases, the 'nchoosek' function is faster for small scalar integer arguments. It also warns about loss of precision for big arguments. See also: *note nchoosek: XREFnchoosek. -- Function File: commutation_matrix (M, N) Return the commutation matrix K(m,n) which is the unique M*N by M*N matrix such that K(m,n) * vec(A) = vec(A') for all m by n matrices A. If only one argument M is given, K(m,m) is returned. See Magnus and Neudecker (1988), 'Matrix Differential Calculus with Applications in Statistics and Econometrics.' -- Function File: duplication_matrix (N) Return the duplication matrix Dn which is the unique n^2 by n*(n+1)/2 matrix such that Dn vech (A) = vec (A) for all symmetric n by n matrices A. See Magnus and Neudecker (1988), 'Matrix Differential Calculus with Applications in Statistics and Econometrics.' -- Mapping Function: dawson (Z) Compute the Dawson (scaled imaginary error) function. The Dawson function is defined as (sqrt (pi) / 2) * exp (-z^2) * erfi (z) See also: *note erfc: XREFerfc, *note erf: XREFerf, *note erfcx: XREFerfcx, *note erfi: XREFerfi, *note erfinv: XREFerfinv, *note erfcinv: XREFerfcinv. -- Built-in Function: [SN, CN, DN, ERR] = ellipj (U, M) -- Built-in Function: [SN, CN, DN, ERR] = ellipj (U, M, TOL) Compute the Jacobi elliptic functions SN, CN, and DN of complex argument U and real parameter M. If M is a scalar, the results are the same size as U. If U is a scalar, the results are the same size as M. If U is a column vector and M is a row vector, the results are matrices with 'length (U)' rows and 'length (M)' columns. Otherwise, U and M must conform in size and the results will be the same size as the inputs. The value of U may be complex. The value of M must be 0 <= M <= 1. The optional input TOL is currently ignored (MATLAB uses this to allow faster, less accurate approximation). If requested, ERR contains the following status information and is the same size as the result. 0. Normal return. 1. Error--no computation, algorithm termination condition not met, return 'NaN'. Reference: Milton Abramowitz and Irene A Stegun, 'Handbook of Mathematical Functions', Chapter 16 (Sections 16.4, 16.13, and 16.15), Dover, 1965. See also: *note ellipke: XREFellipke. -- Function File: K = ellipke (M) -- Function File: K = ellipke (M, TOL) -- Function File: [K, E] = ellipke (...) Compute complete elliptic integrals of the first K(M) and second E(M) kind. M must be a scalar or real array with -Inf <= M <= 1. The optional input TOL controls the stopping tolerance of the algorithm and defaults to 'eps (class (M))'. The tolerance can be increased to compute a faster, less accurate approximation. When called with one output only elliptic integrals of the first kind are returned. Mathematical Note: Elliptic integrals of the first kind are defined as 1 / dt K (m) = | ------------------------------ / sqrt ((1 - t^2)*(1 - m^2*t^2)) 0 Elliptic integrals of the second kind are defined as 1 / sqrt (1 - m^2*t^2) E (m) = | ------------------ dt / sqrt (1 - t^2) 0 Reference: Milton Abramowitz and Irene A. Stegun, 'Handbook of Mathematical Functions', Chapter 17, Dover, 1965. See also: *note ellipj: XREFellipj. -- Mapping Function: erf (Z) Compute the error function. The error function is defined as z 2 / erf (z) = --------- * | e^(-t^2) dt sqrt (pi) / t=0 See also: *note erfc: XREFerfc, *note erfcx: XREFerfcx, *note erfi: XREFerfi, *note dawson: XREFdawson, *note erfinv: XREFerfinv, *note erfcinv: XREFerfcinv. -- Mapping Function: erfc (Z) Compute the complementary error function. The complementary error function is defined as '1 - erf (Z)'. See also: *note erfcinv: XREFerfcinv, *note erfcx: XREFerfcx, *note erfi: XREFerfi, *note dawson: XREFdawson, *note erf: XREFerf, *note erfinv: XREFerfinv. -- Mapping Function: erfcx (Z) Compute the scaled complementary error function. The scaled complementary error function is defined as exp (z^2) * erfc (z) See also: *note erfc: XREFerfc, *note erf: XREFerf, *note erfi: XREFerfi, *note dawson: XREFdawson, *note erfinv: XREFerfinv, *note erfcinv: XREFerfcinv. -- Mapping Function: erfi (Z) Compute the imaginary error function. The imaginary error function is defined as -i * erf (i*z) See also: *note erfc: XREFerfc, *note erf: XREFerf, *note erfcx: XREFerfcx, *note dawson: XREFdawson, *note erfinv: XREFerfinv, *note erfcinv: XREFerfcinv. -- Mapping Function: erfinv (X) Compute the inverse error function. The inverse error function is defined such that erf (Y) == X See also: *note erf: XREFerf, *note erfc: XREFerfc, *note erfcx: XREFerfcx, *note erfi: XREFerfi, *note dawson: XREFdawson, *note erfcinv: XREFerfcinv. -- Mapping Function: erfcinv (X) Compute the inverse complementary error function. The inverse complementary error function is defined such that erfc (Y) == X See also: *note erfc: XREFerfc, *note erf: XREFerf, *note erfcx: XREFerfcx, *note erfi: XREFerfi, *note dawson: XREFdawson, *note erfinv: XREFerfinv. -- Function File: expint (X) Compute the exponential integral: infinity / E_1 (x) = | exp (-t)/t dt / x Note: For compatibility, this functions uses the MATLAB definition of the exponential integral. Most other sources refer to this particular value as E_1 (x), and the exponential integral as infinity / Ei (x) = - | exp (-t)/t dt / -x The two definitions are related, for positive real values of X, by 'E_1 (-x) = -Ei (x) - i*pi'. -- Mapping Function: gamma (Z) Compute the Gamma function. The Gamma function is defined as infinity / gamma (z) = | t^(z-1) exp (-t) dt. / t=0 Programming Note: The gamma function can grow quite large even for small input values. In many cases it may be preferable to use the natural logarithm of the gamma function ('gammaln') in calculations to minimize loss of precision. The final result is then 'exp (RESULT_USING_GAMMALN).' See also: *note gammainc: XREFgammainc, *note gammaln: XREFgammaln, *note factorial: XREFfactorial. -- Mapping Function: gammainc (X, A) -- Mapping Function: gammainc (X, A, "lower") -- Mapping Function: gammainc (X, A, "upper") Compute the normalized incomplete gamma function. This is defined as x 1 / gammainc (x, a) = --------- | exp (-t) t^(a-1) dt gamma (a) / t=0 with the limiting value of 1 as X approaches infinity. The standard notation is P(a,x), e.g., Abramowitz and Stegun (6.5.1). If A is scalar, then 'gammainc (X, A)' is returned for each element of X and vice versa. If neither X nor A is scalar, the sizes of X and A must agree, and 'gammainc' is applied element-by-element. By default the incomplete gamma function integrated from 0 to X is computed. If "upper" is given then the complementary function integrated from X to infinity is calculated. It should be noted that gammainc (X, A) == 1 - gammainc (X, A, "upper") See also: *note gamma: XREFgamma, *note gammaln: XREFgammaln. -- Function File: L = legendre (N, X) -- Function File: L = legendre (N, X, NORMALIZATION) Compute the Legendre function of degree N and order M = 0 ... N. The value N must be a real non-negative integer. X is a vector with real-valued elements in the range [-1, 1]. The optional argument NORMALIZATION may be one of "unnorm", "sch", or "norm". The default if no normalization is given is "unnorm". When the optional argument NORMALIZATION is "unnorm", compute the Legendre function of degree N and order M and return all values for M = 0 ... N. The return value has one dimension more than X. The Legendre Function of degree N and order M: m m 2 m/2 d^m P(x) = (-1) * (1-x ) * ---- P(x) n dx^m n with Legendre polynomial of degree N: 1 d^n 2 n P(x) = ------ [----(x - 1) ] n 2^n n! dx^n 'legendre (3, [-1.0, -0.9, -0.8])' returns the matrix: x | -1.0 | -0.9 | -0.8 ------------------------------------ m=0 | -1.00000 | -0.47250 | -0.08000 m=1 | 0.00000 | -1.99420 | -1.98000 m=2 | 0.00000 | -2.56500 | -4.32000 m=3 | 0.00000 | -1.24229 | -3.24000 When the optional argument 'normalization' is "sch", compute the Schmidt semi-normalized associated Legendre function. The Schmidt semi-normalized associated Legendre function is related to the unnormalized Legendre functions by the following: For Legendre functions of degree N and order 0: 0 0 SP(x) = P(x) n n For Legendre functions of degree n and order m: m m m 2(n-m)! 0.5 SP(x) = P(x) * (-1) * [-------] n n (n+m)! When the optional argument NORMALIZATION is "norm", compute the fully normalized associated Legendre function. The fully normalized associated Legendre function is related to the unnormalized Legendre functions by the following: For Legendre functions of degree N and order M m m m (n+0.5)(n-m)! 0.5 NP(x) = P(x) * (-1) * [-------------] n n (n+m)! -- Mapping Function: gammaln (X) -- Mapping Function: lgamma (X) Return the natural logarithm of the gamma function of X. See also: *note gamma: XREFgamma, *note gammainc: XREFgammainc.  File: octave.info, Node: Rational Approximations, Next: Coordinate Transformations, Prev: Special Functions, Up: Arithmetic 17.7 Rational Approximations ============================ -- Function File: S = rat (X, TOL) -- Function File: [N, D] = rat (X, TOL) Find a rational approximation to X within the tolerance defined by TOL using a continued fraction expansion. For example: rat (pi) = 3 + 1/(7 + 1/16) = 355/113 rat (e) = 3 + 1/(-4 + 1/(2 + 1/(5 + 1/(-2 + 1/(-7))))) = 1457/536 When called with two output arguments return the numerator and denominator separately as two matrices. See also: *note rats: XREFrats. -- Built-in Function: rats (X, LEN) Convert X into a rational approximation represented as a string. The string can be converted back into a matrix as follows: r = rats (hilb (4)); x = str2num (r) The optional second argument defines the maximum length of the string representing the elements of X. By default LEN is 9. If the length of the smallest possible rational approximation exceeds LEN, an asterisk (*) padded with spaces will be returned instead. See also: *note format: XREFformat, *note rat: XREFrat.  File: octave.info, Node: Coordinate Transformations, Next: Mathematical Constants, Prev: Rational Approximations, Up: Arithmetic 17.8 Coordinate Transformations =============================== -- Function File: [THETA, R] = cart2pol (X, Y) -- Function File: [THETA, R, Z] = cart2pol (X, Y, Z) -- Function File: [THETA, R] = cart2pol (C) -- Function File: [THETA, R, Z] = cart2pol (C) -- Function File: P = cart2pol (...) Transform Cartesian coordinates to polar or cylindrical coordinates. The inputs X, Y (, and Z) must be the same shape, or scalar. If called with a single matrix argument then each row of C represents the Cartesian coordinate (X, Y (, Z)). THETA describes the angle relative to the positive x-axis. R is the distance to the z-axis (0, 0, z). If only a single return argument is requested then return a matrix P where each row represents one polar/(cylindrical) coordinate (THETA, PHI (, Z)). See also: *note pol2cart: XREFpol2cart, *note cart2sph: XREFcart2sph, *note sph2cart: XREFsph2cart. -- Function File: [X, Y] = pol2cart (THETA, R) -- Function File: [X, Y, Z] = pol2cart (THETA, R, Z) -- Function File: [X, Y] = pol2cart (P) -- Function File: [X, Y, Z] = pol2cart (P) -- Function File: C = pol2cart (...) Transform polar or cylindrical coordinates to Cartesian coordinates. The inputs THETA, R, (and Z) must be the same shape, or scalar. If called with a single matrix argument then each row of P represents the polar/(cylindrical) coordinate (THETA, R (, Z)). THETA describes the angle relative to the positive x-axis. R is the distance to the z-axis (0, 0, z). If only a single return argument is requested then return a matrix C where each row represents one Cartesian coordinate (X, Y (, Z)). See also: *note cart2pol: XREFcart2pol, *note sph2cart: XREFsph2cart, *note cart2sph: XREFcart2sph. -- Function File: [THETA, PHI, R] = cart2sph (X, Y, Z) -- Function File: [THETA, PHI, R] = cart2sph (C) -- Function File: S = cart2sph (...) Transform Cartesian coordinates to spherical coordinates. The inputs X, Y, and Z must be the same shape, or scalar. If called with a single matrix argument then each row of C represents the Cartesian coordinate (X, Y, Z). THETA describes the angle relative to the positive x-axis. PHI is the angle relative to the xy-plane. R is the distance to the origin (0, 0, 0). If only a single return argument is requested then return a matrix S where each row represents one spherical coordinate (THETA, PHI, R). See also: *note sph2cart: XREFsph2cart, *note cart2pol: XREFcart2pol, *note pol2cart: XREFpol2cart. -- Function File: [X, Y, Z] = sph2cart (THETA, PHI, R) -- Function File: [X, Y, Z] = sph2cart (S) -- Function File: C = sph2cart (...) Transform spherical coordinates to Cartesian coordinates. The inputs THETA, PHI, and R must be the same shape, or scalar. If called with a single matrix argument then each row of S represents the spherical coordinate (THETA, PHI, R). THETA describes the angle relative to the positive x-axis. PHI is the angle relative to the xy-plane. R is the distance to the origin (0, 0, 0). If only a single return argument is requested then return a matrix C where each row represents one Cartesian coordinate (X, Y, Z). See also: *note cart2sph: XREFcart2sph, *note pol2cart: XREFpol2cart, *note cart2pol: XREFcart2pol.  File: octave.info, Node: Mathematical Constants, Prev: Coordinate Transformations, Up: Arithmetic 17.9 Mathematical Constants =========================== -- Built-in Function: e -- Built-in Function: e (N) -- Built-in Function: e (N, M) -- Built-in Function: e (N, M, K, ...) -- Built-in Function: e (..., CLASS) Return a scalar, matrix, or N-dimensional array whose elements are all equal to the base of natural logarithms. The constant 'e' satisfies the equation 'log' (e) = 1. When called with no arguments, return a scalar with the value e. When called with a single argument, return a square matrix with the dimension specified. When called with more than one scalar argument the first two arguments are taken as the number of rows and columns and any further arguments specify additional matrix dimensions. The optional argument CLASS specifies the return type and may be either "double" or "single". See also: *note log: XREFlog, *note exp: XREFexp, *note pi: XREFpi, *note I: XREFI. -- Built-in Function: pi -- Built-in Function: pi (N) -- Built-in Function: pi (N, M) -- Built-in Function: pi (N, M, K, ...) -- Built-in Function: pi (..., CLASS) Return a scalar, matrix, or N-dimensional array whose elements are all equal to the ratio of the circumference of a circle to its diameter. Internally, 'pi' is computed as '4.0 * atan (1.0)'. When called with no arguments, return a scalar with the value of pi. When called with a single argument, return a square matrix with the dimension specified. When called with more than one scalar argument the first two arguments are taken as the number of rows and columns and any further arguments specify additional matrix dimensions. The optional argument CLASS specifies the return type and may be either "double" or "single". See also: *note e: XREFe, *note I: XREFI. -- Built-in Function: I -- Built-in Function: I (N) -- Built-in Function: I (N, M) -- Built-in Function: I (N, M, K, ...) -- Built-in Function: I (..., CLASS) Return a scalar, matrix, or N-dimensional array whose elements are all equal to the pure imaginary unit, defined as 'sqrt (-1)'. I, and its equivalents i, j, and J, are functions so any of the names may be reused for other purposes (such as i for a counter variable). When called with no arguments, return a scalar with the value i. When called with a single argument, return a square matrix with the dimension specified. When called with more than one scalar argument the first two arguments are taken as the number of rows and columns and any further arguments specify additional matrix dimensions. The optional argument CLASS specifies the return type and may be either "double" or "single". See also: *note e: XREFe, *note pi: XREFpi, *note log: XREFlog, *note exp: XREFexp. -- Built-in Function: Inf -- Built-in Function: Inf (N) -- Built-in Function: Inf (N, M) -- Built-in Function: Inf (N, M, K, ...) -- Built-in Function: Inf (..., CLASS) Return a scalar, matrix or N-dimensional array whose elements are all equal to the IEEE representation for positive infinity. Infinity is produced when results are too large to be represented using the IEEE floating point format for numbers. Two common examples which produce infinity are division by zero and overflow. [ 1/0 e^800 ] => Inf Inf When called with no arguments, return a scalar with the value 'Inf'. When called with a single argument, return a square matrix with the dimension specified. When called with more than one scalar argument the first two arguments are taken as the number of rows and columns and any further arguments specify additional matrix dimensions. The optional argument CLASS specifies the return type and may be either "double" or "single". See also: *note isinf: XREFisinf, *note NaN: XREFNaN. -- Built-in Function: NaN -- Built-in Function: NaN (N) -- Built-in Function: NaN (N, M) -- Built-in Function: NaN (N, M, K, ...) -- Built-in Function: NaN (..., CLASS) Return a scalar, matrix, or N-dimensional array whose elements are all equal to the IEEE symbol NaN (Not a Number). NaN is the result of operations which do not produce a well defined numerical result. Common operations which produce a NaN are arithmetic with infinity (Inf - Inf), zero divided by zero (0/0), and any operation involving another NaN value (5 + NaN). Note that NaN always compares not equal to NaN (NaN != NaN). This behavior is specified by the IEEE standard for floating point arithmetic. To find NaN values, use the 'isnan' function. When called with no arguments, return a scalar with the value 'NaN'. When called with a single argument, return a square matrix with the dimension specified. When called with more than one scalar argument the first two arguments are taken as the number of rows and columns and any further arguments specify additional matrix dimensions. The optional argument CLASS specifies the return type and may be either "double" or "single". See also: *note isnan: XREFisnan, *note Inf: XREFInf. -- Built-in Function: eps -- Built-in Function: eps (X) -- Built-in Function: eps (N, M) -- Built-in Function: eps (N, M, K, ...) -- Built-in Function: eps (..., CLASS) Return a scalar, matrix or N-dimensional array whose elements are all eps, the machine precision. More precisely, 'eps' is the relative spacing between any two adjacent numbers in the machine's floating point system. This number is obviously system dependent. On machines that support IEEE floating point arithmetic, 'eps' is approximately 2.2204e-16 for double precision and 1.1921e-07 for single precision. When called with no arguments, return a scalar with the value 'eps (1.0)'. Given a single argument X, return the distance between X and the next largest value. When called with more than one argument the first two arguments are taken as the number of rows and columns and any further arguments specify additional matrix dimensions. The optional argument CLASS specifies the return type and may be either "double" or "single". See also: *note realmax: XREFrealmax, *note realmin: XREFrealmin, *note intmax: XREFintmax, *note bitmax: XREFbitmax. -- Built-in Function: realmax -- Built-in Function: realmax (N) -- Built-in Function: realmax (N, M) -- Built-in Function: realmax (N, M, K, ...) -- Built-in Function: realmax (..., CLASS) Return a scalar, matrix, or N-dimensional array whose elements are all equal to the largest floating point number that is representable. The actual value is system dependent. On machines that support IEEE floating point arithmetic, 'realmax' is approximately 1.7977e+308 for double precision and 3.4028e+38 for single precision. When called with no arguments, return a scalar with the value 'realmax ("double")'. When called with a single argument, return a square matrix with the dimension specified. When called with more than one scalar argument the first two arguments are taken as the number of rows and columns and any further arguments specify additional matrix dimensions. The optional argument CLASS specifies the return type and may be either "double" or "single". See also: *note realmin: XREFrealmin, *note intmax: XREFintmax, *note bitmax: XREFbitmax, *note eps: XREFeps. -- Built-in Function: realmin -- Built-in Function: realmin (N) -- Built-in Function: realmin (N, M) -- Built-in Function: realmin (N, M, K, ...) -- Built-in Function: realmin (..., CLASS) Return a scalar, matrix, or N-dimensional array whose elements are all equal to the smallest normalized floating point number that is representable. The actual value is system dependent. On machines that support IEEE floating point arithmetic, 'realmin' is approximately 2.2251e-308 for double precision and 1.1755e-38 for single precision. When called with no arguments, return a scalar with the value 'realmin ("double")'. When called with a single argument, return a square matrix with the dimension specified. When called with more than one scalar argument the first two arguments are taken as the number of rows and columns and any further arguments specify additional matrix dimensions. The optional argument CLASS specifies the return type and may be either "double" or "single". See also: *note realmax: XREFrealmax, *note intmin: XREFintmin, *note eps: XREFeps.  File: octave.info, Node: Linear Algebra, Next: Vectorization and Faster Code Execution, Prev: Arithmetic, Up: Top 18 Linear Algebra ***************** This chapter documents the linear algebra functions provided in Octave. Reference material for many of these functions may be found in Golub and Van Loan, 'Matrix Computations, 2nd Ed.', Johns Hopkins, 1989, and in the 'LAPACK Users' Guide', SIAM, 1992. The 'LAPACK Users' Guide' is available at: 'http://www.netlib.org/lapack/lug/' A common text for engineering courses is G. Strang, 'Linear Algebra and Its Applications, 4th Edition'. It has become a widespread reference for linear algebra. An alternative is P. Lax 'Linear Algebra and Its Applications', and also is a good choice. It claims to be suitable for high school students with substantial mathematical interests as well as first-year undergraduates. * Menu: * Techniques Used for Linear Algebra:: * Basic Matrix Functions:: * Matrix Factorizations:: * Functions of a Matrix:: * Specialized Solvers::  File: octave.info, Node: Techniques Used for Linear Algebra, Next: Basic Matrix Functions, Up: Linear Algebra 18.1 Techniques Used for Linear Algebra ======================================= Octave includes a polymorphic solver that selects an appropriate matrix factorization depending on the properties of the matrix itself. Generally, the cost of determining the matrix type is small relative to the cost of factorizing the matrix itself. In any case the matrix type is cached once it is calculated so that it is not re-determined each time it is used in a linear equation. The selection tree for how the linear equation is solved or a matrix inverse is formed is given by: 1. If the matrix is upper or lower triangular sparse use a forward or backward substitution using the LAPACK xTRTRS function, and goto 4. 2. If the matrix is square, Hermitian with a real positive diagonal, attempt Cholesky factorization using the LAPACK xPOTRF function. 3. If the Cholesky factorization failed or the matrix is not Hermitian with a real positive diagonal, and the matrix is square, factorize using the LAPACK xGETRF function. 4. If the matrix is not square, or any of the previous solvers flags a singular or near singular matrix, find a least squares solution using the LAPACK xGELSD function. The user can force the type of the matrix with the 'matrix_type' function. This overcomes the cost of discovering the type of the matrix. However, it should be noted that identifying the type of the matrix incorrectly will lead to unpredictable results, and so 'matrix_type' should be used with care. It should be noted that the test for whether a matrix is a candidate for Cholesky factorization, performed above, and by the 'matrix_type' function, does not make certain that the matrix is Hermitian. However, the attempt to factorize the matrix will quickly detect a non-Hermitian matrix.  File: octave.info, Node: Basic Matrix Functions, Next: Matrix Factorizations, Prev: Techniques Used for Linear Algebra, Up: Linear Algebra 18.2 Basic Matrix Functions =========================== -- Built-in Function: AA = balance (A) -- Built-in Function: AA = balance (A, OPT) -- Built-in Function: [DD, AA] = balance (A, OPT) -- Built-in Function: [D, P, AA] = balance (A, OPT) -- Built-in Function: [CC, DD, AA, BB] = balance (A, B, OPT) Balance the matrix A to reduce numerical errors in future calculations. Compute 'AA = DD \ A * DD' in which AA is a matrix whose row and column norms are roughly equal in magnitude, and 'DD = P * D', in which P is a permutation matrix and D is a diagonal matrix of powers of two. This allows the equilibration to be computed without round-off. Results of eigenvalue calculation are typically improved by balancing first. If two output values are requested, 'balance' returns the diagonal D and the permutation P separately as vectors. In this case, 'DD = eye(n)(:,P) * diag (D)', where n is the matrix size. If four output values are requested, compute 'AA = CC*A*DD' and 'BB = CC*B*DD', in which AA and BB have nonzero elements of approximately the same magnitude and CC and DD are permuted diagonal matrices as in DD for the algebraic eigenvalue problem. The eigenvalue balancing option OPT may be one of: "noperm", "S" Scale only; do not permute. "noscal", "P" Permute only; do not scale. Algebraic eigenvalue balancing uses standard LAPACK routines. Generalized eigenvalue problem balancing uses Ward's algorithm (SIAM Journal on Scientific and Statistical Computing, 1981). -- Function File: BW = bandwidth (A, TYPE) -- Function File: [LOWER, UPPER] = bandwidth (A) Compute the bandwidth of A. The TYPE argument is the string "lower" for the lower bandwidth and "upper" for the upper bandwidth. If no TYPE is specified return both the lower and upper bandwidth of A. The lower/upper bandwidth of a matrix is the number of subdiagonals/superdiagonals with nonzero entries. See also: *note isbanded: XREFisbanded, *note isdiag: XREFisdiag, *note istril: XREFistril, *note istriu: XREFistriu. -- Function File: cond (A) -- Function File: cond (A, P) Compute the P-norm condition number of a matrix. 'cond (A)' is defined as 'norm (A, P) * norm (inv (A), P)'. By default, 'P = 2' is used which implies a (relatively slow) singular value decomposition. Other possible selections are 'P = 1, Inf, "fro"' which are generally faster. See 'norm' for a full discussion of possible P values. The condition number of a matrix quantifies the sensitivity of the matrix inversion operation when small changes are made to matrix elements. Ideally the condition number will be close to 1. When the number is large this indicates small changes (such as underflow or round-off error) will produce large changes in the resulting output. In such cases the solution results from numerical computing are not likely to be accurate. See also: *note condest: XREFcondest, *note rcond: XREFrcond, *note norm: XREFnorm, *note svd: XREFsvd. -- Built-in Function: det (A) -- Built-in Function: [D, RCOND] = det (A) Compute the determinant of A. Return an estimate of the reciprocal condition number if requested. Programming Notes: Routines from LAPACK are used for full matrices and code from UMFPACK is used for sparse matrices. The determinant should not be used to check a matrix for singularity. For that, use any of the condition number functions: 'cond', 'condest', 'rcond'. See also: *note cond: XREFcond, *note condest: XREFcondest, *note rcond: XREFrcond. -- Built-in Function: LAMBDA = eig (A) -- Built-in Function: LAMBDA = eig (A, B) -- Built-in Function: [V, LAMBDA] = eig (A) -- Built-in Function: [V, LAMBDA] = eig (A, B) Compute the eigenvalues (and optionally the eigenvectors) of a matrix or a pair of matrices The algorithm used depends on whether there are one or two input matrices, if they are real or complex, and if they are symmetric (Hermitian if complex) or non-symmetric. The eigenvalues returned by 'eig' are not ordered. See also: *note eigs: XREFeigs, *note svd: XREFsvd. -- Built-in Function: G = givens (X, Y) -- Built-in Function: [C, S] = givens (X, Y) Compute the Givens rotation matrix G. The Givens matrix is a 2 by 2 orthogonal matrix 'G = [C S; -S' C]' such that 'G [X; Y] = [*; 0]' with X and Y scalars. If two output arguments are requested, return the factors C and S rather than the Givens rotation matrix. For example: givens (1, 1) => 0.70711 0.70711 -0.70711 0.70711 See also: *note planerot: XREFplanerot. -- Function File: [G, Y] = planerot (X) Given a two-element column vector, return the 2 by 2 orthogonal matrix G such that 'Y = G * X' and 'Y(2) = 0'. See also: *note givens: XREFgivens. -- Built-in Function: X = inv (A) -- Built-in Function: [X, RCOND] = inv (A) Compute the inverse of the square matrix A. Return an estimate of the reciprocal condition number if requested, otherwise warn of an ill-conditioned matrix if the reciprocal condition number is small. In general it is best to avoid calculating the inverse of a matrix directly. For example, it is both faster and more accurate to solve systems of equations (A*x = b) with 'Y = A \ b', rather than 'Y = inv (A) * b'. If called with a sparse matrix, then in general X will be a full matrix requiring significantly more storage. Avoid forming the inverse of a sparse matrix if possible. See also: *note ldivide: XREFldivide, *note rdivide: XREFrdivide. -- Function File: X = linsolve (A, B) -- Function File: X = linsolve (A, B, OPTS) -- Function File: [X, R] = linsolve (...) Solve the linear system 'A*x = b'. With no options, this function is equivalent to the left division operator ('x = A \ b') or the matrix-left-divide function ('x = mldivide (A, b)'). Octave ordinarily examines the properties of the matrix A and chooses a solver that best matches the matrix. By passing a structure OPTS to 'linsolve' you can inform Octave directly about the matrix A. In this case Octave will skip the matrix examination and proceed directly to solving the linear system. *Warning:* If the matrix A does not have the properties listed in the OPTS structure then the result will not be accurate AND no warning will be given. When in doubt, let Octave examine the matrix and choose the appropriate solver as this step takes little time and the result is cached so that it is only done once per linear system. Possible OPTS fields (set value to true/false): LT A is lower triangular UT A is upper triangular UHESS A is upper Hessenberg (currently makes no difference) SYM A is symmetric or complex Hermitian (currently makes no difference) POSDEF A is positive definite RECT A is general rectangular (currently makes no difference) TRANSA Solve 'A'*x = b' by 'transpose (A) \ b' The optional second output R is the inverse condition number of A (zero if matrix is singular). See also: *note mldivide: XREFmldivide, *note matrix_type: XREFmatrix_type, *note rcond: XREFrcond. -- Built-in Function: TYPE = matrix_type (A) -- Built-in Function: TYPE = matrix_type (A, "nocompute") -- Built-in Function: A = matrix_type (A, TYPE) -- Built-in Function: A = matrix_type (A, "upper", PERM) -- Built-in Function: A = matrix_type (A, "lower", PERM) -- Built-in Function: A = matrix_type (A, "banded", NL, NU) Identify the matrix type or mark a matrix as a particular type. This allows more rapid solutions of linear equations involving A to be performed. Called with a single argument, 'matrix_type' returns the type of the matrix and caches it for future use. Called with more than one argument, 'matrix_type' allows the type of the matrix to be defined. If the option "nocompute" is given, the function will not attempt to guess the type if it is still unknown. This is useful for debugging purposes. The possible matrix types depend on whether the matrix is full or sparse, and can be one of the following "unknown" Remove any previously cached matrix type, and mark type as unknown. "full" Mark the matrix as full. "positive definite" Probable full positive definite matrix. "diagonal" Diagonal matrix. (Sparse matrices only) "permuted diagonal" Permuted Diagonal matrix. The permutation does not need to be specifically indicated, as the structure of the matrix explicitly gives this. (Sparse matrices only) "upper" Upper triangular. If the optional third argument PERM is given, the matrix is assumed to be a permuted upper triangular with the permutations defined by the vector PERM. "lower" Lower triangular. If the optional third argument PERM is given, the matrix is assumed to be a permuted lower triangular with the permutations defined by the vector PERM. "banded" "banded positive definite" Banded matrix with the band size of NL below the diagonal and NU above it. If NL and NU are 1, then the matrix is tridiagonal and treated with specialized code. In addition the matrix can be marked as probably a positive definite. (Sparse matrices only) "singular" The matrix is assumed to be singular and will be treated with a minimum norm solution. Note that the matrix type will be discovered automatically on the first attempt to solve a linear equation involving A. Therefore 'matrix_type' is only useful to give Octave hints of the matrix type. Incorrectly defining the matrix type will result in incorrect results from solutions of linear equations; it is entirely *the responsibility of the user* to correctly identify the matrix type. Also, the test for positive definiteness is a low-cost test for a Hermitian matrix with a real positive diagonal. This does not guarantee that the matrix is positive definite, but only that it is a probable candidate. When such a matrix is factorized, a Cholesky factorization is first attempted, and if that fails the matrix is then treated with an LU factorization. Once the matrix has been factorized, 'matrix_type' will return the correct classification of the matrix. -- Built-in Function: norm (A) -- Built-in Function: norm (A, P) -- Built-in Function: norm (A, P, OPT) Compute the p-norm of the matrix A. If the second argument is missing, 'p = 2' is assumed. If A is a matrix (or sparse matrix): P = '1' 1-norm, the largest column sum of the absolute values of A. P = '2' Largest singular value of A. P = 'Inf' or "inf" Infinity norm, the largest row sum of the absolute values of A. P = "fro" Frobenius norm of A, 'sqrt (sum (diag (A' * A)))'. other P, 'P > 1' maximum 'norm (A*x, p)' such that 'norm (x, p) == 1' If A is a vector or a scalar: P = 'Inf' or "inf" 'max (abs (A))'. P = '-Inf' 'min (abs (A))'. P = "fro" Frobenius norm of A, 'sqrt (sumsq (abs (A)))'. P = 0 Hamming norm - the number of nonzero elements. other P, 'P > 1' p-norm of A, '(sum (abs (A) .^ P)) ^ (1/P)'. other P 'P < 1' the p-pseudonorm defined as above. If OPT is the value "rows", treat each row as a vector and compute its norm. The result is returned as a column vector. Similarly, if OPT is "columns" or "cols" then compute the norms of each column and return a row vector. See also: *note cond: XREFcond, *note svd: XREFsvd. -- Function File: null (A) -- Function File: null (A, TOL) Return an orthonormal basis of the null space of A. The dimension of the null space is taken as the number of singular values of A not greater than TOL. If the argument TOL is missing, it is computed as max (size (A)) * max (svd (A)) * eps See also: *note orth: XREForth. -- Function File: orth (A) -- Function File: orth (A, TOL) Return an orthonormal basis of the range space of A. The dimension of the range space is taken as the number of singular values of A greater than TOL. If the argument TOL is missing, it is computed as max (size (A)) * max (svd (A)) * eps See also: *note null: XREFnull. -- Built-in Function: [Y, H] = mgorth (X, V) Orthogonalize a given column vector X with respect to a set of orthonormal vectors comprising the columns of V using the modified Gram-Schmidt method. On exit, Y is a unit vector such that: norm (Y) = 1 V' * Y = 0 X = [V, Y]*H' -- Built-in Function: pinv (X) -- Built-in Function: pinv (X, TOL) Return the pseudoinverse of X. Singular values less than TOL are ignored. If the second argument is omitted, it is taken to be tol = max (size (X)) * sigma_max (X) * eps, where 'sigma_max (X)' is the maximal singular value of X. -- Function File: rank (A) -- Function File: rank (A, TOL) Compute the rank of matrix A, using the singular value decomposition. The rank is taken to be the number of singular values of A that are greater than the specified tolerance TOL. If the second argument is omitted, it is taken to be tol = max (size (A)) * sigma(1) * eps; where 'eps' is machine precision and 'sigma(1)' is the largest singular value of A. The rank of a matrix is the number of linearly independent rows or columns and determines how many particular solutions exist to a system of equations. Use 'null' for finding the remaining homogenous solutions. Example: x = [1 2 3 4 5 6 7 8 9]; rank (x) => 2 The number of linearly independent rows is only 2 because the final row is a linear combination of -1*row1 + 2*row2. See also: *note null: XREFnull, *note sprank: XREFsprank, *note svd: XREFsvd. -- Built-in Function: C = rcond (A) Compute the 1-norm estimate of the reciprocal condition number as returned by LAPACK. If the matrix is well-conditioned then C will be near 1 and if the matrix is poorly conditioned it will be close to 0. The matrix A must not be sparse. If the matrix is sparse then 'condest (A)' or 'rcond (full (A))' should be used instead. See also: *note cond: XREFcond, *note condest: XREFcondest. -- Function File: trace (A) Compute the trace of A, the sum of the elements along the main diagonal. The implementation is straightforward: 'sum (diag (A))'. See also: *note eig: XREFeig. -- Function File: rref (A) -- Function File: rref (A, TOL) -- Function File: [R, K] = rref (...) Return the reduced row echelon form of A. TOL defaults to 'eps * max (size (A)) * norm (A, inf)'. The optional return argument K contains the vector of "bound variables", which are those columns on which elimination has been performed.  File: octave.info, Node: Matrix Factorizations, Next: Functions of a Matrix, Prev: Basic Matrix Functions, Up: Linear Algebra 18.3 Matrix Factorizations ========================== -- Loadable Function: R = chol (A) -- Loadable Function: [R, P] = chol (A) -- Loadable Function: [R, P, Q] = chol (S) -- Loadable Function: [R, P, Q] = chol (S, "vector") -- Loadable Function: [L, ...] = chol (..., "lower") -- Loadable Function: [L, ...] = chol (..., "upper") Compute the Cholesky factor, R, of the symmetric positive definite matrix A. The Cholesky factor is defined by R' * R = A. Called with one output argument 'chol' fails if A or S is not positive definite. With two or more output arguments P flags whether the matrix was positive definite and 'chol' does not fail. A zero value indicated that the matrix was positive definite and the R gives the factorization, and P will have a positive value otherwise. If called with 3 outputs then a sparsity preserving row/column permutation is applied to A prior to the factorization. That is R is the factorization of 'A(Q,Q)' such that R' * R = Q' * A * Q. The sparsity preserving permutation is generally returned as a matrix. However, given the flag "vector", Q will be returned as a vector such that R' * R = A(Q, Q). Called with either a sparse or full matrix and using the "lower" flag, 'chol' returns the lower triangular factorization such that L * L' = A. For full matrices, if the "lower" flag is set only the lower triangular part of the matrix is used for the factorization, otherwise the upper triangular part is used. In general the lower triangular factorization is significantly faster for sparse matrices. See also: *note hess: XREFhess, *note lu: XREFlu, *note qr: XREFqr, *note qz: XREFqz, *note schur: XREFschur, *note svd: XREFsvd, *note ichol: XREFichol, *note cholinv: XREFcholinv, *note chol2inv: XREFchol2inv, *note cholupdate: XREFcholupdate, *note cholinsert: XREFcholinsert, *note choldelete: XREFcholdelete, *note cholshift: XREFcholshift. -- Loadable Function: cholinv (A) Compute the inverse of the symmetric positive definite matrix A using the Cholesky factorization. See also: *note chol: XREFchol, *note chol2inv: XREFchol2inv, *note inv: XREFinv. -- Loadable Function: chol2inv (U) Invert a symmetric, positive definite square matrix from its Cholesky decomposition, U. Note that U should be an upper-triangular matrix with positive diagonal elements. 'chol2inv (U)' provides 'inv (U'*U)' but it is much faster than using 'inv'. See also: *note chol: XREFchol, *note cholinv: XREFcholinv, *note inv: XREFinv. -- Loadable Function: [R1, INFO] = cholupdate (R, U, OP) Update or downdate a Cholesky factorization. Given an upper triangular matrix R and a column vector U, attempt to determine another upper triangular matrix R1 such that * R1'*R1 = R'*R + U*U' if OP is "+" * R1'*R1 = R'*R - U*U' if OP is "-" If OP is "-", INFO is set to * 0 if the downdate was successful, * 1 if R'*R - U*U' is not positive definite, * 2 if R is singular. If INFO is not present, an error message is printed in cases 1 and 2. See also: *note chol: XREFchol, *note cholinsert: XREFcholinsert, *note choldelete: XREFcholdelete, *note cholshift: XREFcholshift. -- Loadable Function: R1 = cholinsert (R, J, U) -- Loadable Function: [R1, INFO] = cholinsert (R, J, U) Given a Cholesky factorization of a real symmetric or complex Hermitian positive definite matrix A = R'*R, R upper triangular, return the Cholesky factorization of A1, where A1(p,p) = A, A1(:,j) = A1(j,:)' = u and p = [1:j-1,j+1:n+1]. u(j) should be positive. On return, INFO is set to * 0 if the insertion was successful, * 1 if A1 is not positive definite, * 2 if R is singular. If INFO is not present, an error message is printed in cases 1 and 2. See also: *note chol: XREFchol, *note cholupdate: XREFcholupdate, *note choldelete: XREFcholdelete, *note cholshift: XREFcholshift. -- Loadable Function: R1 = choldelete (R, J) Given a Cholesky factorization of a real symmetric or complex Hermitian positive definite matrix A = R'*R, R upper triangular, return the Cholesky factorization of A(p,p), where p = [1:j-1,j+1:n+1]. See also: *note chol: XREFchol, *note cholupdate: XREFcholupdate, *note cholinsert: XREFcholinsert, *note cholshift: XREFcholshift. -- Loadable Function: R1 = cholshift (R, I, J) Given a Cholesky factorization of a real symmetric or complex Hermitian positive definite matrix A = R'*R, R upper triangular, return the Cholesky factorization of A(p,p), where p is the permutation 'p = [1:i-1, shift(i:j, 1), j+1:n]' if I < J or 'p = [1:j-1, shift(j:i,-1), i+1:n]' if J < I. See also: *note chol: XREFchol, *note cholupdate: XREFcholupdate, *note cholinsert: XREFcholinsert, *note choldelete: XREFcholdelete. -- Built-in Function: H = hess (A) -- Built-in Function: [P, H] = hess (A) Compute the Hessenberg decomposition of the matrix A. The Hessenberg decomposition is 'P * H * P' = A' where P is a square unitary matrix ('P' * P = I', using complex-conjugate transposition) and H is upper Hessenberg ('H(i, j) = 0 forall i >= j+1)'. The Hessenberg decomposition is usually used as the first step in an eigenvalue computation, but has other applications as well (see Golub, Nash, and Van Loan, IEEE Transactions on Automatic Control, 1979). See also: *note eig: XREFeig, *note chol: XREFchol, *note lu: XREFlu, *note qr: XREFqr, *note qz: XREFqz, *note schur: XREFschur, *note svd: XREFsvd. -- Built-in Function: [L, U] = lu (A) -- Built-in Function: [L, U, P] = lu (A) -- Built-in Function: [L, U, P, Q] = lu (S) -- Built-in Function: [L, U, P, Q, R] = lu (S) -- Built-in Function: [...] = lu (S, THRES) -- Built-in Function: Y = lu (...) -- Built-in Function: [...] = lu (..., "vector") Compute the LU decomposition of A. If A is full subroutines from LAPACK are used and if A is sparse then UMFPACK is used. The result is returned in a permuted form, according to the optional return value P. For example, given the matrix 'a = [1, 2; 3, 4]', [l, u, p] = lu (A) returns l = 1.00000 0.00000 0.33333 1.00000 u = 3.00000 4.00000 0.00000 0.66667 p = 0 1 1 0 The matrix is not required to be square. When called with two or three output arguments and a spare input matrix, 'lu' does not attempt to perform sparsity preserving column permutations. Called with a fourth output argument, the sparsity preserving column transformation Q is returned, such that 'P * A * Q = L * U'. Called with a fifth output argument and a sparse input matrix, 'lu' attempts to use a scaling factor R on the input matrix such that 'P * (R \ A) * Q = L * U'. This typically leads to a sparser and more stable factorization. An additional input argument THRES, that defines the pivoting threshold can be given. THRES can be a scalar, in which case it defines the UMFPACK pivoting tolerance for both symmetric and unsymmetric cases. If THRES is a 2-element vector, then the first element defines the pivoting tolerance for the unsymmetric UMFPACK pivoting strategy and the second for the symmetric strategy. By default, the values defined by 'spparms' are used ([0.1, 0.001]). Given the string argument "vector", 'lu' returns the values of P and Q as vector values, such that for full matrix, 'A (P,:) = L * U', and 'R(P,:) * A (:, Q) = L * U'. With two output arguments, returns the permuted forms of the upper and lower triangular matrices, such that 'A = L * U'. With one output argument Y, then the matrix returned by the LAPACK routines is returned. If the input matrix is sparse then the matrix L is embedded into U to give a return value similar to the full case. For both full and sparse matrices, 'lu' loses the permutation information. See also: *note luupdate: XREFluupdate, *note ilu: XREFilu, *note chol: XREFchol, *note hess: XREFhess, *note qr: XREFqr, *note qz: XREFqz, *note schur: XREFschur, *note svd: XREFsvd. -- Built-in Function: [L, U] = luupdate (L, U, X, Y) -- Built-in Function: [L, U, P] = luupdate (L, U, P, X, Y) Given an LU factorization of a real or complex matrix A = L*U, L lower unit trapezoidal and U upper trapezoidal, return the LU factorization of A + X*Y.', where X and Y are column vectors (rank-1 update) or matrices with equal number of columns (rank-k update). Optionally, row-pivoted updating can be used by supplying a row permutation (pivoting) matrix P; in that case, an updated permutation matrix is returned. Note that if L, U, P is a pivoted LU factorization as obtained by 'lu': [L, U, P] = lu (A); then a factorization of A+X*Y.' can be obtained either as [L1, U1] = lu (L, U, P*X, Y) or [L1, U1, P1] = lu (L, U, P, X, Y) The first form uses the unpivoted algorithm, which is faster, but less stable. The second form uses a slower pivoted algorithm, which is more stable. The matrix case is done as a sequence of rank-1 updates; thus, for large enough k, it will be both faster and more accurate to recompute the factorization from scratch. See also: *note lu: XREFlu, *note cholupdate: XREFcholupdate, *note qrupdate: XREFqrupdate. -- Loadable Function: [Q, R, P] = qr (A) -- Loadable Function: [Q, R, P] = qr (A, '0') -- Loadable Function: [C, R] = qr (A, B) -- Loadable Function: [C, R] = qr (A, B, '0') Compute the QR factorization of A, using standard LAPACK subroutines. For example, given the matrix 'A = [1, 2; 3, 4]', [Q, R] = qr (A) returns Q = -0.31623 -0.94868 -0.94868 0.31623 R = -3.16228 -4.42719 0.00000 -0.63246 The 'qr' factorization has applications in the solution of least squares problems min norm(A x - b) for overdetermined systems of equations (i.e., A is a tall, thin matrix). The QR factorization is 'Q * R = A' where Q is an orthogonal matrix and R is upper triangular. If given a second argument of '0', 'qr' returns an economy-sized QR factorization, omitting zero rows of R and the corresponding columns of Q. If the matrix A is full, the permuted QR factorization '[Q, R, P] = qr (A)' forms the QR factorization such that the diagonal entries of R are decreasing in magnitude order. For example, given the matrix 'a = [1, 2; 3, 4]', [Q, R, P] = qr (A) returns Q = -0.44721 -0.89443 -0.89443 0.44721 R = -4.47214 -3.13050 0.00000 0.44721 P = 0 1 1 0 The permuted 'qr' factorization '[Q, R, P] = qr (A)' factorization allows the construction of an orthogonal basis of 'span (A)'. If the matrix A is sparse, then compute the sparse QR factorization of A, using CSPARSE. As the matrix Q is in general a full matrix, this function returns the Q-less factorization R of A, such that 'R = chol (A' * A)'. If the final argument is the scalar '0' and the number of rows is larger than the number of columns, then an economy factorization is returned. That is R will have only 'size (A,1)' rows. If an additional matrix B is supplied, then 'qr' returns C, where 'C = Q' * B'. This allows the least squares approximation of 'A \ B' to be calculated as [C, R] = qr (A, B) x = R \ C See also: *note chol: XREFchol, *note hess: XREFhess, *note lu: XREFlu, *note qz: XREFqz, *note schur: XREFschur, *note svd: XREFsvd, *note qrupdate: XREFqrupdate, *note qrinsert: XREFqrinsert, *note qrdelete: XREFqrdelete, *note qrshift: XREFqrshift. -- Loadable Function: [Q1, R1] = qrupdate (Q, R, U, V) Given a QR factorization of a real or complex matrix A = Q*R, Q unitary and R upper trapezoidal, return the QR factorization of A + U*V', where U and V are column vectors (rank-1 update) or matrices with equal number of columns (rank-k update). Notice that the latter case is done as a sequence of rank-1 updates; thus, for k large enough, it will be both faster and more accurate to recompute the factorization from scratch. The QR factorization supplied may be either full (Q is square) or economized (R is square). See also: *note qr: XREFqr, *note qrinsert: XREFqrinsert, *note qrdelete: XREFqrdelete, *note qrshift: XREFqrshift. -- Loadable Function: [Q1, R1] = qrinsert (Q, R, J, X, ORIENT) Given a QR factorization of a real or complex matrix A = Q*R, Q unitary and R upper trapezoidal, return the QR factorization of [A(:,1:j-1) x A(:,j:n)], where U is a column vector to be inserted into A (if ORIENT is "col"), or the QR factorization of [A(1:j-1,:);x;A(:,j:n)], where X is a row vector to be inserted into A (if ORIENT is "row"). The default value of ORIENT is "col". If ORIENT is "col", U may be a matrix and J an index vector resulting in the QR factorization of a matrix B such that B(:,J) gives U and B(:,J) = [] gives A. Notice that the latter case is done as a sequence of k insertions; thus, for k large enough, it will be both faster and more accurate to recompute the factorization from scratch. If ORIENT is "col", the QR factorization supplied may be either full (Q is square) or economized (R is square). If ORIENT is "row", full factorization is needed. See also: *note qr: XREFqr, *note qrupdate: XREFqrupdate, *note qrdelete: XREFqrdelete, *note qrshift: XREFqrshift. -- Loadable Function: [Q1, R1] = qrdelete (Q, R, J, ORIENT) Given a QR factorization of a real or complex matrix A = Q*R, Q unitary and R upper trapezoidal, return the QR factorization of [A(:,1:j-1) A(:,j+1:n)], i.e., A with one column deleted (if ORIENT is "col"), or the QR factorization of [A(1:j-1,:);A(j+1:n,:)], i.e., A with one row deleted (if ORIENT is "row"). The default value of ORIENT is "col". If ORIENT is "col", J may be an index vector resulting in the QR factorization of a matrix B such that A(:,J) = [] gives B. Notice that the latter case is done as a sequence of k deletions; thus, for k large enough, it will be both faster and more accurate to recompute the factorization from scratch. If ORIENT is "col", the QR factorization supplied may be either full (Q is square) or economized (R is square). If ORIENT is "row", full factorization is needed. See also: *note qr: XREFqr, *note qrupdate: XREFqrupdate, *note qrinsert: XREFqrinsert, *note qrshift: XREFqrshift. -- Loadable Function: [Q1, R1] = qrshift (Q, R, I, J) Given a QR factorization of a real or complex matrix A = Q*R, Q unitary and R upper trapezoidal, return the QR factorization of A(:,p), where p is the permutation 'p = [1:i-1, shift(i:j, 1), j+1:n]' if I < J or 'p = [1:j-1, shift(j:i,-1), i+1:n]' if J < I. See also: *note qr: XREFqr, *note qrupdate: XREFqrupdate, *note qrinsert: XREFqrinsert, *note qrdelete: XREFqrdelete. -- Built-in Function: LAMBDA = qz (A, B) -- Built-in Function: LAMBDA = qz (A, B, OPT) QZ decomposition of the generalized eigenvalue problem (A x = s B x). There are three ways to call this function: 1. 'LAMBDA = qz (A, B)' Computes the generalized eigenvalues LAMBDA of (A - s B). 2. '[AA, BB, Q, Z, V, W, LAMBDA] = qz (A, B)' Computes QZ decomposition, generalized eigenvectors, and generalized eigenvalues of (A - s B) A * V = B * V * diag (LAMBDA) W' * A = diag (LAMBDA) * W' * B AA = Q * A * Z, BB = Q * B * Z with Q and Z orthogonal (unitary)= I 3. '[AA,BB,Z{, LAMBDA}] = qz (A, B, OPT)' As in form [2], but allows ordering of generalized eigenpairs for, e.g., solution of discrete time algebraic Riccati equations. Form 3 is not available for complex matrices, and does not compute the generalized eigenvectors V, W, nor the orthogonal matrix Q. OPT for ordering eigenvalues of the GEP pencil. The leading block of the revised pencil contains all eigenvalues that satisfy: "N" = unordered (default) "S" = small: leading block has all |lambda| <= 1 "B" = big: leading block has all |lambda| >= 1 "-" = negative real part: leading block has all eigenvalues in the open left half-plane "+" = non-negative real part: leading block has all eigenvalues in the closed right half-plane Note: 'qz' performs permutation balancing, but not scaling (*note XREFbalance::). The order of output arguments was selected for compatibility with MATLAB. See also: *note eig: XREFeig, *note balance: XREFbalance, *note lu: XREFlu, *note chol: XREFchol, *note hess: XREFhess, *note qr: XREFqr, *note qzhess: XREFqzhess, *note schur: XREFschur, *note svd: XREFsvd. -- Function File: [AA, BB, Q, Z] = qzhess (A, B) Compute the Hessenberg-triangular decomposition of the matrix pencil '(A, B)', returning 'AA = Q * A * Z', 'BB = Q * B * Z', with Q and Z orthogonal. For example: [aa, bb, q, z] = qzhess ([1, 2; 3, 4], [5, 6; 7, 8]) => aa = [ -3.02244, -4.41741; 0.92998, 0.69749 ] => bb = [ -8.60233, -9.99730; 0.00000, -0.23250 ] => q = [ -0.58124, -0.81373; -0.81373, 0.58124 ] => z = [ 1, 0; 0, 1 ] The Hessenberg-triangular decomposition is the first step in Moler and Stewart's QZ decomposition algorithm. Algorithm taken from Golub and Van Loan, 'Matrix Computations, 2nd edition'. See also: *note lu: XREFlu, *note chol: XREFchol, *note hess: XREFhess, *note qr: XREFqr, *note qz: XREFqz, *note schur: XREFschur, *note svd: XREFsvd. -- Built-in Function: S = schur (A) -- Built-in Function: S = schur (A, "real") -- Built-in Function: S = schur (A, "complex") -- Built-in Function: S = schur (A, OPT) -- Built-in Function: [U, S] = schur (...) Compute the Schur decomposition of A. The Schur decomposition is defined as S = U' * A * U where U is a unitary matrix ('U'* U' is identity) and S is upper triangular. The eigenvalues of A (and S) are the diagonal elements of S. If the matrix A is real, then the real Schur decomposition is computed, in which the matrix U is orthogonal and S is block upper triangular with blocks of size at most '2 x 2' along the diagonal. The diagonal elements of S (or the eigenvalues of the '2 x 2' blocks, when appropriate) are the eigenvalues of A and S. The default for real matrices is a real Schur decomposition. A complex decomposition may be forced by passing the flag "complex". The eigenvalues are optionally ordered along the diagonal according to the value of OPT. 'OPT = "a"' indicates that all eigenvalues with negative real parts should be moved to the leading block of S (used in 'are'), 'OPT = "d"' indicates that all eigenvalues with magnitude less than one should be moved to the leading block of S (used in 'dare'), and 'OPT = "u"', the default, indicates that no ordering of eigenvalues should occur. The leading K columns of U always span the A-invariant subspace corresponding to the K leading eigenvalues of S. The Schur decomposition is used to compute eigenvalues of a square matrix, and has applications in the solution of algebraic Riccati equations in control (see 'are' and 'dare'). See also: *note rsf2csf: XREFrsf2csf, *note ordschur: XREFordschur, *note lu: XREFlu, *note chol: XREFchol, *note hess: XREFhess, *note qr: XREFqr, *note qz: XREFqz, *note svd: XREFsvd. -- Function File: [U, T] = rsf2csf (UR, TR) Convert a real, upper quasi-triangular Schur form TR to a complex, upper triangular Schur form T. Note that the following relations hold: UR * TR * UR' = U * T * U' and 'U' * U' is the identity matrix I. Note also that U and T are not unique. See also: *note schur: XREFschur. -- Loadable Function: [UR, SR] = ordschur (U, S, SELECT) Reorders the real Schur factorization (U,S) obtained with the 'schur' function, so that selected eigenvalues appear in the upper left diagonal blocks of the quasi triangular Schur matrix. The logical vector SELECT specifies the selected eigenvalues as they appear along S's diagonal. For example, given the matrix 'A = [1, 2; 3, 4]', and its Schur decomposition [U, S] = schur (A) which returns U = -0.82456 -0.56577 0.56577 -0.82456 S = -0.37228 -1.00000 0.00000 5.37228 It is possible to reorder the decomposition so that the positive eigenvalue is in the upper left corner, by doing: [U, S] = ordschur (U, S, [0,1]) See also: *note schur: XREFschur. -- Function File: ANGLE = subspace (A, B) Determine the largest principal angle between two subspaces spanned by the columns of matrices A and B. -- Built-in Function: S = svd (A) -- Built-in Function: [U, S, V] = svd (A) -- Built-in Function: [U, S, V] = svd (A, ECON) Compute the singular value decomposition of A A = U*S*V' The function 'svd' normally returns only the vector of singular values. When called with three return values, it computes U, S, and V. For example, svd (hilb (3)) returns ans = 1.4083189 0.1223271 0.0026873 and [u, s, v] = svd (hilb (3)) returns u = -0.82704 0.54745 0.12766 -0.45986 -0.52829 -0.71375 -0.32330 -0.64901 0.68867 s = 1.40832 0.00000 0.00000 0.00000 0.12233 0.00000 0.00000 0.00000 0.00269 v = -0.82704 0.54745 0.12766 -0.45986 -0.52829 -0.71375 -0.32330 -0.64901 0.68867 If given a second argument, 'svd' returns an economy-sized decomposition, eliminating the unnecessary rows or columns of U or V. See also: *note svd_driver: XREFsvd_driver, *note svds: XREFsvds, *note eig: XREFeig, *note lu: XREFlu, *note chol: XREFchol, *note hess: XREFhess, *note qr: XREFqr, *note qz: XREFqz. -- Built-in Function: VAL = svd_driver () -- Built-in Function: OLD_VAL = svd_driver (NEW_VAL) -- Built-in Function: svd_driver (NEW_VAL, "local") Query or set the underlying LAPACK driver used by 'svd'. Currently recognized values are "gesvd" and "gesdd". The default is "gesvd". When called from inside a function with the "local" option, the variable is changed locally for the function and any subroutines it calls. The original variable value is restored when exiting the function. See also: *note svd: XREFsvd. -- Function File: [HOUSV, BETA, ZER] = housh (X, J, Z) Compute Householder reflection vector HOUSV to reflect X to be the j-th column of identity, i.e., (I - beta*housv*housv')x = norm (x)*e(j) if x(j) < 0, (I - beta*housv*housv')x = -norm (x)*e(j) if x(j) >= 0 Inputs X vector J index into vector Z threshold for zero (usually should be the number 0) Outputs (see Golub and Van Loan): BETA If beta = 0, then no reflection need be applied (zer set to 0) HOUSV householder vector -- Function File: [U, H, NU] = krylov (A, V, K, EPS1, PFLG) Construct an orthogonal basis U of block Krylov subspace [v a*v a^2*v ... a^(k+1)*v] using Householder reflections to guard against loss of orthogonality. If V is a vector, then H contains the Hessenberg matrix such that a*u == u*h+rk*ek', in which 'rk = a*u(:,k)-u*h(:,k)', and ek' is the vector '[0, 0, ..., 1]' of length 'k'. Otherwise, H is meaningless. If V is a vector and K is greater than 'length (A) - 1', then H contains the Hessenberg matrix such that 'a*u == u*h'. The value of NU is the dimension of the span of the Krylov subspace (based on EPS1). If B is a vector and K is greater than M-1, then H contains the Hessenberg decomposition of A. The optional parameter EPS1 is the threshold for zero. The default value is 1e-12. If the optional parameter PFLG is nonzero, row pivoting is used to improve numerical behavior. The default value is 0. Reference: A. Hodel, P. Misra, 'Partial Pivoting in the Computation of Krylov Subspaces of Large Sparse Systems', Proceedings of the 42nd IEEE Conference on Decision and Control, December 2003.  File: octave.info, Node: Functions of a Matrix, Next: Specialized Solvers, Prev: Matrix Factorizations, Up: Linear Algebra 18.4 Functions of a Matrix ========================== -- Function File: expm (A) Return the exponential of a matrix. The matrix exponential is defined as the infinite Taylor series expm (A) = I + A + A^2/2! + A^3/3! + ... However, the Taylor series is _not_ the way to compute the matrix exponential; see Moler and Van Loan, 'Nineteen Dubious Ways to Compute the Exponential of a Matrix', SIAM Review, 1978. This routine uses Ward's diagonal Pade' approximation method with three step preconditioning (SIAM Journal on Numerical Analysis, 1977). Diagonal Pade' approximations are rational polynomials of matrices -1 D (A) N (A) whose Taylor series matches the first '2q+1' terms of the Taylor series above; direct evaluation of the Taylor series (with the same preconditioning steps) may be desirable in lieu of the Pade' approximation when 'Dq(A)' is ill-conditioned. See also: *note logm: XREFlogm, *note sqrtm: XREFsqrtm. -- Function File: S = logm (A) -- Function File: S = logm (A, OPT_ITERS) -- Function File: [S, ITERS] = logm (...) Compute the matrix logarithm of the square matrix A. The implementation utilizes a Pade' approximant and the identity logm (A) = 2^k * logm (A^(1 / 2^k)) The optional input OPT_ITERS is the maximum number of square roots to compute and defaults to 100. The optional output ITERS is the number of square roots actually computed. See also: *note expm: XREFexpm, *note sqrtm: XREFsqrtm. -- Built-in Function: S = sqrtm (A) -- Built-in Function: [S, ERROR_ESTIMATE] = sqrtm (A) Compute the matrix square root of the square matrix A. Ref: N.J. Higham. 'A New sqrtm for MATLAB'. Numerical Analysis Report No. 336, Manchester Centre for Computational Mathematics, Manchester, England, January 1999. See also: *note expm: XREFexpm, *note logm: XREFlogm. -- Built-in Function: kron (A, B) -- Built-in Function: kron (A1, A2, ...) Form the Kronecker product of two or more matrices. This is defined block by block as x = [ a(i,j)*b ] For example: kron (1:4, ones (3, 1)) => 1 2 3 4 1 2 3 4 1 2 3 4 If there are more than two input arguments A1, A2, ..., AN the Kronecker product is computed as kron (kron (A1, A2), ..., AN) Since the Kronecker product is associative, this is well-defined. -- Built-in Function: blkmm (A, B) Compute products of matrix blocks. The blocks are given as 2-dimensional subarrays of the arrays A, B. The size of A must have the form '[m,k,...]' and size of B must be '[k,n,...]'. The result is then of size '[m,n,...]' and is computed as follows: for i = 1:prod (size (A)(3:end)) C(:,:,i) = A(:,:,i) * B(:,:,i) endfor -- Built-in Function: X = syl (A, B, C) Solve the Sylvester equation A X + X B = C using standard LAPACK subroutines. For example: sylvester ([1, 2; 3, 4], [5, 6; 7, 8], [9, 10; 11, 12]) => [ 0.50000, 0.66667; 0.66667, 0.50000 ]  File: octave.info, Node: Specialized Solvers, Prev: Functions of a Matrix, Up: Linear Algebra 18.5 Specialized Solvers ======================== -- Function File: X = bicg (A, B, RTOL, MAXIT, M1, M2, X0) -- Function File: X = bicg (A, B, RTOL, MAXIT, P) -- Function File: [X, FLAG, RELRES, ITER, RESVEC] = bicg (A, B, ...) Solve 'A x = b' using the Bi-conjugate gradient iterative method. - RTOL is the relative tolerance, if not given or set to [] the default value 1e-6 is used. - MAXIT the maximum number of outer iterations, if not given or set to [] the default value 'min (20, numel (b))' is used. - X0 the initial guess, if not given or set to [] the default value 'zeros (size (b))' is used. A can be passed as a matrix or as a function handle or inline function 'f' such that 'f(x, "notransp") = A*x' and 'f(x, "transp") = A'*x'. The preconditioner P is given as 'P = M1 * M2'. Both M1 and M2 can be passed as a matrix or as a function handle or inline function 'g' such that 'g(x, "notransp") = M1 \ x' or 'g(x, "notransp") = M2 \ x' and 'g(x, "transp") = M1' \ x' or 'g(x, "transp") = M2' \ x'. If called with more than one output parameter - FLAG indicates the exit status: - 0: iteration converged to the within the chosen tolerance - 1: the maximum number of iterations was reached before convergence - 3: the algorithm reached stagnation (the value 2 is unused but skipped for compatibility). - RELRES is the final value of the relative residual. - ITER is the number of iterations performed. - RESVEC is a vector containing the relative residual at each iteration. See also: *note bicgstab: XREFbicgstab, *note cgs: XREFcgs, *note gmres: XREFgmres, *note pcg: XREFpcg, *note qmr: XREFqmr. -- Function File: X = bicgstab (A, B, RTOL, MAXIT, M1, M2, X0) -- Function File: X = bicgstab (A, B, RTOL, MAXIT, P) -- Function File: [X, FLAG, RELRES, ITER, RESVEC] = bicgstab (A, B, ...) Solve 'A x = b' using the stabilizied Bi-conjugate gradient iterative method. - RTOL is the relative tolerance, if not given or set to [] the default value 1e-6 is used. - MAXIT the maximum number of outer iterations, if not given or set to [] the default value 'min (20, numel (b))' is used. - X0 the initial guess, if not given or set to [] the default value 'zeros (size (b))' is used. A can be passed as a matrix or as a function handle or inline function 'f' such that 'f(x) = A*x'. The preconditioner P is given as 'P = M1 * M2'. Both M1 and M2 can be passed as a matrix or as a function handle or inline function 'g' such that 'g(x) = M1 \ x' or 'g(x) = M2 \ x'. If called with more than one output parameter - FLAG indicates the exit status: - 0: iteration converged to the within the chosen tolerance - 1: the maximum number of iterations was reached before convergence - 3: the algorithm reached stagnation (the value 2 is unused but skipped for compatibility). - RELRES is the final value of the relative residual. - ITER is the number of iterations performed. - RESVEC is a vector containing the relative residual at each iteration. See also: *note bicg: XREFbicg, *note cgs: XREFcgs, *note gmres: XREFgmres, *note pcg: XREFpcg, *note qmr: XREFqmr. -- Function File: X = cgs (A, B, RTOL, MAXIT, M1, M2, X0) -- Function File: X = cgs (A, B, RTOL, MAXIT, P) -- Function File: [X, FLAG, RELRES, ITER, RESVEC] = cgs (A, B, ...) Solve 'A x = b', where A is a square matrix, using the Conjugate Gradients Squared method. - RTOL is the relative tolerance, if not given or set to [] the default value 1e-6 is used. - MAXIT the maximum number of outer iterations, if not given or set to [] the default value 'min (20, numel (b))' is used. - X0 the initial guess, if not given or set to [] the default value 'zeros (size (b))' is used. A can be passed as a matrix or as a function handle or inline function 'f' such that 'f(x) = A*x'. The preconditioner P is given as 'P = M1 * M2'. Both M1 and M2 can be passed as a matrix or as a function handle or inline function 'g' such that 'g(x) = M1 \ x' or 'g(x) = M2 \ x'. If called with more than one output parameter - FLAG indicates the exit status: - 0: iteration converged to the within the chosen tolerance - 1: the maximum number of iterations was reached before convergence - 3: the algorithm reached stagnation (the value 2 is unused but skipped for compatibility). - RELRES is the final value of the relative residual. - ITER is the number of iterations performed. - RESVEC is a vector containing the relative residual at each iteration. See also: *note pcg: XREFpcg, *note bicgstab: XREFbicgstab, *note bicg: XREFbicg, *note gmres: XREFgmres, *note qmr: XREFqmr. -- Function File: X = gmres (A, B, M, RTOL, MAXIT, M1, M2, X0) -- Function File: X = gmres (A, B, M, RTOL, MAXIT, P) -- Function File: [X, FLAG, RELRES, ITER, RESVEC] = gmres (...) Solve 'A x = b' using the Preconditioned GMRES iterative method with restart, a.k.a. PGMRES(m). - RTOL is the relative tolerance, if not given or set to [] the default value 1e-6 is used. - MAXIT is the maximum number of outer iterations, if not given or set to [] the default value 'min (10, numel (b) / restart)' is used. - X0 is the initial guess, if not given or set to [] the default value 'zeros (size (b))' is used. - M is the restart parameter, if not given or set to [] the default value 'numel (b)' is used. Argument A can be passed as a matrix, function handle, or inline function 'f' such that 'f(x) = A*x'. The preconditioner P is given as 'P = M1 * M2'. Both M1 and M2 can be passed as a matrix, function handle, or inline function 'g' such that 'g(x) = M1\x' or 'g(x) = M2\x'. Besides the vector X, additional outputs are: - FLAG indicates the exit status: 0 : iteration converged to within the specified tolerance 1 : maximum number of iterations exceeded 2 : unused, but skipped for compatibility 3 : algorithm reached stagnation (no change between iterations) - RELRES is the final value of the relative residual. - ITER is a vector containing the number of outer iterations and total iterations performed. - RESVEC is a vector containing the relative residual at each iteration. See also: *note bicg: XREFbicg, *note bicgstab: XREFbicgstab, *note cgs: XREFcgs, *note pcg: XREFpcg, *note pcr: XREFpcr, *note qmr: XREFqmr. -- Function File: X = qmr (A, B, RTOL, MAXIT, M1, M2, X0) -- Function File: X = qmr (A, B, RTOL, MAXIT, P) -- Function File: [X, FLAG, RELRES, ITER, RESVEC] = qmr (A, B, ...) Solve 'A x = b' using the Quasi-Minimal Residual iterative method (without look-ahead). - RTOL is the relative tolerance, if not given or set to [] the default value 1e-6 is used. - MAXIT the maximum number of outer iterations, if not given or set to [] the default value 'min (20, numel (b))' is used. - X0 the initial guess, if not given or set to [] the default value 'zeros (size (b))' is used. A can be passed as a matrix or as a function handle or inline function 'f' such that 'f(x, "notransp") = A*x' and 'f(x, "transp") = A'*x'. The preconditioner P is given as 'P = M1 * M2'. Both M1 and M2 can be passed as a matrix or as a function handle or inline function 'g' such that 'g(x, "notransp") = M1 \ x' or 'g(x, "notransp") = M2 \ x' and 'g(x, "transp") = M1' \ x' or 'g(x, "transp") = M2' \ x'. If called with more than one output parameter - FLAG indicates the exit status: - 0: iteration converged to the within the chosen tolerance - 1: the maximum number of iterations was reached before convergence - 3: the algorithm reached stagnation (the value 2 is unused but skipped for compatibility). - RELRES is the final value of the relative residual. - ITER is the number of iterations performed. - RESVEC is a vector containing the residual norms at each iteration. References: 1. R. Freund and N. Nachtigal, 'QMR: a quasi-minimal residual method for non-Hermitian linear systems', Numerische Mathematik, 1991, 60, pp. 315-339. 2. R. Barrett, M. Berry, T. Chan, J. Demmel, J. Donato, J. Dongarra, V. Eijkhour, R. Pozo, C. Romine, and H. van der Vorst, 'Templates for the solution of linear systems: Building blocks for iterative methods', SIAM, 2nd ed., 1994. See also: *note bicg: XREFbicg, *note bicgstab: XREFbicgstab, *note cgs: XREFcgs, *note gmres: XREFgmres, *note pcg: XREFpcg.  File: octave.info, Node: Vectorization and Faster Code Execution, Next: Nonlinear Equations, Prev: Linear Algebra, Up: Top 19 Vectorization and Faster Code Execution ****************************************** Vectorization is a programming technique that uses vector operations instead of element-by-element loop-based operations. Besides frequently producing more succinct Octave code, vectorization also allows for better optimization in the subsequent implementation. The optimizations may occur either in Octave's own Fortran, C, or C++ internal implementation, or even at a lower level depending on the compiler and external numerical libraries used to build Octave. The ultimate goal is to make use of your hardware's vector instructions if possible or to perform other optimizations in software. Vectorization is not a concept unique to Octave, but it is particularly important because Octave is a matrix-oriented language. Vectorized Octave code will see a dramatic speed up (10X-100X) in most cases. This chapter discusses vectorization and other techniques for writing faster code. * Menu: * Basic Vectorization:: Basic techniques for code optimization * Broadcasting:: Broadcasting operations * Function Application:: Applying functions to arrays, cells, and structs * Accumulation:: Accumulation functions * JIT Compiler:: Just-In-Time Compiler for loops * Miscellaneous Techniques:: Other techniques for speeding up code * Examples::  File: octave.info, Node: Basic Vectorization, Next: Broadcasting, Up: Vectorization and Faster Code Execution 19.1 Basic Vectorization ======================== To a very good first approximation, the goal in vectorization is to write code that avoids loops and uses whole-array operations. As a trivial example, consider for i = 1:n for j = 1:m c(i,j) = a(i,j) + b(i,j); endfor endfor compared to the much simpler c = a + b; This isn't merely easier to write; it is also internally much easier to optimize. Octave delegates this operation to an underlying implementation which, among other optimizations, may use special vector hardware instructions or could conceivably even perform the additions in parallel. In general, if the code is vectorized, the underlying implementation has more freedom about the assumptions it can make in order to achieve faster execution. This is especially important for loops with "cheap" bodies. Often it suffices to vectorize just the innermost loop to get acceptable performance. A general rule of thumb is that the "order" of the vectorized body should be greater or equal to the "order" of the enclosing loop. As a less trivial example, instead of for i = 1:n-1 a(i) = b(i+1) - b(i); endfor write a = b(2:n) - b(1:n-1); This shows an important general concept about using arrays for indexing instead of looping over an index variable. *Note Index Expressions::. Also use boolean indexing generously. If a condition needs to be tested, this condition can also be written as a boolean index. For instance, instead of for i = 1:n if (a(i) > 5) a(i) -= 20 endif endfor write a(a>5) -= 20; which exploits the fact that 'a > 5' produces a boolean index. Use elementwise vector operators whenever possible to avoid looping (operators like '.*' and '.^'). *Note Arithmetic Ops::. For simple inline functions, the 'vectorize' function can do this automatically. -- Built-in Function: vectorize (FUN) Create a vectorized version of the inline function FUN by replacing all occurrences of '*', '/', etc., with '.*', './', etc. This may be useful, for example, when using inline functions with numerical integration or optimization where a vector-valued function is expected. fcn = vectorize (inline ("x^2 - 1")) => fcn = f(x) = x.^2 - 1 quadv (fcn, 0, 3) => 6 See also: *note inline: XREFinline, *note formula: XREFformula, *note argnames: XREFargnames. Also exploit broadcasting in these elementwise operators both to avoid looping and unnecessary intermediate memory allocations. *Note Broadcasting::. Use built-in and library functions if possible. Built-in and compiled functions are very fast. Even with an m-file library function, chances are good that it is already optimized, or will be optimized more in a future release. For instance, even better than a = b(2:n) - b(1:n-1); is a = diff (b); Most Octave functions are written with vector and array arguments in mind. If you find yourself writing a loop with a very simple operation, chances are that such a function already exists. The following functions occur frequently in vectorized code: * Index manipulation * find * sub2ind * ind2sub * sort * unique * lookup * ifelse / merge * Repetition * repmat * repelems * Vectorized arithmetic * sum * prod * cumsum * cumprod * sumsq * diff * dot * cummax * cummin * Shape of higher dimensional arrays * reshape * resize * permute * squeeze * deal  File: octave.info, Node: Broadcasting, Next: Function Application, Prev: Basic Vectorization, Up: Vectorization and Faster Code Execution 19.2 Broadcasting ================= Broadcasting refers to how Octave binary operators and functions behave when their matrix or array operands or arguments differ in size. Since version 3.6.0, Octave now automatically broadcasts vectors, matrices, and arrays when using elementwise binary operators and functions. Broadly speaking, smaller arrays are "broadcast" across the larger one, until they have a compatible shape. The rule is that corresponding array dimensions must either 1. be equal, or 2. one of them must be 1. In case all dimensions are equal, no broadcasting occurs and ordinary element-by-element arithmetic takes place. For arrays of higher dimensions, if the number of dimensions isn't the same, then missing trailing dimensions are treated as 1. When one of the dimensions is 1, the array with that singleton dimension gets copied along that dimension until it matches the dimension of the other array. For example, consider x = [1 2 3; 4 5 6; 7 8 9]; y = [10 20 30]; x + y Without broadcasting, 'x + y' would be an error because the dimensions do not agree. However, with broadcasting it is as if the following operation were performed: x = [1 2 3 4 5 6 7 8 9]; y = [10 20 30 10 20 30 10 20 30]; x + y => 11 22 33 14 25 36 17 28 39 That is, the smaller array of size '[1 3]' gets copied along the singleton dimension (the number of rows) until it is '[3 3]'. No actual copying takes place, however. The internal implementation reuses elements along the necessary dimension in order to achieve the desired effect without copying in memory. Both arrays can be broadcast across each other, for example, all pairwise differences of the elements of a vector with itself: y - y' => 0 10 20 -10 0 10 -20 -10 0 Here the vectors of size '[1 3]' and '[3 1]' both get broadcast into matrices of size '[3 3]' before ordinary matrix subtraction takes place. A special case of broadcasting that may be familiar is when all dimensions of the array being broadcast are 1, i.e., the array is a scalar. Thus for example, operations like 'x - 42' and 'max (x, 2)' are basic examples of broadcasting. For a higher-dimensional example, suppose 'img' is an RGB image of size '[m n 3]' and we wish to multiply each color by a different scalar. The following code accomplishes this with broadcasting, img .*= permute ([0.8, 0.9, 1.2], [1, 3, 2]); Note the usage of permute to match the dimensions of the '[0.8, 0.9, 1.2]' vector with 'img'. For functions that are not written with broadcasting semantics, 'bsxfun' can be useful for coercing them to broadcast. -- Built-in Function: bsxfun (F, A, B) The binary singleton expansion function performs broadcasting, that is, it applies a binary function F element-by-element to two array arguments A and B, and expands as necessary singleton dimensions in either input argument. F is a function handle, inline function, or string containing the name of the function to evaluate. The function F must be capable of accepting two column-vector arguments of equal length, or one column vector argument and a scalar. The dimensions of A and B must be equal or singleton. The singleton dimensions of the arrays will be expanded to the same dimensionality as the other array. See also: *note arrayfun: XREFarrayfun, *note cellfun: XREFcellfun. Broadcasting is only applied if either of the two broadcasting conditions hold. As usual, however, broadcasting does not apply when two dimensions differ and neither is 1: x = [1 2 3 4 5 6]; y = [10 20 30 40]; x + y This will produce an error about nonconformant arguments. Besides common arithmetic operations, several functions of two arguments also broadcast. The full list of functions and operators that broadcast is plus + .+ minus - .- times .* rdivide ./ ldivide .\ power .^ .** lt < le <= eq == gt > ge >= ne != ~= and & or | atan2 hypot max min mod rem xor += -= .+= .-= .*= ./= .\= .^= .**= &= |= Beware of resorting to broadcasting if a simpler operation will suffice. For matrices A and B, consider the following: C = sum (permute (A, [1, 3, 2]) .* permute (B, [3, 2, 1]), 3); This operation broadcasts the two matrices with permuted dimensions across each other during elementwise multiplication in order to obtain a larger 3-D array, and this array is then summed along the third dimension. A moment of thought will prove that this operation is simply the much faster ordinary matrix multiplication, 'C = A*B;'. A note on terminology: "broadcasting" is the term popularized by the Numpy numerical environment in the Python programming language. In other programming languages and environments, broadcasting may also be known as _binary singleton expansion_ (BSX, in MATLAB, and the origin of the name of the 'bsxfun' function), _recycling_ (R programming language), _single-instruction multiple data_ (SIMD), or _replication_. 19.2.1 Broadcasting and Legacy Code ----------------------------------- The new broadcasting semantics almost never affect code that worked in previous versions of Octave. Consequently, all code inherited from MATLAB that worked in previous versions of Octave should still work without change in Octave. The only exception is code such as try c = a.*b; catch c = a.*a; end_try_catch that may have relied on matrices of different size producing an error. Due to how broadcasting changes semantics with older versions of Octave, by default Octave warns if a broadcasting operation is performed. To disable this warning, refer to its ID (*note warning_ids: XREFwarning_ids.): warning ("off", "Octave:broadcast"); If you want to recover the old behavior and produce an error, turn this warning into an error: warning ("error", "Octave:broadcast"); For broadcasting on scalars that worked in previous versions of Octave, this warning will not be emitted.  File: octave.info, Node: Function Application, Next: Accumulation, Prev: Broadcasting, Up: Vectorization and Faster Code Execution 19.3 Function Application ========================= As a general rule, functions should already be written with matrix arguments in mind and should consider whole matrix operations in a vectorized manner. Sometimes, writing functions in this way appears difficult or impossible for various reasons. For those situations, Octave provides facilities for applying a function to each element of an array, cell, or struct. -- Function File: arrayfun (FUNC, A) -- Function File: X = arrayfun (FUNC, A) -- Function File: X = arrayfun (FUNC, A, B, ...) -- Function File: [X, Y, ...] = arrayfun (FUNC, A, ...) -- Function File: arrayfun (..., "UniformOutput", VAL) -- Function File: arrayfun (..., "ErrorHandler", ERRFUNC) Execute a function on each element of an array. This is useful for functions that do not accept array arguments. If the function does accept array arguments it is better to call the function directly. The first input argument FUNC can be a string, a function handle, an inline function, or an anonymous function. The input argument A can be a logic array, a numeric array, a string array, a structure array, or a cell array. By a call of the function 'arrayfun' all elements of A are passed on to the named function FUNC individually. The named function can also take more than two input arguments, with the input arguments given as third input argument B, fourth input argument C, ... If given more than one array input argument then all input arguments must have the same sizes, for example: arrayfun (@atan2, [1, 0], [0, 1]) => [ 1.5708 0.0000 ] If the parameter VAL after a further string input argument "UniformOutput" is set 'true' (the default), then the named function FUNC must return a single element which then will be concatenated into the return value and is of type matrix. Otherwise, if that parameter is set to 'false', then the outputs are concatenated in a cell array. For example: arrayfun (@(x,y) x:y, "abc", "def", "UniformOutput", false) => { [1,1] = abcd [1,2] = bcde [1,3] = cdef } If more than one output arguments are given then the named function must return the number of return values that also are expected, for example: [A, B, C] = arrayfun (@find, [10; 0], "UniformOutput", false) => A = { [1,1] = 1 [2,1] = [](0x0) } B = { [1,1] = 1 [2,1] = [](0x0) } C = { [1,1] = 10 [2,1] = [](0x0) } If the parameter ERRFUNC after a further string input argument "ErrorHandler" is another string, a function handle, an inline function, or an anonymous function, then ERRFUNC defines a function to call in the case that FUNC generates an error. The definition of the function must be of the form function [...] = errfunc (S, ...) where there is an additional input argument to ERRFUNC relative to FUNC, given by S. This is a structure with the elements "identifier", "message", and "index" giving, respectively, the error identifier, the error message, and the index of the array elements that caused the error. The size of the output argument of ERRFUNC must have the same size as the output argument of FUNC, otherwise a real error is thrown. For example: function y = ferr (s, x), y = "MyString"; endfunction arrayfun (@str2num, [1234], "UniformOutput", false, "ErrorHandler", @ferr) => { [1,1] = MyString } See also: *note spfun: XREFspfun, *note cellfun: XREFcellfun, *note structfun: XREFstructfun. -- Function File: Y = spfun (F, S) Compute 'f(S)' for the nonzero values of S. This results in a sparse matrix with the same structure as S. The function F can be passed as a string, a function handle, or an inline function. See also: *note arrayfun: XREFarrayfun, *note cellfun: XREFcellfun, *note structfun: XREFstructfun. -- Built-in Function: cellfun (NAME, C) -- Built-in Function: cellfun ("size", C, K) -- Built-in Function: cellfun ("isclass", C, CLASS) -- Built-in Function: cellfun (FUNC, C) -- Built-in Function: cellfun (FUNC, C, D) -- Built-in Function: [A, ...] = cellfun (...) -- Built-in Function: cellfun (..., "ErrorHandler", ERRFUNC) -- Built-in Function: cellfun (..., "UniformOutput", VAL) Evaluate the function named NAME on the elements of the cell array C. Elements in C are passed on to the named function individually. The function NAME can be one of the functions 'isempty' Return 1 for empty elements. 'islogical' Return 1 for logical elements. 'isnumeric' Return 1 for numeric elements. 'isreal' Return 1 for real elements. 'length' Return a vector of the lengths of cell elements. 'ndims' Return the number of dimensions of each element. 'numel' 'prodofsize' Return the number of elements contained within each cell element. The number is the product of the dimensions of the object at each cell element. 'size' Return the size along the K-th dimension. 'isclass' Return 1 for elements of CLASS. Additionally, 'cellfun' accepts an arbitrary function FUNC in the form of an inline function, function handle, or the name of a function (in a character string). The function can take one or more arguments, with the inputs arguments given by C, D, etc. Equally the function can return one or more output arguments. For example: cellfun ("atan2", {1, 0}, {0, 1}) => [ 1.57080 0.00000 ] The number of output arguments of 'cellfun' matches the number of output arguments of the function. The outputs of the function will be collected into the output arguments of 'cellfun' like this: function [a, b] = twoouts (x) a = x; b = x*x; endfunction [aa, bb] = cellfun (@twoouts, {1, 2, 3}) => aa = 1 2 3 bb = 1 4 9 Note that per default the output argument(s) are arrays of the same size as the input arguments. Input arguments that are singleton (1x1) cells will be automatically expanded to the size of the other arguments. If the parameter "UniformOutput" is set to true (the default), then the function must return scalars which will be concatenated into the return array(s). If "UniformOutput" is false, the outputs are concatenated into a cell array (or cell arrays). For example: cellfun ("tolower", {"Foo", "Bar", "FooBar"}, "UniformOutput", false) => {"foo", "bar", "foobar"} Given the parameter "ErrorHandler", then ERRFUNC defines a function to call in case FUNC generates an error. The form of the function is function [...] = errfunc (S, ...) where there is an additional input argument to ERRFUNC relative to FUNC, given by S. This is a structure with the elements "identifier", "message" and "index", giving respectively the error identifier, the error message, and the index into the input arguments of the element that caused the error. For example: function y = foo (s, x), y = NaN; endfunction cellfun ("factorial", {-1,2}, "ErrorHandler", @foo) => [NaN 2] Use 'cellfun' intelligently. The 'cellfun' function is a useful tool for avoiding loops. It is often used with anonymous function handles; however, calling an anonymous function involves an overhead quite comparable to the overhead of an m-file function. Passing a handle to a built-in function is faster, because the interpreter is not involved in the internal loop. For example: a = {...} v = cellfun (@(x) det (x), a); # compute determinants v = cellfun (@det, a); # faster See also: *note arrayfun: XREFarrayfun, *note structfun: XREFstructfun, *note spfun: XREFspfun. -- Function File: structfun (FUNC, S) -- Function File: [A, ...] = structfun (...) -- Function File: structfun (..., "ErrorHandler", ERRFUNC) -- Function File: structfun (..., "UniformOutput", VAL) Evaluate the function named NAME on the fields of the structure S. The fields of S are passed to the function FUNC individually. 'structfun' accepts an arbitrary function FUNC in the form of an inline function, function handle, or the name of a function (in a character string). In the case of a character string argument, the function must accept a single argument named X, and it must return a string value. If the function returns more than one argument, they are returned as separate output variables. If the parameter "UniformOutput" is set to true (the default), then the function must return a single element which will be concatenated into the return value. If "UniformOutput" is false, the outputs are placed into a structure with the same fieldnames as the input structure. s.name1 = "John Smith"; s.name2 = "Jill Jones"; structfun (@(x) regexp (x, '(\w+)$', "matches"){1}, s, "UniformOutput", false) => { name1 = Smith name2 = Jones } Given the parameter "ErrorHandler", ERRFUNC defines a function to call in case FUNC generates an error. The form of the function is function [...] = errfunc (SE, ...) where there is an additional input argument to ERRFUNC relative to FUNC, given by SE. This is a structure with the elements "identifier", "message" and "index", giving respectively the error identifier, the error message, and the index into the input arguments of the element that caused the error. For an example on how to use an error handler, *note cellfun: XREFcellfun. See also: *note cellfun: XREFcellfun, *note arrayfun: XREFarrayfun, *note spfun: XREFspfun. Consistent with earlier advice, seek to use Octave built-in functions whenever possible for the best performance. This advice applies especially to the four functions above. For example, when adding two arrays together element-by-element one could use a handle to the built-in addition function '@plus' or define an anonymous function '@(x,y) x + y'. But, the anonymous function is 60% slower than the first method. *Note Operator Overloading::, for a list of basic functions which might be used in place of anonymous ones.  File: octave.info, Node: Accumulation, Next: JIT Compiler, Prev: Function Application, Up: Vectorization and Faster Code Execution 19.4 Accumulation ================= Whenever it's possible to categorize according to indices the elements of an array when performing a computation, accumulation functions can be useful. -- Function File: accumarray (SUBS, VALS, SZ, FUNC, FILLVAL, ISSPARSE) -- Function File: accumarray (SUBS, VALS, ...) Create an array by accumulating the elements of a vector into the positions defined by their subscripts. The subscripts are defined by the rows of the matrix SUBS and the values by VALS. Each row of SUBS corresponds to one of the values in VALS. If VALS is a scalar, it will be used for each of the row of SUBS. If SUBS is a cell array of vectors, all vectors must be of the same length, and the subscripts in the Kth vector must correspond to the Kth dimension of the result. The size of the matrix will be determined by the subscripts themselves. However, if SZ is defined it determines the matrix size. The length of SZ must correspond to the number of columns in SUBS. An exception is if SUBS has only one column, in which case SZ may be the dimensions of a vector and the subscripts of SUBS are taken as the indices into it. The default action of 'accumarray' is to sum the elements with the same subscripts. This behavior can be modified by defining the FUNC function. This should be a function or function handle that accepts a column vector and returns a scalar. The result of the function should not depend on the order of the subscripts. The elements of the returned array that have no subscripts associated with them are set to zero. Defining FILLVAL to some other value allows these values to be defined. This behavior changes, however, for certain values of FUNC. If FUNC is 'min' (respectively, 'max') then the result will be filled with the minimum (respectively, maximum) integer if VALS is of integral type, logical false (respectively, logical true) if VALS is of logical type, zero if FILLVAL is zero and all values are non-positive (respectively, non-negative), and NaN otherwise. By default 'accumarray' returns a full matrix. If ISSPARSE is logically true, then a sparse matrix is returned instead. The following 'accumarray' example constructs a frequency table that in the first column counts how many occurrences each number in the second column has, taken from the vector X. Note the usage of 'unique' for assigning to all repeated elements of X the same index (*note unique: XREFunique.). X = [91, 92, 90, 92, 90, 89, 91, 89, 90, 100, 100, 100]; [U, ~, J] = unique (X); [accumarray(J', 1), U'] => 2 89 3 90 2 91 2 92 3 100 Another example, where the result is a multi-dimensional 3-D array and the default value (zero) appears in the output: accumarray ([1, 1, 1; 2, 1, 2; 2, 3, 2; 2, 1, 2; 2, 3, 2], 101:105) => ans(:,:,1) = [101, 0, 0; 0, 0, 0] => ans(:,:,2) = [0, 0, 0; 206, 0, 208] The sparse option can be used as an alternative to the 'sparse' constructor (*note sparse: XREFsparse.). Thus sparse (I, J, SV) can be written with 'accumarray' as accumarray ([I, J], SV', [], [], 0, true) For repeated indices, 'sparse' adds the corresponding value. To take the minimum instead, use 'min' as an accumulator function: accumarray ([I, J], SV', [], @min, 0, true) The complexity of accumarray in general for the non-sparse case is generally O(M+N), where N is the number of subscripts and M is the maximum subscript (linearized in multi-dimensional case). If FUNC is one of '@sum' (default), '@max', '@min' or '@(x) {x}', an optimized code path is used. Note that for general reduction function the interpreter overhead can play a major part and it may be more efficient to do multiple accumarray calls and compute the results in a vectorized manner. See also: *note accumdim: XREFaccumdim, *note unique: XREFunique, *note sparse: XREFsparse. -- Function File: accumdim (SUBS, VALS, DIM, N, FUNC, FILLVAL) Create an array by accumulating the slices of an array into the positions defined by their subscripts along a specified dimension. The subscripts are defined by the index vector SUBS. The dimension is specified by DIM. If not given, it defaults to the first non-singleton dimension. The length of SUBS must be equal to 'size (VALS, DIM)'. The extent of the result matrix in the working dimension will be determined by the subscripts themselves. However, if N is defined it determines this extent. The default action of 'accumdim' is to sum the subarrays with the same subscripts. This behavior can be modified by defining the FUNC function. This should be a function or function handle that accepts an array and a dimension, and reduces the array along this dimension. As a special exception, the built-in 'min' and 'max' functions can be used directly, and 'accumdim' accounts for the middle empty argument that is used in their calling. The slices of the returned array that have no subscripts associated with them are set to zero. Defining FILLVAL to some other value allows these values to be defined. An example of the use of 'accumdim' is: accumdim ([1, 2, 1, 2, 1], [ 7, -10, 4; -5, -12, 8; -12, 2, 8; -10, 9, -3; -5, -3, -13]) => [-10,-11,-1;-15,-3,5] See also: *note accumarray: XREFaccumarray.  File: octave.info, Node: JIT Compiler, Next: Miscellaneous Techniques, Prev: Accumulation, Up: Vectorization and Faster Code Execution 19.5 JIT Compiler ================= Vectorization is the preferred technique for eliminating loops and speeding up code. Nevertheless, it is not always possible to replace every loop. In such situations it may be worth trying Octave's *experimental* Just-In-Time (JIT) compiler. A JIT compiler works by analyzing the body of a loop, translating the Octave statements into another language, compiling the new code segment into an executable, and then running the executable and collecting any results. The process is not simple and there is a significant amount of work to perform for each step. It can still make sense, however, if the number of loop iterations is large. Because Octave is an interpreted language every time through a loop Octave must parse the statements in the loop body before executing them. With a JIT compiler this is done just once when the body is translated to another language. The JIT compiler is a very new feature in Octave and not all valid Octave statements can currently be accelerated. However, if no other technique is available it may be worth benchmarking the code with JIT enabled. The function 'jit_enable' is used to turn compilation on or off. The function 'jit_startcnt' sets the threshold for acceleration. Loops with iteration counts above 'jit_startcnt' will be accelerated. The functions 'jit_failcnt' and 'debug_jit' are not likely to be of use to anyone not working directly on the implementation of the JIT compiler. -- Built-in Function: VAL = jit_enable () -- Built-in Function: OLD_VAL = jit_enable (NEW_VAL) -- Built-in Function: jit_enable (NEW_VAL, "local") Query or set the internal variable that enables Octave's JIT compiler. When called from inside a function with the "local" option, the variable is changed locally for the function and any subroutines it calls. The original variable value is restored when exiting the function. See also: *note jit_startcnt: XREFjit_startcnt, *note debug_jit: XREFdebug_jit. -- Built-in Function: VAL = jit_startcnt () -- Built-in Function: OLD_VAL = jit_startcnt (NEW_VAL) -- Built-in Function: jit_startcnt (NEW_VAL, "local") Query or set the internal variable that determines whether JIT compilation will take place for a specific loop. Because compilation is a costly operation it does not make sense to employ JIT when the loop count is low. By default only loops with greater than 1000 iterations will be accelerated. When called from inside a function with the "local" option, the variable is changed locally for the function and any subroutines it calls. The original variable value is restored when exiting the function. See also: *note jit_enable: XREFjit_enable, *note jit_failcnt: XREFjit_failcnt, *note debug_jit: XREFdebug_jit. -- Built-in Function: VAL = jit_failcnt () -- Built-in Function: OLD_VAL = jit_failcnt (NEW_VAL) -- Built-in Function: jit_failcnt (NEW_VAL, "local") Query or set the internal variable that counts the number of JIT fail exceptions for Octave's JIT compiler. When called from inside a function with the "local" option, the variable is changed locally for the function and any subroutines it calls. The original variable value is restored when exiting the function. See also: *note jit_enable: XREFjit_enable, *note jit_startcnt: XREFjit_startcnt, *note debug_jit: XREFdebug_jit. -- Built-in Function: VAL = debug_jit () -- Built-in Function: OLD_VAL = debug_jit (NEW_VAL) -- Built-in Function: debug_jit (NEW_VAL, "local") Query or set the internal variable that determines whether debugging/tracing is enabled for Octave's JIT compiler. When called from inside a function with the "local" option, the variable is changed locally for the function and any subroutines it calls. The original variable value is restored when exiting the function. See also: *note jit_enable: XREFjit_enable, *note jit_startcnt: XREFjit_startcnt.  File: octave.info, Node: Miscellaneous Techniques, Next: Examples, Prev: JIT Compiler, Up: Vectorization and Faster Code Execution 19.6 Miscellaneous Techniques ============================= Here are some other ways of improving the execution speed of Octave programs. * Avoid computing costly intermediate results multiple times. Octave currently does not eliminate common subexpressions. Also, certain internal computation results are cached for variables. For instance, if a matrix variable is used multiple times as an index, checking the indices (and internal conversion to integers) is only done once. * Be aware of lazy copies (copy-on-write). When a copy of an object is created, the data is not immediately copied, but rather shared. The actual copying is postponed until the copied data needs to be modified. For example: a = zeros (1000); # create a 1000x1000 matrix b = a; # no copying done here b(1) = 1; # copying done here Lazy copying applies to whole Octave objects such as matrices, cells, struct, and also individual cell or struct elements (not array elements). Additionally, index expressions also use lazy copying when Octave can determine that the indexed portion is contiguous in memory. For example: a = zeros (1000); # create a 1000x1000 matrix b = a(:,10:100); # no copying done here b = a(10:100,:); # copying done here This applies to arrays (matrices), cell arrays, and structs indexed using '()'. Index expressions generating comma-separated lists can also benefit from shallow copying in some cases. In particular, when A is a struct array, expressions like '{a.x}, {a(:,2).x}' will use lazy copying, so that data can be shared between a struct array and a cell array. Most indexing expressions do not live longer than their parent objects. In rare cases, however, a lazily copied slice outlasts its parent, in which case it becomes orphaned, still occupying unnecessarily more memory than needed. To provide a remedy working in most real cases, Octave checks for orphaned lazy slices at certain situations, when a value is stored into a "permanent" location, such as a named variable or cell or struct element, and possibly economizes them. For example: a = zeros (1000); # create a 1000x1000 matrix b = a(:,10:100); # lazy slice a = []; # the original "a" array is still allocated c{1} = b; # b is reallocated at this point * Avoid deep recursion. Function calls to m-file functions carry a relatively significant overhead, so rewriting a recursion as a loop often helps. Also, note that the maximum level of recursion is limited. * Avoid resizing matrices unnecessarily. When building a single result matrix from a series of calculations, set the size of the result matrix first, then insert values into it. Write result = zeros (big_n, big_m) for i = over:and_over ridx = ... cidx = ... result(ridx, cidx) = new_value (); endfor instead of result = []; for i = ever:and_ever result = [ result, new_value() ]; endfor Sometimes the number of items can not be computed in advance, and stack-like operations are needed. When elements are being repeatedly inserted or removed from the end of an array, Octave detects it as stack usage and attempts to use a smarter memory management strategy by pre-allocating the array in bigger chunks. This strategy is also applied to cell and struct arrays. a = []; while (condition) ... a(end+1) = value; # "push" operation ... a(end) = []; # "pop" operation ... endwhile * Avoid calling 'eval' or 'feval' excessively. Parsing input or looking up the name of a function in the symbol table are relatively expensive operations. If you are using 'eval' merely as an exception handling mechanism, and not because you need to execute some arbitrary text, use the 'try' statement instead. *Note The try Statement::. * Use 'ignore_function_time_stamp' when appropriate. If you are calling lots of functions, and none of them will need to change during your run, set the variable 'ignore_function_time_stamp' to "all". This will stop Octave from checking the time stamp of a function file to see if it has been updated while the program is being run.  File: octave.info, Node: Examples, Prev: Miscellaneous Techniques, Up: Vectorization and Faster Code Execution 19.7 Examples ============= The following are examples of vectorization questions asked by actual users of Octave and their solutions. * For a vector 'A', the following loop n = length (A); B = zeros (n, 2); for i = 1:length (A) ## this will be two columns, the first is the difference and ## the second the mean of the two elements used for the diff. B(i,:) = [A(i+1)-A(i), (A(i+1) + A(i))/2]; endfor can be turned into the following one-liner: B = [diff(A)(:), 0.5*(A(1:end-1)+A(2:end))(:)] Note the usage of colon indexing to flatten an intermediate result into a column vector. This is a common vectorization trick.  File: octave.info, Node: Nonlinear Equations, Next: Diagonal and Permutation Matrices, Prev: Vectorization and Faster Code Execution, Up: Top 20 Nonlinear Equations ********************** * Menu: * Solvers:: * Minimizers::  File: octave.info, Node: Solvers, Next: Minimizers, Up: Nonlinear Equations 20.1 Solvers ============ Octave can solve sets of nonlinear equations of the form F (x) = 0 using the function 'fsolve', which is based on the MINPACK subroutine 'hybrd'. This is an iterative technique so a starting point must be provided. This also has the consequence that convergence is not guaranteed even if a solution exists. -- Function File: fsolve (FCN, X0, OPTIONS) -- Function File: [X, FVEC, INFO, OUTPUT, FJAC] = fsolve (FCN, ...) Solve a system of nonlinear equations defined by the function FCN. FCN should accept a vector (array) defining the unknown variables, and return a vector of left-hand sides of the equations. Right-hand sides are defined to be zeros. In other words, this function attempts to determine a vector X such that 'FCN (X)' gives (approximately) all zeros. X0 determines a starting guess. The shape of X0 is preserved in all calls to FCN, but otherwise it is treated as a column vector. OPTIONS is a structure specifying additional options. Currently, 'fsolve' recognizes these options: "FunValCheck", "OutputFcn", "TolX", "TolFun", "MaxIter", "MaxFunEvals", "Jacobian", "Updating", "ComplexEqn" "TypicalX", "AutoScaling" and "FinDiffType". If "Jacobian" is "on", it specifies that FCN, called with 2 output arguments also returns the Jacobian matrix of right-hand sides at the requested point. "TolX" specifies the termination tolerance in the unknown variables, while "TolFun" is a tolerance for equations. Default is '1e-7' for both "TolX" and "TolFun". If "AutoScaling" is on, the variables will be automatically scaled according to the column norms of the (estimated) Jacobian. As a result, TolF becomes scaling-independent. By default, this option is off because it may sometimes deliver unexpected (though mathematically correct) results. If "Updating" is "on", the function will attempt to use Broyden updates to update the Jacobian, in order to reduce the amount of Jacobian calculations. If your user function always calculates the Jacobian (regardless of number of output arguments) then this option provides no advantage and should be set to false. "ComplexEqn" is "on", 'fsolve' will attempt to solve complex equations in complex variables, assuming that the equations possess a complex derivative (i.e., are holomorphic). If this is not what you want, you should unpack the real and imaginary parts of the system to get a real system. For description of the other options, see 'optimset'. On return, FVAL contains the value of the function FCN evaluated at X. INFO may be one of the following values: 1 Converged to a solution point. Relative residual error is less than specified by TolFun. 2 Last relative step size was less that TolX. 3 Last relative decrease in residual was less than TolF. 0 Iteration limit exceeded. -3 The trust region radius became excessively small. Note: If you only have a single nonlinear equation of one variable, using 'fzero' is usually a much better idea. Note about user-supplied Jacobians: As an inherent property of the algorithm, a Jacobian is always requested for a solution vector whose residual vector is already known, and it is the last accepted successful step. Often this will be one of the last two calls, but not always. If the savings by reusing intermediate results from residual calculation in Jacobian calculation are significant, the best strategy is to employ OutputFcn: After a vector is evaluated for residuals, if OutputFcn is called with that vector, then the intermediate results should be saved for future Jacobian evaluation, and should be kept until a Jacobian evaluation is requested or until OutputFcn is called with a different vector, in which case they should be dropped in favor of this most recent vector. A short example how this can be achieved follows: function [fvec, fjac] = user_func (x, optimvalues, state) persistent sav = [], sav0 = []; if (nargin == 1) ## evaluation call if (nargout == 1) sav0.x = x; # mark saved vector ## calculate fvec, save results to sav0. elseif (nargout == 2) ## calculate fjac using sav. endif else ## outputfcn call. if (all (x == sav0.x)) sav = sav0; endif ## maybe output iteration status, etc. endif endfunction ## ... fsolve (@user_func, x0, optimset ("OutputFcn", @user_func, ...)) See also: *note fzero: XREFfzero, *note optimset: XREFoptimset. The following is a complete example. To solve the set of equations -2x^2 + 3xy + 4 sin(y) = 6 3x^2 - 2xy^2 + 3 cos(x) = -4 you first need to write a function to compute the value of the given function. For example: function y = f (x) y = zeros (2, 1); y(1) = -2*x(1)^2 + 3*x(1)*x(2) + 4*sin(x(2)) - 6; y(2) = 3*x(1)^2 - 2*x(1)*x(2)^2 + 3*cos(x(1)) + 4; endfunction Then, call 'fsolve' with a specified initial condition to find the roots of the system of equations. For example, given the function 'f' defined above, [x, fval, info] = fsolve (@f, [1; 2]) results in the solution x = 0.57983 2.54621 fval = -5.7184e-10 5.5460e-10 info = 1 A value of 'info = 1' indicates that the solution has converged. When no Jacobian is supplied (as in the example above) it is approximated numerically. This requires more function evaluations, and hence is less efficient. In the example above we could compute the Jacobian analytically as function [y, jac] = f (x) y = zeros (2, 1); y(1) = -2*x(1)^2 + 3*x(1)*x(2) + 4*sin(x(2)) - 6; y(2) = 3*x(1)^2 - 2*x(1)*x(2)^2 + 3*cos(x(1)) + 4; if (nargout == 2) jac = zeros (2, 2); jac(1,1) = 3*x(2) - 4*x(1); jac(1,2) = 4*cos(x(2)) + 3*x(1); jac(2,1) = -2*x(2)^2 - 3*sin(x(1)) + 6*x(1); jac(2,2) = -4*x(1)*x(2); endif endfunction The Jacobian can then be used with the following call to 'fsolve': [x, fval, info] = fsolve (@f, [1; 2], optimset ("jacobian", "on")); which gives the same solution as before. -- Function File: fzero (FUN, X0) -- Function File: fzero (FUN, X0, OPTIONS) -- Function File: [X, FVAL, INFO, OUTPUT] = fzero (...) Find a zero of a univariate function. FUN is a function handle, inline function, or string containing the name of the function to evaluate. X0 should be a two-element vector specifying two points which bracket a zero. In other words, there must be a change in sign of the function between X0(1) and X0(2). More mathematically, the following must hold sign (FUN(X0(1))) * sign (FUN(X0(2))) <= 0 If X0 is a single scalar then several nearby and distant values are probed in an attempt to obtain a valid bracketing. If this is not successful, the function fails. OPTIONS is a structure specifying additional options. Currently, 'fzero' recognizes these options: "FunValCheck", "OutputFcn", "TolX", "MaxIter", "MaxFunEvals". For a description of these options, see *note optimset: XREFoptimset. On exit, the function returns X, the approximate zero point and FVAL, the function value thereof. INFO is an exit flag that can have these values: * 1 The algorithm converged to a solution. * 0 Maximum number of iterations or function evaluations has been reached. * -1 The algorithm has been terminated from user output function. * -5 The algorithm may have converged to a singular point. OUTPUT is a structure containing runtime information about the 'fzero' algorithm. Fields in the structure are: * iterations Number of iterations through loop. * nfev Number of function evaluations. * bracketx A two-element vector with the final bracketing of the zero along the x-axis. * brackety A two-element vector with the final bracketing of the zero along the y-axis. See also: *note optimset: XREFoptimset, *note fsolve: XREFfsolve.  File: octave.info, Node: Minimizers, Prev: Solvers, Up: Nonlinear Equations 20.2 Minimizers =============== Often it is useful to find the minimum value of a function rather than just the zeroes where it crosses the x-axis. 'fminbnd' is designed for the simpler, but very common, case of a univariate function where the interval to search is bounded. For unbounded minimization of a function with potentially many variables use 'fminunc' or 'fminsearch'. The two functions use different internal algorithms and some knowledge of the objective function is required. For functions which can be differentiated, 'fminunc' is appropriate. For functions with discontinuities, or for which a gradient search would fail, use 'fminsearch'. *Note Optimization::, for minimization with the presence of constraint functions. Note that searches can be made for maxima by simply inverting the objective function ('Fto_max = -Fto_min'). -- Function File: [X, FVAL, INFO, OUTPUT] = fminbnd (FUN, A, B, OPTIONS) Find a minimum point of a univariate function. FUN should be a function handle or name. A, B specify a starting interval. OPTIONS is a structure specifying additional options. Currently, 'fminbnd' recognizes these options: "FunValCheck", "OutputFcn", "TolX", "MaxIter", "MaxFunEvals". For a description of these options, see *note optimset: XREFoptimset. On exit, the function returns X, the approximate minimum point and FVAL, the function value thereof. INFO is an exit flag that can have these values: * 1 The algorithm converged to a solution. * 0 Maximum number of iterations or function evaluations has been exhausted. * -1 The algorithm has been terminated from user output function. Notes: The search for a minimum is restricted to be in the interval bound by A and B. If you only have an initial point to begin searching from you will need to use an unconstrained minimization algorithm such as 'fminunc' or 'fminsearch'. 'fminbnd' internally uses a Golden Section search strategy. See also: *note fzero: XREFfzero, *note fminunc: XREFfminunc, *note fminsearch: XREFfminsearch, *note optimset: XREFoptimset. -- Function File: fminunc (FCN, X0) -- Function File: fminunc (FCN, X0, OPTIONS) -- Function File: [X, FVAL, INFO, OUTPUT, GRAD, HESS] = fminunc (FCN, ...) Solve an unconstrained optimization problem defined by the function FCN. FCN should accept a vector (array) defining the unknown variables, and return the objective function value, optionally with gradient. 'fminunc' attempts to determine a vector X such that 'FCN (X)' is a local minimum. X0 determines a starting guess. The shape of X0 is preserved in all calls to FCN, but otherwise is treated as a column vector. OPTIONS is a structure specifying additional options. Currently, 'fminunc' recognizes these options: "FunValCheck", "OutputFcn", "TolX", "TolFun", "MaxIter", "MaxFunEvals", "GradObj", "FinDiffType", "TypicalX", "AutoScaling". If "GradObj" is "on", it specifies that FCN, when called with 2 output arguments, also returns the Jacobian matrix of partial first derivatives at the requested point. 'TolX' specifies the termination tolerance for the unknown variables X, while 'TolFun' is a tolerance for the objective function value FVAL. The default is '1e-7' for both options. For a description of the other options, see 'optimset'. On return, X is the location of the minimum and FVAL contains the value of the objective function at X. INFO may be one of the following values: 1 Converged to a solution point. Relative gradient error is less than specified by 'TolFun'. 2 Last relative step size was less than 'TolX'. 3 Last relative change in function value was less than 'TolFun'. 0 Iteration limit exceeded--either maximum number of algorithm iterations 'MaxIter' or maximum number of function evaluations 'MaxFunEvals'. -1 Algorithm terminated by 'OutputFcn'. -3 The trust region radius became excessively small. Optionally, 'fminunc' can return a structure with convergence statistics (OUTPUT), the output gradient (GRAD) at the solution X, and approximate Hessian (HESS) at the solution X. Application Notes: If have only a single nonlinear equation of one variable then using 'fminbnd' is usually a better choice. The algorithm used by 'fminsearch' is a gradient search which depends on the objective function being differentiable. If the function has discontinuities it may be better to use a derivative-free algorithm such as 'fminsearch'. See also: *note fminbnd: XREFfminbnd, *note fminsearch: XREFfminsearch, *note optimset: XREFoptimset. -- Function File: X = fminsearch (FUN, X0) -- Function File: X = fminsearch (FUN, X0, OPTIONS) -- Function File: [X, FVAL] = fminsearch (...) Find a value of X which minimizes the function FUN. The search begins at the point X0 and iterates using the Nelder & Mead Simplex algorithm (a derivative-free method). This algorithm is better-suited to functions which have discontinuities or for which a gradient-based search such as 'fminunc' fails. Options for the search are provided in the parameter OPTIONS using the function 'optimset'. Currently, 'fminsearch' accepts the options: "TolX", "MaxFunEvals", "MaxIter", "Display". For a description of these options, see 'optimset'. On exit, the function returns X, the minimum point, and FVAL, the function value thereof. Example usages: fminsearch (@(x) (x(1)-5).^2+(x(2)-8).^4, [0;0]) fminsearch (inline ("(x(1)-5).^2+(x(2)-8).^4", "x"), [0;0]) See also: *note fminbnd: XREFfminbnd, *note fminunc: XREFfminunc, *note optimset: XREFoptimset.  File: octave.info, Node: Diagonal and Permutation Matrices, Next: Sparse Matrices, Prev: Nonlinear Equations, Up: Top 21 Diagonal and Permutation Matrices ************************************ * Menu: * Basic Usage:: Creation and Manipulation of Diagonal/Permutation Matrices * Matrix Algebra:: Linear Algebra with Diagonal/Permutation Matrices * Function Support:: Functions That Are Aware of These Matrices * Example Code:: Examples of Usage * Zeros Treatment:: Differences in Treatment of Zero Elements  File: octave.info, Node: Basic Usage, Next: Matrix Algebra, Up: Diagonal and Permutation Matrices 21.1 Creating and Manipulating Diagonal/Permutation Matrices ============================================================ A diagonal matrix is defined as a matrix that has zero entries outside the main diagonal; that is, 'D(i,j) == 0' if 'i != j'. Most often, square diagonal matrices are considered; however, the definition can equally be applied to non-square matrices, in which case we usually speak of a rectangular diagonal matrix. A permutation matrix is defined as a square matrix that has a single element equal to unity in each row and each column; all other elements are zero. That is, there exists a permutation (vector) 'p' such that 'P(i,j) == 1' if 'j == p(i)' and 'P(i,j) == 0' otherwise. Octave provides special treatment of real and complex rectangular diagonal matrices, as well as permutation matrices. They are stored as special objects, using efficient storage and algorithms, facilitating writing both readable and efficient matrix algebra expressions in the Octave language. The special treatment may be disabled by using the functions "disable_diagonal_matrix" and "disable_permutation_matrix". -- Built-in Function: VAL = disable_diagonal_matrix () -- Built-in Function: OLD_VAL = disable_diagonal_matrix (NEW_VAL) -- Built-in Function: disable_diagonal_matrix (NEW_VAL, "local") Query or set the internal variable that controls whether diagonal matrices are stored in a special space-efficient format. The default value is true. If this option is disabled Octave will store diagonal matrices as full matrices. When called from inside a function with the "local" option, the variable is changed locally for the function and any subroutines it calls. The original variable value is restored when exiting the function. See also: *note disable_range: XREFdisable_range, *note disable_permutation_matrix: XREFdisable_permutation_matrix. -- Built-in Function: VAL = disable_permutation_matrix () -- Built-in Function: OLD_VAL = disable_permutation_matrix (NEW_VAL) -- Built-in Function: disable_permutation_matrix (NEW_VAL, "local") Query or set the internal variable that controls whether permutation matrices are stored in a special space-efficient format. The default value is true. If this option is disabled Octave will store permutation matrices as full matrices. When called from inside a function with the "local" option, the variable is changed locally for the function and any subroutines it calls. The original variable value is restored when exiting the function. See also: *note disable_range: XREFdisable_range, *note disable_diagonal_matrix: XREFdisable_diagonal_matrix. The space savings are significant as demonstrated by the following code. x = diag (rand (10, 1)); xf = full (x); sizeof (x) => 80 sizeof (xf) => 800 * Menu: * Creating Diagonal Matrices:: * Creating Permutation Matrices:: * Explicit and Implicit Conversions::  File: octave.info, Node: Creating Diagonal Matrices, Next: Creating Permutation Matrices, Up: Basic Usage 21.1.1 Creating Diagonal Matrices --------------------------------- The most common and easiest way to create a diagonal matrix is using the built-in function "diag". The expression 'diag (v)', with V a vector, will create a square diagonal matrix with elements on the main diagonal given by the elements of V, and size equal to the length of V. 'diag (v, m, n)' can be used to construct a rectangular diagonal matrix. The result of these expressions will be a special diagonal matrix object, rather than a general matrix object. Diagonal matrix with unit elements can be created using "eye". Some other built-in functions can also return diagonal matrices. Examples include "balance" or "inv". Example: diag (1:4) => Diagonal Matrix 1 0 0 0 0 2 0 0 0 0 3 0 0 0 0 4 diag (1:3,5,3) => Diagonal Matrix 1 0 0 0 2 0 0 0 3 0 0 0 0 0 0  File: octave.info, Node: Creating Permutation Matrices, Next: Explicit and Implicit Conversions, Prev: Creating Diagonal Matrices, Up: Basic Usage 21.1.2 Creating Permutation Matrices ------------------------------------ For creating permutation matrices, Octave does not introduce a new function, but rather overrides an existing syntax: permutation matrices can be conveniently created by indexing an identity matrix by permutation vectors. That is, if Q is a permutation vector of length N, the expression P = eye (n) (:, q); will create a permutation matrix - a special matrix object. eye (n) (q, :) will also work (and create a row permutation matrix), as well as eye (n) (q1, q2). For example: eye (4) ([1,3,2,4],:) => Permutation Matrix 1 0 0 0 0 0 1 0 0 1 0 0 0 0 0 1 eye (4) (:,[1,3,2,4]) => Permutation Matrix 1 0 0 0 0 0 1 0 0 1 0 0 0 0 0 1 Mathematically, an identity matrix is both diagonal and permutation matrix. In Octave, 'eye (n)' returns a diagonal matrix, because a matrix can only have one class. You can convert this diagonal matrix to a permutation matrix by indexing it by an identity permutation, as shown below. This is a special property of the identity matrix; indexing other diagonal matrices generally produces a full matrix. eye (3) => Diagonal Matrix 1 0 0 0 1 0 0 0 1 eye(3)(1:3,:) => Permutation Matrix 1 0 0 0 1 0 0 0 1 Some other built-in functions can also return permutation matrices. Examples include "inv" or "lu".