Thursday, August 13, 2009

gmres


Generalized Minimum Residual method (with restarts)

Syntax

  • x = gmres(A,b)
    gmres(A,b,restart)
    gmres(A,b,restart,tol)
    gmres(A,b,restart,tol,maxit)
    gmres(A,b,restart,tol,maxit,M)
    gmres(A,b,restart,tol,maxit,M1,M2)
    gmres(A,b,restart,tol,maxit,M1,M2,x0)
    gmres(afun,b,restart,tol,maxit,m1fun,m2fun,x0,p1,p2,...)
    [x,flag] = gmres(A,b,...)
    [x,flag,relres] = gmres(A,b,...)
    [x,flag,relres,iter] = gmres(A,b,...)
    [x,flag,relres,iter,resvec] = gmres(A,b,...)

Description

x = gmres(A,b) attempts to solve the system of linear equations A*x = b for x. The n-by-n coefficient matrix A must be square and should be large and sparse. The column vector b must have length n. A can be a function afun such that afun(x) returns A*x. For this syntax, gmres does not restart; the maximum number of iterations is min(n,10).

If gmres converges, a message to that effect is displayed. If gmres fails to converge after the maximum number of iterations or halts for any reason, a warning message is printed displaying the relative residual norm(b-A*x)/norm(b) and the iteration number at which the method stopped or failed.

gmres(A,b,restart) restarts the method every restart inner iterations. The maximum number of outer iterations is min(n/restart,10). The maximum number of total iterations is restart*min(n/restart,10). If restart is n or [], then gmres does not restart and the maximum number of total iterations is min(n,10).

gmres(A,b,restart,tol) specifies the tolerance of the method. If tol is [], then gmres uses the default, 1e-6.

gmres(A,b,restart,tol,maxit) specifies the maximum number of outer iterations, i.e., the total number of iterations does not exceed restart*maxit. If maxit is [] then gmres uses the default, min(n/restart,10). If restart is n or [], then the maximum number of total iterations is maxit (instead of restart*maxit).

gmres(A,b,restart,tol,maxit,M) and gmres(A,b,restart,tol,maxit,M1,M2) use preconditioner M or M = M1*M2 and effectively solve the system inv(M)*A*x = inv(M)*b for x. If M is [] then gmres applies no preconditioner. M can be a function that returns M\x.

gmres(A,b,restart,tol,maxit,M1,M2,x0) specifies the first initial guess. If x0 is [], then gmres uses the default, an all-zero vector.

gmres(afun,b,restart,tol,maxit,m1fun,m2fun,x0,p1,p2,...) passes parameters to functions afun(x,p1,p2,...), m1fun(x,p1,p2,...), and m2fun(x,p1,p2,...).

[x,flag] = gmres(A,b,...) also returns a convergence flag:




flag = 0
gmres converged to the desired tolerance tol within maxit outer iterations.
flag = 1
gmres iterated maxit times but did not converge.
flag = 2
Preconditioner M was ill-conditioned.
flag = 3
gmres stagnated. (Two consecutive iterates were the same.)

Whenever flag is not 0, the solution x returned is that with minimal norm residual computed over all the iterations. No messages are displayed if the flag output is specified.

[x,flag,relres] = gmres(A,b,...) also returns the relative residual norm(b-A*x)/norm(b). If flag is 0, relres <= tol.

[x,flag,relres,iter] = gmres(A,b,...) also returns both the outer and inner iteration numbers at which x was computed, where 0 <= iter(1) <= maxit and 0 <= iter(2) <= restart.

[x,flag,relres,iter,resvec] = gmres(A,b,...) also returns a vector of the residual norms at each inner iteration, including norm(b-A*x0).

Examples

Example 1.

  • A = gallery('wilk',21); 
    b = sum(A,2);
    tol = 1e-12;
    maxit = 15;
    M1 = diag([10:-1:1 1 1:10]);

    x = gmres(A,b,10,tol,maxit,M1,[],[]);
    gmres(10) converged at iteration 2(10) to a solution with relative
    residual 1.9e-013

Alternatively, use this matrix-vector product function

  • function y = afun(x,n)
    y = [0;
    x(1:n-1)] + [((n-1)/2:-1:0)';
    (1:(n-1)/2)'] .*x + [x(2:n);
    0];

and this preconditioner backsolve function

  • function y = mfun(r,n)
    y = r ./ [((n-1)/2:-1:1)'; 1; (1:(n-1)/2)'];

as inputs to gmres

  • x1 = gmres(@afun,b,10,tol,maxit,@mfun,[],[],21);

Note that both afun and mfun must accept the gmres extra input n=21.

Example 2.

  • load west0479
    A = west0479
    b = sum(A,2)
    [x,flag] = gmres(A,b,5)

flag is 1 because gmres does not converge to the default tolerance 1e-6 within the default 10 outer iterations.

  • [L1,U1] = luinc(A,1e-5);
    [x1,flag1] = gmres(A,b,5,1e-6,5,L1,U1);

flag1 is 2 because the upper triangular U1 has a zero on its diagonal, and gmres fails in the first iteration when it tries to solve a system such as U1*y = r for y using backslash.

  • [L2,U2] = luinc(A,1e-6);
    tol = 1e-15;
    [x4,flag4,relres4,iter4,resvec4] = gmres(A,b,4,tol,5,L2,U2);
    [x6,flag6,relres6,iter6,resvec6] = gmres(A,b,6,tol,3,L2,U2);
    [x8,flag8,relres8,iter8,resvec8] = gmres(A,b,8,tol,3,L2,U2);

flag4, flag6, and flag8 are all 0 because gmres converged when restarted at iterations 4, 6, and 8 while preconditioned by the incomplete LU factorization with a drop tolerance of 1e-6. This is verified by the plots of outer iteration number against relative residual. A combined plot of all three clearly shows the restarting at iterations 4 and 6. The total number of iterations computed may be more for lower values of restart, but the number of length n vectors stored is fewer, and the amount of work done in the method decreases proportionally.

global


Define a global variable

Syntax

  • global X Y Z

Description

global X Y Z defines X, Y, and Z as global in scope.

Ordinarily, each MATLAB function, defined by an M-file, has its own local variables, which are separate from those of other functions, and from those of the base workspace. However, if several functions, and possibly the base workspace, all declare a particular name as global, they all share a single copy of that variable. Any assignment to that variable, in any function, is available to all the functions declaring it global.

If the global variable does not exist the first time you issue the global statement, it is initialized to the empty matrix.

If a variable with the same name as the global variable already exists in the current workspace, MATLAB issues a warning and changes the value of that variable to match the global.

Remarks

Use clear global variable to clear a global variable from the global workspace. Use clear variable to clear the global link from the current workspace without affecting the value of the global.

To use a global within a callback, declare the global, use it, then clear the global link from the workspace. This avoids declaring the global after it has been referenced. For example,

uicontrol('style','pushbutton','CallBack',...

'global MY_GLOBAL,disp(MY_GLOBAL),MY_GLOBAL = MY_GLOBAL+1,clear MY_GLOBAL',...

'string','count')

There is no function form of the global command (i.e., you cannot use parentheses and quote the variable names).

Examples

