functions in math.i -

             _dgecox  
 
     LAPACK dgecon routine, except norm argument not a string.  

             _dgelss  
 
     LAPACK dgelss routine.  

             _dgelx  
 
     LAPACK dgels routine, except trans argument not a string.  

             _dgesv  
 
     LAPACK dgesv routine.  

             _dgesvx  
 
     LAPACK dgesvd routine, except jobu and jobvt are not strings.  

             _dgetrf  
 
     LAPACK dgetrf routine.  Performs LU factorization.  

             _dgtsv  
 
     LAPACK dgtsv routine.  

             fft(x, direction)  
             fft(x, ljdir, rjdir)  
          or fft(x, ljdir, rjdir, setup=workspace)  
 
     returns the complex Fast Fourier Transform of array X.  
     The DIRECTION determines which direction the transform is in --  
     e.g.- from time to frequency or vice-versa -- as follows:  
     DIRECTION    meaning  
     ---------    -------  
         1        "forward" transform (coefficients of exp(+i * 2*pi*kl/N))  
                  on every dimension of X  
        -1        "backward" transform (coefficients of exp(-i * 2*pi*kl/N))  
                  on every dimension of X  
     [1,-1,1]     forward transform on first and third dimensions of X,  
                  backward transform on second dimension of X (any other  
                  dimensions remain untransformed)  
     [-1,0,0,1]   backward transform on first dimension of X, forward  
                  transform on fourth dimension of X  
        etc.  
     The third positional argument, if present, allows the direction  
     of dimensions of X to be specified relative to the final dimension  
     of X, instead of relative to the first dimension of X.  In this  
     case, both LJDIR and RJDIR must be vectors of integers -- the  
     scalar form is illegal:  
        LJDIR    RJDIR      meaning  
        -----    -----      -------  
        []        [1]       forward transform last dimension of X  
        [1]        []       forward transform first dimension of X  
        []        [-1,-1]   backward transform last two dimensions of X,  
                            leaving any other dimensions untransformed  
     [-1,0,0,1]    []       backward transform on first dimension of X,  
                            forward transform on fourth dimension of X  
        []      [-1,0,0,1]  backward transform on 4th to last dimension of X,  
                            forward transform on last dimension of X  
        etc.  
     Note that the final element of RJDIR corresponds to the last dimension  
     of X, while the initial element of LJDIR corresponds to the first  
     dimension of X.  
     The explicit meaning of "forward" transform -- the coefficients of  
     exp(+i * 2*pi*kl/N) -- is:  
     result for j=1,...,n  
                result(j)=the sum from k=1,...,n of  
                      x(k)*exp(-i*(j-1)*(k-1)*2*pi/n)  
                            where i=sqrt(-1)  
     Note that the result is unnormalized.  Applying the "backward"  
     transform to the result of a "forward" transform returns N times  
     the original vector of length N.  Equivalently, applying either  
     the "forward" or "backward" transform four times in succession  
     yields N^2 times the original vector of length N.  
     Performing the transform requires some WORKSPACE, which can be  
     set up beforehand by calling fft_setup, if fft is to be called  
     more than once with arrays X of the same shape.  If no setup  
     keyword argument is supplied, the workspace allocation and setup  
     must be repeated for each call.  

             fft_braw, n, c, wsave  
 
     Swarztrauber's cfftb.  You can use this to avoid the additional  
     2*N storage incurred by fft_setup.  

             fft_fraw, n, c, wsave  
 
     Swarztrauber's cfftf.  You can use this to avoid the additional  
     2*N storage incurred by fft_setup.  

             fft_init, n, wsave  
 
     Swarztrauber's cffti.  This actually requires wsave=array(0.0, 4*n+15),  
     instead of the 6*n+15 doubles of storage used by fft_raw to handle the  
     possibility of multidimensional arrays.  If the storage matters, you  
     can call cfftf and/or cfftb as the Yorick functions fft_fraw and/or  
     fft_braw.  

             fft_inplace, x, direction  
          or fft_inplace, x, ljdir, rjdir  
          or fft_inplace, x, ljdir, rjdir, setup=workspace  
 
     is the same as the fft function, except that the transform is  
     performed "in_place" on the array X, which must be of type complex.  

             workspace= fft_setup(dimsof(x))  
          or workspace= fft_setup(dimsof(x), direction)  
          or workspace= fft_setup(dimsof(x), ljdir, rjdir)  
 
     allocates and sets up the workspace for a subsequent call to  
            fft(X, DIRECTION, setup=WORKSPACE)  
     or  
            fft(X, LJDIR, RJDIR, setup=WORKSPACE)  
     The DIRECTION or LJDIR, RJDIR arguments compute WORKSPACE only for  
     the dimensions which will actually be transformed.  If only the  
     dimsof(x) argument is supplied, then WORKSPACE will be enough to  
     transform any or all dimensions of X.  With DIRECTION or LJDIR, RJDIR  
     supplied, WORKSPACE will only be enough to compute the dimensions  
     which are actually to be transformed.  The WORKSPACE does not  
     depend on the sign of any element in the DIRECTION (or LJDIR, RJDIR),  
     so you can use the same WORKSPACE for both "forward" and "backward"  
     transforms.  
     Furthermore, as long as the length of any dimensions of the array  
     X to be transformed are present in WORKSPACE, it may be used in  
     a call to fft with the array.  Thus, if X were a 25-by-64 array,  
     and Y were a 64-vector, the following sequence is legal:  
          ws= fft_setup(dimsof(x));  
          xf= fft(x, 1, setup=ws);  
          yf= fft(y, -1, setup=ws);  
     The WORKSPACE required for a dimension of length N is 6*N+15 doubles.  

             LUrcond(a)  
          or LUrcond(a, one_norm=1)  
 
     returns the reciprocal condition number of the N-by-N matrix A.  
     If the ONE_NORM argument is non-nil and non-zero, the 1-norm  
     condition number is returned, otherwise the infinity-norm condition  
     number is returned.  
     The condition number is the ratio of the largest to the smallest  
     singular value, max(singular_values)*max(1/singular_values) (or  
     sum(abs(singular_values)*sum(abs(1/singular_values)) if ONE_NORM  
     is selected?).  If the reciprocal condition number is near zero  
     then A is numerically singular; specifically, if  
          1.0+LUrcond(a) == 1.0  
     then A is numerically singular.  

             LUsolve(a, b)  
          or LUsolve(a, b, which=which)  
          or a_inverse= LUsolve(a)  
 
     returns the solution x of the matrix equation:  
        A(,+)*x(+) = B  
     If A is an n-by-n matrix then B must have length n, and the returned  
     x will also have length n.  
     B may have additional dimensions, in which case the returned x  
     will have the same additional dimensions.  The WHICH dimension of B,  
     and of the returned x is the one of length n which participates  
     in the matrix solve.  By default, WHICH=1, so that the equations  
     being solved are:  
        A(,+)*x(+,..) = B  
     Non-positive WHICH counts from the final dimension (as for the  
     sort and transpose functions), so that WHICH=0 solves:  
        x(..,+)*A(,+) = B  
     Other examples:  
        A_ij X_jklm = B_iklm   (WHICH=1)  
        A_ij X_kjlm = B_kilm   (WHICH=2)  
        A_ij X_klmj = B_klmi   (WHICH=4 or WHICH=0)  
     If the B argument is omitted, the inverse of A is returned:  
     A(,+)*x(+,) and A(,+)*x(,+) will be unit matrices.  
     LUsolve works by LU decomposition using Gaussian elimination with  
     pivoting.  It is the fastest way to solve square matrices.  QRsolve  
     handles non-square matrices, as does SVsolve.  SVsolve is slowest,  
     but can deal with highly singular matrices sensibly.  

             QRsolve(a, b)  
          or QRsolve(a, b, which=which)  
 
     returns the solution x (in a least squares or least norm sense  
     described below) of the matrix equation:  
        A(,+)*x(+) = B  
     If A is an m-by-n matrix (i.e.- m equations in n unknowns), then B  
     must have length m, and the returned x will have length n.  
     If nm, the system is underdetermined -- many solutions are possible  
             -- the returned x has minimum L2 norm among all solutions  
     B may have additional dimensions, in which case the returned x  
     will have the same additional dimensions also have those dimensions.  
     The WHICH dimension of B and the returned x is the one of length m  
     or n which participates in the matrix solve.  By default, WHICH=1,  
     so that the equations being solved are:  
        A(,+)*x(+,..) = B  
     Non-positive WHICH counts from the final dimension (as for the  
     sort and transpose functions), so that WHICH=0 solves:  
        A(,+)*x(..,+) = B  
     QRsolve works by QR factorization if nm.  
     QRsolve is slower than LUsolve.  Its main attraction is that it can  
     handle overdetermined or underdetermined systems of equations  
     (nonsquare matrices).  QRsolve may fail for singular systems; try  
     SVsolve in this case.  

             roll(x, ljoff, rjoff)  
          or roll, x, ljoff, rjoff  
          or roll(x)  
          or roll, x  
 
     "rolls" selected dimensions of the array X.  The roll offsets  
     LJOFF and RJOFF (both optional) work in the same fashion as the  
     LJDIR and RJDIR arguments to the fft function:  
        A scalar LJDIR (and nil RJDIR) rolls all dimensions of X by  
        the specified offset.  
        Otherwise, the elements of the LJDIR vector [ljoff1, ljoff2, ...]  
        are used as the roll offsets for the first, second, etc.  
        dimensions of X.  
        Similarly, the elements of the RJDIR vector [..., rjoff1, rjoff0]  
        are matched to the final dimensions of X, so the next to last  
        dimension is rolled by rjoff1 and the last dimension by rjoff0.  
        As a special case (mostly for use with the fft function), if  
        both LJDIR and RJDIR are nil, every dimension is rolled by  
        half of its length.  Thus,  
           roll(x)  
        it equivalent to  
           roll(x, dimsof(x)(2:0)/2)  
     The result of the roll function is complex if X is complex, otherwise  
     double (i.e.- any other array type is promoted to type double).  If  
     roll is invoked as a subroutine, the operation is performed in place.  

             s= SVdec(a, u, vt)  
          or s= SVdec(a, u, vt, full=1)  
 
     performs the singular value decomposition of the m-by-n matrix A:  
        A = (U(,+) * SIGMA(+,))(,+) * VT(+,)  
     where U is an m-by-m orthogonal matrix, VT is an n-by-n orthogonal  
     matrix, and SIGMA is an m-by-n matrix which is zero except for its  
     min(m,n) diagonal elements.  These diagonal elements are the return  
     value of the function, S.  The returned S is always arranged in  
     order of descending absolute value.  U(,1:min(m,n)) are the left  
     singular vectors corresponding to the min(m,n) elements of S;  
     VT(1:min(m,n),) are the right singular vectors.  (The original A  
     matrix maps a right singular vector onto the corresponding left  
     singular vector, stretched by a factor of the singular value.)  
     Note that U and VT are strictly outputs; if you don't need them,  
     they need not be present in the calling sequence.  
     By default, U will be an m-by-min(m,n) matrix, and V will be  
     a min(m,n)-by-n matrix (i.e.- only the singular vextors are returned,  
     not the full orthogonal matrices).  Set the FULL keyword to a  
     non-zero value to get the full m-by-m and n-by-n matrices.  
     On rare occasions, the routine may fail; if it does, the  
     first SVinfo values of the returned S are incorrect.  Hence,  
     the external variable SVinfo will be 0 after a successful call  
     to SVdec.  If SVinfo>0, then external SVe contains the superdiagonal  
     elements of the bidiagonal matrix whose diagonal is the returned  
     S, and that bidiagonal matrix is equal to (U(+,)*A(+,))(,+) * V(+,).  
     Numerical Recipes (Press, et. al. Cambridge University Press 1988)  
     has a good discussion of how to use the SVD -- see section 2.9.  

             SVsolve(a, b)  
          or SVsolve(a, b, rcond)  
          or SVsolve(a, b, rcond, which=which)  
 
     returns the solution x (in a least squares sense described below) of  
     the matrix equation:  
        A(,+)*x(+) = B  
     If A is an m-by-n matrix (i.e.- m equations in n unknowns), then B  
     must have length m, and the returned x will have length n.  
     If nm, the system is underdetermined -- many solutions are possible  
             -- the returned x has minimum L2 norm among all solutions  
     SVsolve works by singular value decomposition, therefore it is  
     immune to failure due to singularity of the A matrix.  The optional  
     RCOND argument defaults to 1.0e-9; singular values less than RCOND  
     times the largest singular value (absolute value) will be set to 0.0.  
     If RCOND<=0.0, machine precision is used.  The effective rank of the  
     matrix is returned as the external variable SVrank.  
     You can examine the details of the SVD by calling the SVdec routine,  
     which returns the singular vectors as well as the singular values.  
     Numerical Recipes (Press, et. al. Cambridge University Press 1988)  
     has a good discussion of how to use the SVD -- see section 2.9.  
     B may have additional dimensions, in which case the returned x  
     will have the same additional dimensions.  The WHICH argument  
     (default 1) controls which dimension of B takes part in the matrix  
     solve.  See QRsolve or LUsolve for a complete discussion.  

             TDsolve(c, d, e, b)  
          or TDsolve(c, d, e, b, which=which)  
 
     returns the solution to the tridiagonal system:  
        D(1)*x(1)       + E(1)*x(2)                       = B(1)  
        C(1:-1)*x(1:-2) + D(2:-1)*x(2:-1) + E(2:0)*x(3:0) = B(2:-1)  
                          C(0)*x(-1)      + D(0)*x(0)     = B(0)  
     (i.e.- C is the subdiagonal, D the diagonal, and E the superdiagonal;  
     C and E have one fewer element than D, which is the same length as  
     both B and x)  
     B may have additional dimensions, in which case the returned x  
     will have the same additional dimensions.  The WHICH dimension of B,  
     and of the returned x is the one of length n which participates  
     in the matrix solve.  By default, WHICH=1, so that the equations  
     being solved involve B(,..) and x(+,..).  
     Non-positive WHICH counts from the final dimension (as for the  
     sort and transpose functions), so that WHICH=0 involves B(..,)  
     and x(..,+).  
     The C, D, and E arguments may be either scalars or vectors; they  
     will be broadcast as appropriate.  

             unit(n)  
          or unit(n, m)  
 
     returns n-by-n (or n-by-m) unit matrix, i.e.- matrix with diagonal  
     elements all 1.0, off diagonal elements 0.0