Scilab Bag Of Tricks: The Scilab-2.5 IAQ (Infrequently Asked Questions) | ||
---|---|---|

Prev | Chapter 7. Scilab Core | Next |

In the following two sections we shall treat the
"anatomy" of native, i.e. low-level Scilab functions.
This will confront us with all the gory details of the stack, the
low-level API, and the calling conventions. Having the
"Guide for Developers",
`Internals.ps` (see also the section called *Local Documents* in Chapter 8) ready is a good idea. Where the
developer guide is at the end of its wits, a study of the source
code is appropriate, especially the file
`SCI``/routines/interf/stack1.f`

We start out discussing simple functions. Simple in the sense that they are self-contained and only take non-function parameters as their arguments. In the second part we shall consider functions that take other functions (either Scilab functions or externals) as arguments.

A typical native Scilab function proceeds as follows:

Check the number of input and output parameters.

Get the "pointers" to all actual input parameters; supply default values for optional parameters; issue warnings or errors as appropriate for too many or too few parameters.

Allocate space for all temporary variables, "workspaces", etc.

It might be necessary to translate the input variables which are in Scilab format into the appropriate format for the worker routine. This is necessary for example if the worker routine uses

`DOUBLE COMPLEX`(or`COMPLEX*16`) variables.Perform the calculations or transformations that

*really*make up the procedure.As in Step 4, it might be necessary to transform the results, now form the worker routine's format back into Scilab format.

If necessary, allocate space on the Scilab stack and copy results to this space.

Now that the outline is clear we are ready to dissect a simple
function: `ortho`. The function takes exactly
one argument ` a`, that is a real or
complex

type Complex is record Re, Im : Float'Base; end record; type FloatVector is array (Positive range <>) of Float; type ComplexVector is array (Positive range <>) of Complex; type FloatMatrix is array (1..Lda, Positive range <>) of Float; type ComplexMatrix is array (1..Lda, Positive range <>) of Complex; proceduredgeqrf(: in Natural;M: in Natural;N: in out FloatMatrix;A: in Natural;Lda: out FloatVector;Tau: out FloatVector;Work: in Integer;Lwork: out Integer); procedureInfodorgqr(: in Natural;M: in Natural;N: in Natural;K: in out FloatMatrix;A: in Natural;Lda: out FloatVector;Tau: out FloatVector;Work: in Integer;Lwork: out Integer); procedureInfozgeqrf(: in Natural;M: in Natural;N: in out ComplexMatrix;A: in Natural;Lda: out ComplexVector;Tau: out ComplexVector;Work: in Integer;Lwork: out Integer); procedureInfozungqr(: in Natural;M: in Natural;N: in Natural;K: in out ComplexMatrix;A: in Natural;Lda: out ComplexVector;Tau: out ComplexVector;Work: in Integer;Lwork: out Integer); procedureInfodcopy(: in Natural;N: in FloatVector;X: in Integer;IncX: out FloatVector;Y: in Integer);IncY

The two [d|z]geqrf-functions compute a QR-factorization of a real
or complex
` m`-by-

Example 7-1 is one of the longest examples in the running text, but don't be scared as we will explain line-by-line and variable-by-variable what is where and why.

**Example 7-1. Simple native Scilab function
**

subroutine ortho -- Native functions are parameterless implicit none -- Switch into weeny mode :-) * CONSTANTS integer realtype parameter (realtype = 0) -- See Table 7-1 for type association * LOCAL VARIABLES character*6 fname -- The name of the routine as string logical checklhs, checkrhs, cremat, getmat -- Scilab API functions integer topk integer n, m, mattyp integer tausz, worksz, info integer areadr, aimadr, badr, tauadr integer wrkadr, rreadr, rimadr, dumadr * EXTERNAL FUNCTIONS/SUBROUTINES external checklhs, checkrhs, cremat, getmat -- Scilab API functions external error external dcopy, dgeqrf, dorgqr, zgeqrf, zungqr -- LAPACK/BLAS worker subroutines * HEADER include '/site/X11R6/src/scilab/routines/stack.h' -- Scilab API header * TEXT fname = 'ortho' -- Function name (for error messages) topk = top --topis defined instack.hrhs = max(0, rhs) if (.not. checkrhs(fname, 1, 1)) return if (.not. checklhs(fname, 1, 1)) return * fetch input parameters if (.not. getmat(fname, topk, top - rhs + 1, $ mattyp, m, n, areadr, aimadr)) return if (n * m .eq. 0) return -- Quick return on empty matrix tausz = min(m, n) -- Prescribed by man-page worksz = max(1, n) -- ... same here if (mattyp .eq. realtype) then * real case * allocate temporary variables; all are real if (.not. cremat(fname, top + 1, realtype, tausz, 1, $ tauadr, dumadr)) return if (.not. cremat(fname, top + 2, realtype, worksz, 1, $ wrkadr, dumadr)) return if (.not. cremat(fname, top + 3, realtype, m, n, $ badr, dumadr)) return * prepare worker routines' input parameters call dcopy(n * m, stk(areadr), 1, stk(badr), 1) * call worker routines call dgeqrf(m, n, stk(badr), m, stk(tauadr), $ stk(wrkadr), worksz, info) if (info .ne. 0) then -- Any error is considered fatal buf = fname // ' dgeqrf failed' call error(999) return endif call dorgqr(m, n, tausz, stk(badr), m, stk(tauadr), $ stk(wrkadr), worksz, info) if (info .ne. 0) then -- Any error is considered fatal buf = fname // ' dorgqr failed' call error(999) return endif else * complex case; mattyp != realtype * allocate temporary variables, * use two REAL*8 for one COMPLEX*16 if (.not. cremat(fname, top + 1, realtype, 2 * tausz, 1, $ tauadr, dumadr)) return if (.not. cremat(fname, top + 2, realtype, 2 * worksz, 1, $ wrkadr, dumadr)) return if (.not. cremat(fname, top + 3, realtype, 2 * m, 2 * n, $ badr, dumadr)) return * prepare worker routines' input parameters, joining * two REAL*8 arrays into one COMPLEX*16 array call dcopy(n * m, stk(areadr), 1, stk(badr), 2) call dcopy(n * m, stk(aimadr), 1, stk(badr + 1), 2) * call worker routines call zgeqrf(m, n, stk(badr), m, stk(tauadr), $ stk(wrkadr), worksz, info) if (info .ne. 0) then -- Any error is considered fatal buf = fname // ' zgeqrf failed' call error(999) return endif call zungqr(m, n, tausz, stk(badr), m, stk(tauadr), $ stk(wrkadr), worksz, info) if (info .ne. 0) then -- Any error is considered fatal buf = fname // ' zorgqr failed' call error(999) return endif endif * get ready to exit if (lhs .ge. 1) then if (.not. cremat(fname, top, mattyp, m, n, $ rreadr, rimadr)) return if (mattyp .eq. realtype) then call dcopy(m * n, stk(badr), 1, stk(rreadr), 1) else call dcopy(m * n, stk(badr), 2, stk(rreadr), 1) call dcopy(m * n, stk(badr + 1), 2, stk(rimadr), 1) endif endif end