Here is the code for the functions tic and toc (some comments abridged). These functions manipulate a stopwatch-like timer. The global variable TICTOC is shared by the two functions, but it is invisible in the base workspace or in any other functions that do not declare it.

  • function tic
    % TIC Start a stopwatch timer.
    % TIC; any stuff; TOC
    % prints the time required.
    % See also: TOC, CLOCK.
    global TICTOC
    TICTOC = clock;

    function t = toc
    % TOC Read the stopwatch timer.
    % TOC prints the elapsed time since TIC was used.
    % t = TOC; saves elapsed time in t, does not print.
    % See also: TIC, ETIME.
    global TICTOC
    if nargout < 1
    elapsed_time = etime(clock,TICTOC)
    else
    t = etime(clock,TICTOC);
    end
ginput


Input data using the mouse

Syntax

  • [x,y] = ginput(n)
    [x,y] = ginput
    [x,y,button] = ginput(...)

Description

ginput enables you to select points from the figure using the mouse or arrow keys for cursor positioning. The figure must have focus before ginput receives input.

[x,y] = ginput(n) enables you to select n points from the current axes and returns the x- and y-coordinates in the column vectors x and y, respectively. You can press the Return key to terminate the input before entering n points.

[x,y] = ginput gathers an unlimited number of points until you press the Return key.

[x,y,button] = ginput(...) returns the x-coordinates, the y-coordinates, and the button or key designation. button is a vector of integers indicating which mouse buttons you pressed (1 for left, 2 for middle, 3 for right), or ASCII numbers indicating which keys on the keyboard you pressed.

Remarks

If you select points from multiple axes, the results you get are relative to those axes coordinates systems.

Examples

Pick 10 two-dimensional points from the figure window.

  • [x,y] = ginput(10)

Position the cursor with the mouse (or the arrow keys on terminals without a mouse, such as Tektronix emulators). Enter data points by pressing a mouse button or a key on the keyboard. To terminate input before entering 10 points, press the Return key.

ginput


Input data using the mouse

Syntax

  • [x,y] = ginput(n)
    [x,y] = ginput
    [x,y,button] = ginput(...)

Description

ginput enables you to select points from the figure using the mouse or arrow keys for cursor positioning. The figure must have focus before ginput receives input.

[x,y] = ginput(n) enables you to select n points from the current axes and returns the x- and y-coordinates in the column vectors x and y, respectively. You can press the Return key to terminate the input before entering n points.

[x,y] = ginput gathers an unlimited number of points until you press the Return key.

[x,y,button] = ginput(...) returns the x-coordinates, the y-coordinates, and the button or key designation. button is a vector of integers indicating which mouse buttons you pressed (1 for left, 2 for middle, 3 for right), or ASCII numbers indicating which keys on the keyboard you pressed.

Remarks

If you select points from multiple axes, the results you get are relative to those axes coordinates systems.

Examples

Pick 10 two-dimensional points from the figure window.

  • [x,y] = ginput(10)

Position the cursor with the mouse (or the arrow keys on terminals without a mouse, such as Tektronix emulators). Enter data points by pressing a mouse button or a key on the keyboard. To terminate input before entering 10 points, press the Return key.

get (timer)


Display or get timer object properties

Syntax

  • get(obj)
    out = get(obj)
    out = get(obj,'PropertyName')

Description

get(obj) displays all property names and their current values for timer object obj.

V = get(obj) returns a structure, V, where each field name is the name of a property of obj and each field contains the value of that property.

V = get(obj,'PropertyName') returns the value, V, of the timer object property specified in PropertyName.

If PropertyName is a1-by-N or N-by-1 cell array of strings containing property names, V is a 1-by-N cell array of values. If obj is a vector of timer objects, V is an M-by-N cell array of property values where M is equal to the length of obj and N is equal to the number of properties specified.

Example

  • t = timer;
    get(t)
    AveragePeriod: NaN
    BusyMode: 'drop'
    ErrorFcn: []
    ExecutionMode: 'singleShot'
    InstantPeriod: NaN
    LastError: 'none'
    Name: 'timer-1'
    Period: 1
    Running: 'off'
    StartDelay: 0
    StartFcn: []
    StopFcn: []
    Tag: ''
    TasksToExecute: Inf
    TasksExecuted: 0
    TimerFcn: []
    Type: 'timer'
    UserData: []
    get(t, {'StartDelay','Period'})
    ans =

    [0] [1]


getappdata

Get value of application-defined data

Syntax

  • value = getappdata(h,name)
    values = getappdata(h)

Description

value = getappdata(h,name) gets the value of the application-defined data with the name specified by name, in the object with the handle h. If the application-defined data does not exist, MATLAB returns an empty matrix in value.

value = getappdata(h) returns all application-defined data for the object with handle h.




getenv

Get environment variable

Syntax

  • getenv 'name'
    N = getenv('name')

Description

getenv 'name' searches the underlying operating system's environment list for a string of the form name=value, where name is the input string. If found, MATLAB returns the string, value. If the specified name cannot be found, an empty matrix is returned.

N = getenv('name') returns value to the variable, N.

Examples

  • os = getenv('OS')

    os =
    Windows_NT


getfield

Get field of structure array

    Note getfield is obsolete and will be removed in a future release. Please use dynamic field names instead.

Syntax

  • f = getfield(s,'field')
    f = getfield(s,{i,j},'field',{k})

Description

f = getfield(s,'field'), where s is a 1-by-1 structure, returns the contents of the specified field. This is equivalent to the syntax f = s.field.

If s is a structure having dimensions greater than 1-by-1, getfield returns the first of all output values requested in the call. That is, for structure array s(m,n), getfield returns f = s(1,1).field.

f = getfield(s,{i,j},'field',{k}) returns the contents of the specified field. This is equivalent to the syntax f = s(i,j).field(k). All subscripts must be passed as cell arrays--that is, they must be enclosed in curly braces (similar to{i,j} and {k} above). Pass field references as strings.

Examples

Given the structure

  • mystr(1,1).name = 'alice';
    mystr(1,1).ID = 0;
    mystr(2,1).name = 'gertrude';
    mystr(2,1).ID = 1

Then the command f = getfield(mystr,{2,1},'name') yields

  • f =

    gertrude

To list the contents of all name (or other) fields, embed getfield in a loop.

  • for k = 1:2
    name{k} = getfield(mystr,{k,1},'name');
    end
    name

    name =

    'alice' 'gertrude'

The following example starts out by creating a structure using the standard structure syntax. It then reads the fields of the structure using getfield with variable and quoted field names and additional subscripting arguments.

  • class = 5;     student = 'John_Doe';
    grades(class).John_Doe.Math(10,21:30) = ...
    [85, 89, 76, 93, 85, 91, 68, 84, 95, 73];

Use getfield to access the structure fields.

  • getfield(grades,{class}, student, 'Math', {10,21:30})

    ans =

    85 89 76 93 85 91 68 84 95 73




getframe


Get movie frame

Syntax

  • F = getframe
    F = getframe(h)
    F = getframe(h,rect)
    [X,Map] = getframe(...)

Description

getframe returns a movie frame. The frame is a snapshot (pixmap) of the current axes or figure.

F = getframe gets a frame from the current axes.

F = getframe(h) gets a frame from the figure or axes identified by the handle h.

F = getframe(h,rect) specifies a rectangular area from which to copy the pixmap. rect is relative to the lower-left corner of the figure or axes h, in pixel units. rect is a four-element vector in the form [left bottom width height], where width and height define the dimensions of the rectangle.

F = getframe(...) returns a movie frame, which is a structure having two fields:

  • cdata - The image data stored as a matrix of uint8 values. The dimensions of F.cdata are height-by-width-by-3.
  • colormap - The colormap stored as an n-by-3 matrix of doubles. F.colormap is empty on true color systems.

