gradp.src

没有说明

/*
** gradp.src
** (C) Copyright 1988-1998 by Aptech Systems, Inc.
** All Rights Reserved.
**
** This Software Product is PROPRIETARY SOURCE CODE OF APTECH
** SYSTEMS, INC.    This File Header must accompany all files using
** any portion, in whole or in part, of this Source Code.   In
** addition, the right to create such files is strictly limited by
** Section 2.A. of the GAUSS Applications License Agreement
** accompanying this Software Product.
**
** If you wish to distribute any portion of the proprietary Source
** Code, in whole or in part, you must first obtain written
** permission from Aptech Systems.
**
**> gradp
**
**  Purpose:    Computes the gradient vector or matrix (Jacobian) of a
**              vector-valued function that has been defined in a procedure.
**              Single-sided (forward difference) gradients are computed.
**
**  Format:     g = gradp(&f,x0);
**
**  Input:      f    scalar, procedure pointer to a vector-valued function:
**
**                                          f:Kx1 -> Nx1
**
**                   It is acceptable for f(x) to have been defined in terms of
**                   global arguments in addition to x, and thus f can return
**                   an Nx1 vector:
**
**                        proc f(x);
**                           retp( exp(x*b) );
**                        endp;
**
**              x0    Kx1 vector of points at which to compute gradient.
**
**  Output:     g     NxK matrix containing the gradients of f with respect
**                    to the variable x at x0.
**
**  Remarks:    gradp will return a row for every row that is returned by f.
**              For instance, if f returns a 1x1 result, then gradp will
**              return a 1xK row vector. This allows the same function to be
**              used where N is the number of rows in the result returned by f.
**              Thus, for instance, gradp can be used to compute the
**              Jacobian matrix of a set of equations.
**
**  Example:    proc myfunc(x);
**                 retp( x .* 2 .* exp( x .* x ./ 3 ));
**              endp;
**
**              x0 = { 2.5, 3.0, 3.5 };
**              y = gradp(&myfunc,x0);
**
**                           82.98901842    0.00000000    0.00000000
**                  y =       0.00000000  281.19752975    0.00000000
**                            0.00000000    0.00000000 1087.95414117
**
**              It is a 3x3 matrix because we are passing it 3 arguments and
**              myfunc returns 3 results when we do that.  The off-diagonals
**              are zeros because the cross-derivatives of 3 arguments are 0.
**
**  Globals:    None
**
**  See Also:   hessp
*/

proc 1 = gradp(f,x0);
    local f:proc;
    local n,k,grdd,dh,ax0,xdh,arg,dax0,i,f0;

    /* check for complex input */
    if iscplx(x0);
        if hasimag(x0);
            errorlog "ERROR: Not implemented for complex matrices.";
            end;
        else;
            x0 = real(x0);
        endif;
    endif;

    f0 = f(x0);
    n = rows(f0);
    k = rows(x0);
    grdd = zeros(n,k);

/* Computation of stepsize (dh) for gradient */

    ax0 = abs(x0);
    if x0 /= 0;
        dax0 = x0./ax0;
    else;
        dax0 = 1;
    endif;
    dh = (1e-8)*maxc((ax0~(1e-2)*ones(rows(x0),1))').*dax0;
    xdh = x0+dh;
    dh = xdh-x0;    /* This increases precision slightly */
    arg = diagrv(reshape(x0,k,k)',xdh);

    i = 1;
    do until i > k;
        grdd[.,i] = f(arg[.,i]);
        i = i+1;
    endo;

    grdd = (grdd-f0)./(dh');

    retp(grdd);
endp;