- Check the number of input and output parameters. Here the task is easy as we need one and write one. This line and the next correspond to Step 1.
- Get the addresses as mentioned in Step 2 of the real and imaginary part
of the matrix passed as only parameter to
`ortho`. Note that`getmat`will return`False`if the parameter at the given parameter stack position is not a matrix of numbers. `getmat`is called with the second parameter,`topk`holding the value of the parameter stack pointer when the control flow entered`ortho`. This as well as the function name passed in`fname`is necessary for the cleanup and messaging in case of an error.The only parameter we use is on top of the parameter stack as

`top - rhs + 1`equals`top`in our case.On successful return

`getmat`not only sets the data stack addresses`areadr`, and`aimadr`, but also tells us via`mattyp`whether the matrix is real complex, and via`m`, and`n`how large the matrix is.The following lines directly depend on the sizes passed back, calculating the necessary space for two scratch arrays.

- Allocating space for the temporary variables tau,
work, and b on the data stack is Step 3. tau and work are necessary
because of the LAPACK routines used; b is a copy of a
as the LAPACK routine works with the matrix in place,
i.e. would mangle the input variable a. The
temporaries are accessed the same way parameters are
accessed: through indices into the data stack. These
indices are
`tauadr`,`wrkadr`, and`badr`. Their position on the parameter stack is`top + 1`,`top + 2`, and`top + 3`, respectively. We request a purely real storage for each of the three temporary variables with the third parameter being

. Therefore the index for the imaginary part is a dummy index,`realtype`= 0`dumadr`.The sizes of the vectors or matrices have been computed before.

- There is no "translation" to do in the real case. So Step 4 is easy. The input variable – of which we definitely know that it is real – is simply copied into the scratch space we have allocated on the data stack.
Note how the

*address*of the matrices is passed. The idiom is`stk(`, where)*index*has been obtained through a get*-, or cre*-function. The memnonic "stk" means data stack.*index*- Everything is set up correctly and initialized. We have reached Step 5. The worker routines can take over now.
- In the complex case the allocation of the temporaries
variables requires a bit more thought, although it is
again just Step 3. We know that the LAPACK
routines need the complex vectors/matrices in packed
form. Thus, we allocate
*one*real (DOUBLE PRECISION) vector/matrix of twice the size each time thereby accommodating the the storage requirement of complex (DOUBLE COMPLEX) variables. Otherwise this step proceeds as in the real case. - Due to the different handling of complex variables in Scilab and in LAPACK, Step 4 requires two separate calls to the copy function.
The first line says: "Copy

`m`times`n`elements from the first position in the DOUBLE PRECISION variable`stk(areadr)`taking each entry (3rd parameter, read stride: 1) into the DOUBLE COMPLEX output variable`stk(badr)`filling every other entry (5th parameter, write stride: 2)." The second line does almost the same, but starts off writing at the second element`stk(badr + 1)`, therfore filling the imaginary parts into`stk(badr)`. This corresponds to Step 4.- Again we have reached Step 5; everything is set up correctly and initialized, and the worker routines can take over.
- If there is an output variable, we copy the results into it. Otherwise, we skip the expensive copy operation.
- At this point a purely real result,
`stk(rreadr)`, can simply be copied onto the output parameter,`stk(badr)`. The situation is a bit more complicated for a complex result, as we have to split the

`DOUBLE COMPLEX`result from LAPACK into two`DOUBLE PRECISION`matrices. Here are the crucial lines again:The first line says: "Copy

`m`times`n`elements from the first position in the DOUBLE COMPLEX result`stk(badr)`taking every other entry (3rd parameter, read stride: 2) into the DOUBLE PRECISION output variable`stk(rreadr)`filling each entry (5th parameter, write stride: 1)." The second line does almost the same, but starts off at the second element`stk(badr + 1)`, therfore copying the imaginary parts into`stk(rimadr)`. This way we are merging Step 6, and Step 7 into one.

call dcopy(n * m, stk(areadr), 1, stk(badr), 2) call dcopy(n * m, stk(aimadr), 1, stk(badr + 1), 2)

call dcopy(m * n, stk(badr), 2, stk(rreadr), 1) call dcopy(m * n, stk(badr + 1), 2, stk(rimadr), 1)

FIXME: write it!