To capture an image, use this approach:

  • F = getframe(gcf);
    image(F.cdata)
    colormap(F.colormap)

[X,Map] = getframe(...) returns the frame as an indexed image matrix X and a colormap Map. This version is obsolete and is supported only for compatibility with earlier version of MATLAB. Since indexed images cannot always capture true color displays, you should use the single output argument form of getframe. To write code that is compatible with earlier version of MATLAB and that can take advantage of true color support, use the following approach:

  • F = getframe(gcf);
    [X,Map] = frame2im(f);
    imshow(X,Map)

Remarks

Usually, getframe is used in a for loop to assemble an array of movie frames for playback using movie. For example,

  • for j = 1:n
    plotting commands
    F(j) = getframe;
    end
    movie(F)

To create movies that are compatible with earlier versions of MATLAB (before Release 11/MATLAB 5.3) use this approach:

  • M = moviein(n);
    for j = 1:n
    plotting commands
    M(:,j) = getframe;
    end
    movie(M)

Capture Regions

Note that F = getframe; returns the contents of the current axes, exclusive of the axis labels, title, or tick labels. F = getframe(gcf); captures the entire interior of the current figure window. To capture the figure window menu, use the form F = getframe(h,rect) with a rectangle sized to include the menu.

Examples

Make the peaks function vibrate.

  • Z = peaks; surf(Z)
    axis tight
    set(gca,'nextplot','replacechildren');
    for j = 1:20
    surf(sin(2*pi*j/20)*Z,Z)
    F(j) = getframe;
    end
    movie(F,20) % Play the movie twenty times

get

Get object properties

Syntax

  • get(h)
    get(h,'PropertyName')
    = get(H,)
    a = get(h)
    a = get(0,'Factory')
    a = get(0,'FactoryObjectTypePropertyName')
    a = get(h,'Default')
    a = get(h,'DefaultObjectTypePropertyName')

Description

get(h) returns all properties and their current values of the graphics object identified by the handle h.

get(h,'PropertyName') returns the value of the property 'PropertyName' of the graphics object identified by h.

= get(H,pn) returns n property values for m graphics objects in the m-by-n cell array, where m = length(H) and n is equal to the number of property names contained in pn.

a = get(h) returns a structure whose field names are the object's property names and whose values are the current values of the corresponding properties. h must be a scalar. If you do not specify an output argument, MATLAB displays the information on the screen.

a = get(0,'Factory') returns the factory-defined values of all user-settable properties. a is a structure array whose field names are the object property names and whose field values are the values of the corresponding properties. If you do not specify an output argument, MATLAB displays the information on the screen.

a = get(0,'FactoryObjectTypePropertyName') returns the factory-defined value of the named property for the specified object type. The argument, FactoryObjectTypePropertyName, is the word Factory concatenated with the object type (e.g., Figure) and the property name (e.g., Color).

  • FactoryFigureColor

a = get(h,'Default') returns all default values currently defined on object h. a is a structure array whose field names are the object property names and whose field values are the values of the corresponding properties. If you do not specify an output argument, MATLAB displays the information on the screen.

a = get(h,'DefaultObjectTypePropertyName') returns the factory-defined value of the named property for the specified object type. The argument, DefaultObjectTypePropertyName, is the word Default concatenated with the object type (e.g., Figure) and the property name (e.g., Color).

  • DefaultFigureColor

Examples

You can obtain the default value of the LineWidth property for line graphics objects defined on the root level with the statement:

  • get(0,'DefaultLineLineWidth')
    ans =
    0.5000

To query a set of properties on all axes children define a cell array of property names:

  • props = {'HandleVisibility', 'Interruptible';
    'SelectionHighlight', 'Type'};
    output = get(get(gca,'Children'),props);

The variable output is a cell array of dimension
length(get(gca,'Children')-by-4.

For example, type

  • patch;surface;text;line
    output = get(get(gca,'Children'),props)
    output =
    'on' 'on' 'on' 'line'
    'on' 'off' 'on' 'text'
    'on' 'on' 'on' 'surface'
    'on' 'on' 'on' 'patch'


get (COM)



Retrieve a property value from an interface or get a list of properties

Syntax

  • v = get(h[, 'propertyname'])

Arguments

h
Handle for a COM object previously returned from actxcontrol, actxserver, get, or invoke.

propertyname
A string that is the name of the property value to be retrieved.

Description

Returns the value of the property specified by propertyname. If no property is specified, then get returns a list of all properties for the object or interface.

The meaning and type of the return value is dependent upon the specific property being retrieved. The object's documentation should describe the specific meaning of the return value.

Examples

Create a COM server running Microsoft Excel:

  • e = actxserver ('Excel.Application');

Retrieve a single property value:

  • get(e, 'Path')
    ans =
    D:\Applications\MSOffice\Office

Retrieve a list of all properties for the CommandBars interface:

  • c = get(e, 'CommandBars');
    get(c)
    ans =
    Application: [1x1
    Interface.excel.application.CommandBars.Application]
    Creator: 1.4808e+009
    ActionControl: []
    ActiveMenuBar: [1x1
    Interface.excel.application.CommandBars.ActiveMenuBar]
    Count: 94
    DisplayTooltips: 1
    DisplayKeysInTooltips: 0
    LargeButtons: 0
    MenuAnimationStyle: 'msoMenuAnimationNone'
    Parent: [1x1
    Interface.excel.application.CommandBars.Parent]
    AdaptiveMenus: 0
    DisplayFonts: 1


get (serial)


Return serial port object properties

Syntax

  • get(obj)
    out = get(obj)
    out = get(obj,'PropertyName')

Arguments


obj
A serial port object or an array of serial port objects.
'PropertyName'
A property name or a cell array of property names.
out
A single property value, a structure of property values, or a cell array of property values.

Description

get(obj) returns all property names and their current values to the command line for obj.

out = get(obj) returns the structure out where each field name is the name of a property of obj, and each field contains the value of that property.

out = get(obj,'PropertyName') returns the value out of the property specified by PropertyName for obj. If PropertyName is replaced by a 1-by-n or n-by-1 cell array of strings containing property names, then get returns a 1-by-n cell array of values to out. If obj is an array of serial port objects, then out will be a m-by-n cell array of property values where m is equal to the length of obj and n is equal to the number of properties specified.

Remarks

Refer to Displaying Property Names and Property Values for a list of serial port object properties that you can return with get.

When you specify a property name, you can do so without regard to case, and you can make use of property name completion. For example, if s is a serial port object, then these commands are all valid.

  • out = get(s,'BaudRate');
    out = get(s,'baudrate');
    out = get(s,'BAUD');

If you use the help command to display help for get, then you need to supply the pathname shown below.

  • help serial/get

Example

This example illustrates some of the ways you can use get to return property values for the serial port object s.

  • s = serial('COM1');
    out1 = get(s);
    out2 = get(s,{'BaudRate','DataBits'});
    get(s,'Parity')
    ans =
    none
genpath


Generate a path string

Syntax

  • genpath
    genpath directory
    p = genpath('directory')

Description

genpath returns a path string formed by recursively adding all the directories below matlabroot/toolbox. Empty directories are not included.

genpath directory returns a path string formed by recursively adding all the directories below directory. Empty directories are not included.

p = genpath('directory') returns the path string to variable, p.

Examples

You generate a path that includes matlabroot/toolbox/images and all directories below that with the following command:

p = genpath(fullfile(matlabroot,'toolbox','images'))

p =

  • matlabroot\toolbox\images;matlabroot\toolbox\images\images; 
    matlabroot\toolbox\images\images\ja;matlabroot\toolbox\images\
    imdemos;matlabroot\toolbox\images\imdemos\ja;

You can also use genpath in conjunction with addpath to add subdirectories to the path from the command line. The following example adds the /control directory and its subdirectories to the current path.

  • %  Display the current path
    path

    MATLABPATH

    K:\toolbox\matlab\general
    K:\toolbox\matlab\ops
    K:\toolbox\matlab\lang
    K:\toolbox\matlab\elmat
    K:\toolbox\matlab\elfun
    :
    :
    :

    % Use GENPATH to add /control and its subdirectories
    addpath(genpath('K:/toolbox/control'))


    % Display the new path
    path

    MATLABPATH

    K:\toolbox\control
    K:\toolbox\control\ctrlutil
    K:\toolbox\control\control
    K:\toolbox\control\ctrlguis
    K:\toolbox\control\ctrldemos
    K:\toolbox\matlab\general
    K:\toolbox\matlab\ops
    K:\toolbox\matlab\lang
    K:\toolbox\matlab\elmat
    K:\toolbox\matlab\elfun
    :
    :
    :



gcbf


Get handle of figure containing object whose callback is executing

Syntax

  • fig = gcbf

Description

fig = gcbf returns the handle of the figure that contains the object whose callback is currently executing. This object can be the figure itself, in which case, gcbf returns the figure's handle.

When no callback is executing, gcbf returns the empty matrix, [].

The value returned by gcbf is identical to the figure output argument returned by gcbo.


gcbo

Return the handle of the object whose callback is currently executing

Syntax

h = gcbo

[h, figure] = gcbo

Description

h = gcbo returns the handle of the graphics object whose callback is executing.

[h, figure] = gcbo returns the handle of the current callback object and the handle of the figure containing this object.

Remarks

MATLAB stores the handle of the object whose callback is executing in the root CallbackObject property. If a callback interrupts another callback, MATLAB replaces the CallbackObject value with the handle of the object whose callback is interrupting. When that callback completes, MATLAB restores the handle of the object whose callback was interrupted.

The root CallbackObject property is read-only, so its value is always valid at any time during callback execution. The root CurrentFigure property, and the figure CurrentAxes and CurrentObject properties (returned by gcf, gca, and gco respectively) are user settable, so they can change during the execution of a callback, especially if that callback is interrupted by another callback. Therefore, those functions are not reliable indicators of which object's callback is executing.

When you write callback routines for the CreateFcn and DeleteFcn of any object and the figure ResizeFcn, you must use gcbo since those callbacks do not update the root's CurrentFigure property, or the figure's CurrentObject or CurrentAxes properties; they only update the root's CallbackObject property.

When no callbacks are executing, gcbo returns [] (an empty matrix).




gcd


Greatest common divisor

Syntax

  • G = gcd(A,B)
    [G,C,D] = gcd(A,B)

Description

G = gcd(A,B) returns an array containing the greatest common divisors of the corresponding elements of integer arrays A and B. By convention, gcd(0,0) returns a value of 0; all other inputs return positive integers for G.

[G,C,D] = gcd(A,B) returns both the greatest common divisor array G, and the arrays C and D, which satisfy the equation: A(i).*C(i) + B(i).*D(i) = G(i). These are useful for solving Diophantine equations and computing elementary Hermite transformations.

Examples

The first example involves elementary Hermite transformations.

For any two integers a and b there is a 2-by-2 matrix E with integer entries and determinant = 1 (a unimodular matrix) such that:

  • E * [a;b] = [g,0],

where g is the greatest common divisor of a and b as returned by the command
[g,c,d] = gcd(a,b).

The matrix E equals:

  • c       d
    -b/g a/g

In the case where a = 2 and b = 4:

  • [g,c,d] = gcd(2,4)
    g =
    2
    c =
    1
    d =
    0

So that

  • E =
    1 0
    -2 1

In the next example, we solve for x and y in the Diophantine equation 30x + 56y = 8.

  • [g,c,d] = gcd(30,56)
    g =
    2
    c =
    -13
    d =
    7

By the definition, for scalars c and d:

  • 30(-13) + 56(7) = 2,

Multiplying through by 8/2:

  • 30(-13*4) + 56(7*4) = 8

Comparing this to the original equation, a solution can be read by inspection:

  • x = (-13*4) = -52; y = (7*4) = 28

gcf



Get current figure handle

Syntax

  • h = gcf

Description

h = gcf returns the handle of the current figure. The current figure is the figure window in which graphics commands such as plot, title, and surf draw their results. If no figure exists, MATLAB creates one and returns its handle. You can use the statement

  • get(0,'CurrentFigure')

if you do not want MATLAB to create a figure if one does not already exist.



gco



Return handle of current object

Syntax

  • h = gco
    h = gco(figure_handle)

Description

h = gco returns the handle of the current object.

h = gco(figure_handle) returns the value of the current object for the figure specified by figure_handle.

Remarks

The current object is the last object clicked on, excluding uimenus. If the mouse click did not occur over a figure child object, the figure becomes the current object. MATLAB stores the handle of the current object in the figure's CurrentObject property.

The CurrentObject of the CurrentFigure does not always indicate the object whose callback is being executed. Interruptions of callbacks by other callbacks can change the CurrentObject or even the CurrentFigure. Some callbacks, such as CreateFcn and DeleteFcn, and uimenu Callback intentionally do not update CurrentFigure or CurrentObject.

gcbo provides the only completely reliable way to retrieve the handle to the object whose callback is executing, at any point in the callback function, regardless of the type of callback or of any previous interruptions.

Examples

This statement returns the handle to the current object in figure window 2:

  • h = gco(2)




axes

gca

Get current axes handle

Syntax

  • h = gca

Description

h = gca returns the handle to the current axes for the current figure. If no axes exists, MATLAB creates one and returns its handle. You can use the statement

  • get(gcf,'CurrentAxes')

if you do not want MATLAB to create an axes if one does not already exist.

The current axes is the target for graphics output when you create axes children. Graphics commands such as plot, text, and surf draw their results in the current axes. Changing the current figure also changes the current axes.



gama function

gamma, gammainc, gammaln

Gamma functions

Syntax

  • Y = gamma(A)                Gamma function
    Y = gammainc(X,A) Incomplete gamma function
    Y = gammaln(A) Logarithm of gamma function

Definition

The gamma function is defined by the integral:

The gamma function interpolates the factorial function. For integer n:

  • gamma(n+1) = n! = prod(1:n)

The incomplete gamma function is:



Description

Y = gamma(A) returns the gamma function at the elements of A. A must be real.

Y = gammainc(X,A) returns the incomplete gamma function of corresponding elements of X and A. Arguments X and A must be real and the same size (or either can be scalar).

Y = gammaln(A) returns the logarithm of the gamma function, gammaln(A) = log(gamma(A)). The gammaln command avoids the underflow and overflow that may occur if it is computed directly using log(gamma(A)).


matrices

gallery

Test matrices

Syntax

  • [A,B,C,...] = gallery('tmfun',P1,P2,...)
    gallery(3) a badly conditioned 3-by-3 matrix
    gallery(5) an interesting eigenvalue problem

Description

[A,B,C,...] = gallery('tmfun',P1,P2,...) returns the test matrices specified by string tmfun. tmfun is the name of a matrix family selected from the table below. P1, P2,... are input parameters required by the individual matrix family. The number of optional parameters P1,P2,... used in the calling syntax varies from matrix to matrix.The exact calling syntaxes are detailed in the individual matrix descriptions below.

The gallery holds over fifty different test matrix functions useful for testing algorithms and other purposes.


Test Matrices
cauchy
chebspec
chebvand
chow
circul
clement
compar
condex
cycol
dorr
dramadah
fiedler
forsythe
frank
gearmat
grcar
hanowa
house
invhess
invol
ipjfact
jordbloc
kahan
kms
krylov
lauchli
lehmer
leslie
lesp
lotkin
minij
moler
neumann
orthog
parter
pei
poisson
prolate
randcolu
randcorr
rando
randhess
randsvd
redheff
riemann
ris
rosser
smoke
toeppd
tridiag
triw
vander
wathen
wilk



cauchy--Cauchy matrix

C = gallery('cauchy',x,y) returns an n-by-n matrix, C(i,j) = 1/(x(i)+y(j)). Arguments x and y are vectors of length n. If you pass in scalars for x and y, they are interpreted as vectors 1:x and 1:y.

C = gallery('cauchy',x) returns the same as above with y = x. That is, the command returns C(i,j) = 1/(x(i)+x(j)).

Explicit formulas are known for the inverse and determinant of a Cauchy matrix. The determinant det(C) is nonzero if x and y both have distinct elements. C is totally positive if 0 <> and 0 <>.

chebspec--Chebyshev spectral differentiation matrix

C = gallery('chebspec',n,switch) returns a Chebyshev spectral differentiation matrix of order n. Argument switch is a variable that determines the character of the output matrix. By default, switch = 0.

For switch = 0 ("no boundary conditions"), C is nilpotent (C^n = 0) and has the null vector ones(n,1). The matrix C is similar to a Jordan block of size n with eigenvalue zero.

For switch = 1, C is nonsingular and well-conditioned, and its eigenvalues have negative real parts.

The eigenvector matrix of the Chebyshev spectral differentiation matrix is ill-conditioned.

chebvand--Vandermonde-like matrix for the Chebyshev polynomials

C = gallery('chebvand',p) produces the (primal) Chebyshev Vandermonde matrix based on the vector of points p, which define where the Chebyshev polynomial is calculated.

C = gallery('chebvand',m,p) where m is scalar, produces a rectangular version of the above, with m rows.

If p is a vector, then where is the Chebyshev polynomial of degree i-1. If p is a scalar, then p equally spaced points on the interval [0,1] are used to calculate C.

chow--Singular Toeplitz lower Hessenberg matrix

A = gallery('chow',n,alpha,delta) returns A such that A = H(alpha) + delta*eye(n), where and argument n is the order of the Chow matrix. Default value for scalars alpha and delta are 1 and 0, respectively.

H(alpha) has p = floor(n/2) eigenvalues that are equal to zero. The rest of the eigenvalues are equal to 4*alpha*cos(k*pi/(n+2))^2, k=1:n-p.

circul--Circulant matrix

C = gallery('circul',v) returns the circulant matrix whose first row is the vector v.

A circulant matrix has the property that each row is obtained from the previous one by cyclically permuting the entries one step forward. It is a special Toeplitz matrix in which the diagonals "wrap around."

If v is a scalar, then C = gallery('circul',1:v).

The eigensystem of C (n-by-n) is known explicitly: If t is an nth root of unity, then the inner product of v and is an eigenvalue of C and w(n:-1:1) is an eigenvector.

clement--Tridiagonal matrix with zero diagonal entries

A = gallery('clement',n,sym) returns an n-by-n tridiagonal matrix with zeros on its main diagonal and known eigenvalues. It is singular if order n is odd. About 64 percent of the entries of the inverse are zero. The eigenvalues include plus and minus the numbers n-1, n-3, n-5, ..., as well as (for odd n) a final eigenvalue of 1 or 0.

Argument sym determines whether the Clement matrix is symmetric. For sym = 0 (the default) the matrix is nonsymmetric, while for sym = 1, it is symmetric.

compar--Comparison matrices

A = gallery('compar',A,1) returns A with each diagonal element replaced by its absolute value, and each off-diagonal element replaced by minus the absolute value of the largest element in absolute value in its row. However, if A is triangular compar(A,1) is too.

gallery('compar',A) is diag(B) - tril(B,-1) - triu(B,1), where B = abs(A). compar(A) is often denoted by M(A) in the literature.

gallery('compar',A,0) is the same as gallery('compar',A).

condex--Counter-examples to matrix condition number estimators

A = gallery('condex',n,k,theta) returns a "counter-example" matrix to a condition estimator. It has order n and scalar parameter theta (default 100).

The matrix, its natural size, and the estimator to which it applies are specified by k:


k = 1
4-by-4
LINPACK
k = 2
3-by-3
LINPACK
k = 3
arbitrary
LINPACK (rcond) (independent of theta)
k = 4
n >= 4
LAPACK (RCOND) (default). It is the inverse of this matrix that is a counter-example.

If n is not equal to the natural size of the matrix, then the matrix is padded out with an identity matrix to order n.

cycol--Matrix whose columns repeat cyclically

A = gallery('cycol',[m n],k) returns an m-by-n matrix with cyclically repeating columns, where one "cycle" consists of randn(m,k). Thus, the rank of matrix A cannot exceed k, and k must be a scalar.

Argument k defaults to round(n/4), and need not evenly divide n.

A = gallery('cycol',n,k), where n is a scalar, is the same as gallery('cycol',[n n],k).

dorr--Diagonally dominant, ill-conditioned, tridiagonal matrix

[c,d,e] = gallery('dorr',n,theta) returns the vectors defining an n-by-n, row diagonally dominant, tridiagonal matrix that is ill-conditioned for small nonnegative values of theta. The default value of theta is 0.01. The Dorr matrix itself is the same as gallery('tridiag',c,d,e).

A = gallery('dorr',n,theta) returns the matrix itself, rather than the defining vectors.

dramadah--Matrix of zeros and ones whose inverse has large integer entries

A = gallery('dramadah',n,k) returns an n-by-n matrix of 0's and 1's for which mu(A) = norm(inv(A),'fro') is relatively large, although not necessarily maximal. An anti-Hadamard matrix A is a matrix with elements 0 or 1 for which mu(A) is maximal.

n and k must both be scalars. Argument k determines the character of the output matrix:


k = 1
Default. A is Toeplitz, with abs(det(A)) = 1, and mu(A) > c(1.75)^n, where c is a constant. The inverse of A has integer entries.
k = 2
A is upper triangular and Toeplitz. The inverse of A has integer entries.
k = 3
A has maximal determinant among lower Hessenberg (0,1) matrices. det(A) = the nth Fibonacci number. A is Toeplitz. The eigenvalues have an interesting distribution in the complex plane.

fiedler--Symmetric matrix

A = gallery('fiedler',c), where c is a length n vector, returns the n-by-n symmetric matrix with elements abs(n(i)-n(j)). For scalar c, A = gallery('fiedler',1:c).

Matrix A has a dominant positive eigenvalue and all the other eigenvalues are negative.

Explicit formulas for inv(A) and det(A) are given in [Todd, J., Basic Numerical Mathematics, Vol. 2: Numerical Algebra, Birkhauser, Basel, and Academic Press, New York, 1977, p. 159] and attributed to Fiedler. These indicate that inv(A) is tridiagonal except for nonzero (1,n) and (n,1) elements.

forsythe--Perturbed Jordan block

A = gallery('forsythe',n,alpha,lambda) returns the n-by-n matrix equal to the Jordan block with eigenvalue lambda, excepting that A(n,1) = alpha. The default values of scalars alpha and lambda are sqrt(eps) and 0, respectively.

The characteristic polynomial of A is given by:

  • det(A-t*I) = (lambda-t)^N - alpha*(-1)^n.

frank--Matrix with ill-conditioned eigenvalues

F = gallery('frank',n,k) returns the Frank matrix of order n. It is upper Hessenberg with determinant 1. If k = 1, the elements are reflected about the anti-diagonal (1,n)--(n,1). The eigenvalues of F may be obtained in terms of the zeros of the Hermite polynomials. They are positive and occur in reciprocal pairs; thus if n is odd, 1 is an eigenvalue. F has floor(n/2) ill-conditioned eigenvalues--the smaller ones.

gearmat--Gear matrix

A = gallery('gearmat',n,i,j) returns the n-by-n matrix with ones on the sub- and super-diagonals, sign(i) in the (1,abs(i)) position, sign(j) in the (n,n+1-abs(j)) position, and zeros everywhere else. Arguments i and j default to n and -n, respectively.

Matrix A is singular, can have double and triple eigenvalues, and can be defective.

All eigenvalues are of the form 2*cos(a) and the eigenvectors are of the form [sin(w+a), sin(w+2*a), ..., sin(w+n*a)], where a and w are given in Gear, C. W., "A Simple Set of Test Matrices for Eigenvalue Programs", Math. Comp., Vol. 23 (1969), pp. 119-125.

grcar--Toeplitz matrix with sensitive eigenvalues

A = gallery('grcar',n,k) returns an n-by-n Toeplitz matrix with -1s on the subdiagonal, 1s on the diagonal, and k superdiagonals of 1s. The default is k = 3. The eigenvalues are sensitive.

hanowa--Matrix whose eigenvalues lie on a vertical line in the complex plane

A = gallery('hanowa',n,d) returns an n-by-n block 2-by-2 matrix of the form:

  • [d*eye(m)   -diag(1:m)
    diag(1:m) d*eye(m)]

Argument n is an even integer n=2*m. Matrix A has complex eigenvalues of the form d ± k*i, for 1 <= k <= m. The default value of d is -1.

house--Householder matrix

[v,beta,s] = gallery('house',x,k) takes x, an n-element column vector, and returns V and beta such that H*x = s*e1. In this expression, e1 is the first column of eye(n), abs(s) = norm(x), and H = eye(n) - beta*V*V' is a Householder matrix.

k determines the sign of s:


k = 0
sign(s) = -sign(x(1)) (default)
k = 1
sign(s) = sign(x(1))
k = 2
sign(s) = 1 (x must be real)

If x is complex, then sign(x) = x./abs(x) when x is nonzero.

If x = 0, or if x = alpha*e1 (alpha >= 0) and either k = 1 or k = 2, then V = 0, beta = 1, and s = x(1). In this case, H is the identity matrix, which is not strictly a Householder matrix.

invhess--Inverse of an upper Hessenberg matrix

A = gallery('invhess',x,y), where x is a length n vector and y is a length n-1 vector, returns the matrix whose lower triangle agrees with that of ones(n,1)*x' and whose strict upper triangle agrees with that of [1 y]*ones(1,n).

The matrix is nonsingular if x(1) ~= 0 and x(i+1) ~= y(i) for all i, and its inverse is an upper Hessenberg matrix. Argument y defaults to -x(1:n-1).

If x is a scalar, invhess(x) is the same as invhess(1:x).

invol--Involutory matrix

A = gallery('invol',n) returns an n-by-n involutory (A*A = eye(n)) and ill-conditioned matrix. It is a diagonally scaled version of hilb(n).

B = (eye(n)-A)/2 and B = (eye(n)+A)/2 are idempotent (B*B = B).

ipjfact--Hankel matrix with factorial elements

[A,d] = gallery('ipjfact',n,k) returns A, an n-by-n Hankel matrix, and d, the determinant of A, which is known explicitly. If k = 0 (the default), then the elements of A are A(i,j) = (i+j)! If k = 1, then the elements of A are A(i,j) = 1/(i+j).

Note that the inverse of A is also known explicitly.

jordbloc--Jordan block

A = gallery('jordbloc',n,lambda) returns the n-by-n Jordan block with eigenvalue lambda. The default value for lambda is 1.

kahan--Upper trapezoidal matrix

A = gallery('kahan',n,theta,pert) returns an upper trapezoidal matrix that has interesting properties regarding estimation of condition and rank.

If n is a two-element vector, then A is n(1)-by-n(2); otherwise, A is n-by-n. The useful range of theta is 0 <>, with a default value of 1.2.

To ensure that the QR factorization with column pivoting does not interchange columns in the presence of rounding errors, the diagonal is perturbed by pert*eps*diag([n:-1:1]). The default pert is 25, which ensures no interchanges for gallery('kahan',n) up to at least n = 90 in IEEE arithmetic.

kms--Kac-Murdock-Szego Toeplitz matrix

A = gallery('kms',n,rho) returns the n-by-n Kac-Murdock-Szego Toeplitz matrix such that A(i,j) = rho^(abs(i-j)), for real rho.

For complex rho, the same formula holds except that elements below the diagonal are conjugated. rho defaults to 0.5.

The KMS matrix A has these properties:

  • An LDL' factorization with L = inv(gallery('triw',n,-rho,1))', and D(i,i) = (1-abs(rho)^2)*eye(n), except D(1,1) = 1.
  • Positive definite if and only if 0 <>.
  • The inverse inv(A) is tridiagonal.

krylov--Krylov matrix

B = gallery('krylov',A,x,j) returns the Krylov matrix

  • [x, Ax, A^2x, ..., A^(j-1)x]

where A is an n-by-n matrix and x is a length n vector. The defaults are x = ones(n,1), and j = n.

B = gallery('krylov',n) is the same as gallery('krylov',(randn(n)).

lauchli--Rectangular matrix

A = gallery('lauchli',n,mu) returns the (n+1)-by-n matrix

  • [ones(1,n); mu*eye(n)]

The Lauchli matrix is a well-known example in least squares and other problems that indicates the dangers of forming A'*A. Argument mu defaults to sqrt(eps).

lehmer--Symmetric positive definite matrix

A = gallery('lehmer',n) returns the symmetric positive definite n-by-n matrix such that A(i,j) = i/j for j >= i.

The Lehmer matrix A has these properties:

  • A is totally nonnegative.
  • The inverse inv(A) is tridiagonal and explicitly known.
  • The order n <= cond(A) <= 4*n*n.

leslie--

L = gallery('leslie',a,b) is the n-by-n matrix from the Leslie population model with average birth numbers a(1:n) and survival rates b(1:n-1). It is zero, apart from the first row (which contains the a(i)) and the first subdiagonal (which contains the b(i)). For a valid model, the a(i) are nonnegative and the b(i) are positive and bounded by 1, i.e., 0 <>.

L = gallery('leslie',n) generates the Leslie matrix with a = ones(n,1), b = ones(n-1,1).

lesp--Tridiagonal matrix with real, sensitive eigenvalues

A = gallery('lesp',n) returns an n-by-n matrix whose eigenvalues are real and smoothly distributed in the interval approximately [-2*N-3.5, -4.5].

The sensitivities of the eigenvalues increase exponentially as the eigenvalues grow more negative. The matrix is similar to the symmetric tridiagonal matrix with the same diagonal entries and with off-diagonal entries 1, via a similarity transformation with D = diag(1!,2!,...,n!).

lotkin--Lotkin matrix

A = gallery('lotkin',n) returns the Hilbert matrix with its first row altered to all ones. The Lotkin matrix A is nonsymmetric, ill-conditioned, and has many negative eigenvalues of small magnitude. Its inverse has integer entries and is known explicitly.

minij--Symmetric positive definite matrix

A = gallery('minij',n) returns the n-by-n symmetric positive definite matrix with A(i,j) = min(i,j).

The minij matrix has these properties:

  • The inverse inv(A) is tridiagonal and equal to -1 times the second difference matrix, except its (n,n) element is 1.
  • Givens' matrix, 2*A-ones(size(A)), has tridiagonal inverse and eigenvalues 0.5*sec((2*r-1)*pi/(4*n))^2, where r=1:n.
  • (n+1)*ones(size(A))-A has elements that are max(i,j) and a tridiagonal inverse.

moler--Symmetric positive definite matrix

A = gallery('moler',n,alpha) returns the symmetric positive definite n-by-n matrix U'*U, where U = gallery('triw',n,alpha).

For the default alpha = -1, A(i,j) = min(i,j)-2, and A(i,i) = i. One of the eigenvalues of A is small.

neumann--Singular matrix from the discrete Neumann problem (sparse)

C = gallery('neumann',n) returns the sparse n-by-n singular, row diagonally dominant matrix resulting from discretizing the Neumann problem with the usual five-point operator on a regular mesh. Argument n is a perfect square integer or a two-element vector. C is sparse and has a one-dimensional null space with null vector ones(n,1).

orthog--Orthogonal and nearly orthogonal matrices

Q = gallery('orthog',n,k) returns the kth type of matrix of order n, where k > 0 selects exactly orthogonal matrices, and k <> selects diagonal scalings of orthogonal matrices. Available types are:


k = 1
Q(i,j) = sqrt(2/(n+1)) * sin(i*j*pi/(n+1))
Symmetric eigenvector matrix for second difference matrix. This is the default.
k = 2
Q(i,j) = 2/(sqrt(2*n+1)) * sin(2*i*j*pi/(2*n+1))
Symmetric.
k = 3
Q(r,s) = exp(2*pi*i*(r-1)*(s-1)/n) / sqrt(n)
Unitary, the Fourier matrix. Q^4 is the identity. This is essentially the same matrix as fft(eye(n))/sqrt(n)!
k = 4
Helmert matrix: a permutation of a lower Hessenberg matrix, whose first row is ones(1:n)/sqrt(n).
k = 5
Q(i,j) = sin(2*pi*(i-1)*(j-1)/n) + cos(2*pi*(i-1)*(j-1)/n)
Symmetric matrix arising in the Hartley transform.
K = 6
Q(i,j) = sqrt(2/n)*cos((i-1/2)*(j-1/2)*pi/n)
Symmetric matrix arising as a discrete cosine transform.
k = -1
Q(i,j) = cos((i-1)*(j-1)*pi/(n-1))
Chebyshev Vandermonde-like matrix, based on extrema of T(n-1).
k = -2
Q(i,j) = cos((i-1)*(j-1/2)*pi/n))
Chebyshev Vandermonde-like matrix, based on zeros of T(n).

parter--Toeplitz matrix with singular values near pi

C = gallery('parter',n) returns the matrix C such that C(i,j) = 1/(i-j+0.5).

C is a Cauchy matrix and a Toeplitz matrix. Most of the singular values of C are very close to pi.

pei--Pei matrix

A = gallery('pei',n,alpha), where alpha is a scalar, returns the symmetric matrix alpha*eye(n) + ones(n). The default for alpha is 1. The matrix is singular for alpha equal to either 0 or -n.

poisson--Block tridiagonal matrix from Poisson's equation (sparse)

A = gallery('poisson',n) returns the block tridiagonal (sparse) matrix of order n^2 resulting from discretizing Poisson's equation with the 5-point operator on an n-by-n mesh.

prolate--Symmetric, ill-conditioned Toeplitz matrix

A = gallery('prolate',n,w) returns the n-by-n prolate matrix with parameter w. It is a symmetric Toeplitz matrix.

If 0 <> then A is positive definite

  • The eigenvalues of A are distinct, lie in (0,1), and tend to cluster around 0 and 1.
  • The default value of w is 0.25.

randcolu -- Random matrix with normalized cols and specified singular values

A = gallery('randcolu',n) is a random n-by-n matrix with columns of unit 2-norm, with random singular values whose squares are from a uniform distribution.

A'*A is a correlation matrix of the form produced by gallery('randcorr',n).

gallery('randcolu',x) where x is an n-vector (n > 1), produces a random n-by-n matrix having singular values given by the vector x. The vector x must have nonnegative elements whose sum of squares is n.

gallery('randcolu',x,m) where m >= n, produces an m-by-n matrix.

gallery('randcolu',x,m,k) provides a further option:


k = 0
diag(x) is initially subjected to a random two-sided orthogonal transformation, and then a sequence of Givens rotations is applied (default).
k = 1
The initial transformation is omitted. This is much faster, but the resulting matrix may have zero entries.

randcorr -- Random correlation matrix with specified eigenvalues

gallery('randcorr',n) is a random n-by-n correlation matrix with random eigenvalues from a uniform distribution. A correlation matrix is a symmetric positive semidefinite matrix with 1s on the diagonal.

gallery('randcorr',x) produces a random correlation matrix having eigenvalues given by the vector x, where length(x) > 1. The vector x must have nonnegative elements summing to length(x).

gallery('randcorr',x,k) provides a further option:


k = 0
The diagonal matrix of eigenvalues is initially subjected to a random orthogonal similarity transformation, and then a sequence of Givens rotations is applied (default).
k = 1
The initial transformation is omitted. This is much faster, but the resulting matrix may have some zero entries.

For more information, see:

[1] Bendel, R. B. and M. R. Mickey, "Population Correlation Matrices for Sampling Experiments," Commun. Statist. Simulation Comput., B7, 1978, pp. 163-182.

[2] Davies, P. I. and N. J. Higham, "Numerically Stable Generation of Correlation Matrices and Their Factors," BIT, Vol. 40, 2000, pp. 640-651.

randhess--Random, orthogonal upper Hessenberg matrix

H = gallery('randhess',n) returns an n-by-n real, random, orthogonal upper Hessenberg matrix.

H = gallery('randhess',x) if x is an arbitrary, real, length n vector with n > 1, constructs H nonrandomly using the elements of x as parameters.

Matrix H is constructed via a product of n-1 Givens rotations.

rando--Random matrix composed of elements -1, 0 or 1

A = gallery('rando',n,k) returns a random n-by-n matrix with elements from one of the following discrete distributions:


k = 1
A(i,j) = 0 or 1 with equal probability (default).
k = 2
A(i,j) = -1 or 1 with equal probability.
k = 3
A(i,j) = -1, 0 or 1 with equal probability.

Argument n may be a two-element vector, in which case the matrix is n(1)-by-n(2).

randsvd--Random matrix with preassigned singular values

A = gallery('randsvd',n,kappa,mode,kl,ku) returns a banded (multidiagonal) random matrix of order n with cond(A) = kappa and singular values from the distribution mode. If n is a two-element vector, A is n(1)-by-n(2).

Arguments kl and ku specify the number of lower and upper off-diagonals, respectively, in A. If they are omitted, a full matrix is produced. If only kl is present, ku defaults to kl.

Distribution mode can be:


1
One large singular value.
2
One small singular value.
3
Geometrically distributed singular values (default).
4
Arithmetically distributed singular values.
5
Random singular values with uniformly distributed logarithm.
<>
If mode is -1, -2, -3, -4, or -5, then randsvd treats mode as abs(mode), except that in the original matrix of singular values the order of the diagonal entries is reversed: small to large instead of large to small.

Condition number kappa defaults to sqrt(1/eps). In the special case where kappa <>, A is a random, full, symmetric, positive definite matrix with cond(A) = -kappa and eigenvalues distributed according to mode. Arguments kl and ku, if present, are ignored.

A = gallery('randsvd',n,kappa,mode,kl,ku,method) specifies how the computations are carried out. method = 0 is the default, while method = 1 uses an alternative method that is much faster for large dimensions, even though it uses more flops.

redheff--Redheffer's matrix of 1s and 0s

A = gallery('redheff',n) returns an n-by-n matrix of 0's and 1's defined by A(i,j) = 1, if j = 1 or if i divides j, and A(i,j) = 0 otherwise.

The Redheffer matrix has these properties:

  • (n-floor(log2(n)))-1 eigenvalues equal to 1
  • A real eigenvalue (the spectral radius) approximately sqrt(n)
  • A negative eigenvalue approximately -sqrt(n)
  • The remaining eigenvalues are provably "small."
  • The Riemann hypothesis is true if and only if

    for every epsilon > 0.

Barrett and Jarvis conjecture that "the small eigenvalues all lie inside the unit circle abs(Z) = 1," and a proof of this conjecture, together with a proof that some eigenvalue tends to zero as n tends to infinity, would yield a new proof of the prime number theorem.

riemann--Matrix associated with the Riemann hypothesis

A = gallery('riemann',n) returns an n-by-n matrix for which the Riemann hypothesis is true if and only if

for every .

The Riemann matrix is defined by:

  • A = B(2:n+1,2:n+1)

where B(i,j) = i-1 if i divides j, and B(i,j) = -1 otherwise.

The Riemann matrix has these properties:

  • Each eigenvalue e(i) satisfies abs(e(i)) <= m-1/m, where m = n+1.
  • i <= e(i) <= i+1 with at most m-sqrt(m) exceptions.
  • All integers in the interval (m/3, m/2] are eigenvalues.

ris--Symmetric Hankel matrix

A = gallery('ris',n) returns a symmetric n-by-n Hankel matrix with elements

  • A(i,j) = 0.5/(n-i-j+1.5)

The eigenvalues of A cluster around and . This matrix was invented by F.N. Ris.

rosser--Classic symmetric eigenvalue test matrix

A = rosser returns the Rosser matrix. This matrix was a challenge for many matrix eigenvalue algorithms. But the QR algorithm, as perfected by Wilkinson and implemented in MATLAB, has no trouble with it. The matrix is 8-by-8 with integer elements. It has:

  • A double eigenvalue
  • Three nearly equal eigenvalues
  • Dominant eigenvalues of opposite sign
  • A zero eigenvalue
  • A small, nonzero eigenvalue

smoke--Complex matrix with a 'smoke ring' pseudospectrum

A = gallery('smoke',n) returns an n-by-n matrix with 1's on the superdiagonal, 1 in the (n,1) position, and powers of roots of unity along the diagonal.

A = gallery('smoke',n,1) returns the same except that element A(n,1) is zero.

The eigenvalues of gallery('smoke',n,1) are the nth roots of unity; those of gallery('smoke',n) are the nth roots of unity times 2^(1/n).

toeppd--Symmetric positive definite Toeplitz matrix

A = gallery('toeppd',n,m,w,theta) returns an n-by-n symmetric, positive semi-definite (SPD) Toeplitz matrix composed of the sum of m rank 2 (or, for certain theta, rank 1) SPD Toeplitz matrices. Specifically,

  • T = w(1)*T(theta(1)) + ... + w(m)*T(theta(m))

where T(theta(k)) has (i,j) element cos(2*pi*theta(k)*(i-j)).

By default: m = n, w = rand(m,1), and theta = rand(m,1).

toeppen--Pentadiagonal Toeplitz matrix (sparse)

P = gallery('toeppen',n,a,b,c,d,e) returns the n-by-n sparse, pentadiagonal Toeplitz matrix with the diagonals: P(3,1) = a, P(2,1) = b, P(1,1) = c, P(1,2) = d, and P(1,3) = e, where a, b, c, d, and e are scalars.

By default, (a,b,c,d,e) = (1,-10,0,10,1), yielding a matrix of Rutishauser. This matrix has eigenvalues lying approximately on the line segment 2*cos(2*t) + 20*i*sin(t).

tridiag--Tridiagonal matrix (sparse)

A = gallery('tridiag',c,d,e) returns the tridiagonal matrix with subdiagonal c, diagonal d, and superdiagonal e. Vectors c and e must have length(d)-1.

A = gallery('tridiag',n,c,d,e), where c, d, and e are all scalars, yields the Toeplitz tridiagonal matrix of order n with subdiagonal elements c, diagonal elements d, and superdiagonal elements e. This matrix has eigenvalues

  • d + 2*sqrt(c*e)*cos(k*pi/(n+1))

where k = 1:n. (see [1].)

A = gallery('tridiag',n) is the same as A = gallery('tridiag',n,-1,2,-1), which is a symmetric positive definite M-matrix (the negative of the second difference matrix).

triw--Upper triangular matrix discussed by Wilkinson and others

A = gallery('triw',n,alpha,k) returns the upper triangular matrix with ones on the diagonal and alphas on the first k >= 0 superdiagonals.

Order n may be a 2-element vector, in which case the matrix is n(1)-by-n(2) and upper trapezoidal.

Ostrowski ["On the Spectrum of a One-parametric Family of Matrices, J. Reine Angew. Math., 1954] shows that

  • cond(gallery('triw',n,2)) = cot(pi/(4*n))^2,

and, for large abs(alpha), cond(gallery('triw',n,alpha)) is approximately abs(alpha)^n*sin(pi/(4*n-2)).

Adding -2^(2-n) to the (n,1) element makes triw(n) singular, as does adding
-2^(1-n) to all the elements in the first column.

vander--Vandermonde matrix

A = gallery('vander',c) returns the Vandermonde matrix whose second to last column is c. The jth column of a Vandermonde matrix is given by A(:,j) = C^(n-j).

wathen--Finite element matrix (sparse, random entries)

A = gallery('wathen',nx,ny) returns a sparse, random, n-by-n finite element matrix where n = 3*nx*ny + 2*nx + 2*ny + 1.

Matrix A is precisely the "consistent mass matrix" for a regular nx-by-ny grid of 8-node (serendipity) elements in two dimensions. A is symmetric, positive definite for any (positive) values of the "density," rho(nx,ny), which is chosen randomly in this routine.

A = gallery('wathen',nx,ny,1) returns a diagonally scaled matrix such that

  • 0.25 <= eig(inv(D)*A) <= 4.5

where D = diag(diag(A)) for any positive integers nx and ny and any densities rho(nx,ny).

wilk--Various matrices devised or discussed by Wilkinson

[A,b] = gallery('wilk',n) returns a different matrix or linear system depending on the value of n.


n = 3
Upper triangular system Ux=b illustrating inaccurate solution.
n = 4
Lower triangular system Lx=b, ill-conditioned.
n = 5
hilb(6)(1:5,2:6)*1.8144. A symmetric positive definite matrix.
n = 21
W21+, a tridiagonal matrix. Eigenvalue problem.