Previous External Development Guide: Adding System Routines Next

Example: A Complete Numerical Routine Example (FZ_ROOTS2)

The following is a complete implementation of the IDL system function FZ_ROOTS, used to find the roots of an m-degree complex polynomial, using Laguerre's method. The result is an m-element complex vector. We call this version FZ_ROOTS2 to avoid a name clash with the real routine. FZ_ROOTS2 has an additional keyword, TC_INPUT, that is not present in the real routine.

FZ_ROOTS2 uses the routine zroots(), described in section 9.5 of Numerical Recipes in C: The Art of Scientific Computing (Second Edition), published by Cambridge University Press: 

void zroots(fcomplex a[], int m, fcomplex roots[], int polish)

Quoting from the referenced book:

Given the degree m and the m+1 complex coefficients a[0..m] of the polynomial

,

this routine successively calls laguer and finds all m complex roots in roots[1..m]. The boolean variable polish should be input as true (1) if polishing (also by Laguerre's method) is desired, false (0) if the roots will be subsequently polished by other means.

FZ_ROOTS2 will support both single and double precision complex values as well as give the caller control over the error tolerance, which is hard wired into the Numerical Recipes code as a C preprocessor constant named EPS. In order to support these requirements, we have copied the zroots() function given in the book and altered it to support both data types and make EPS a user specified parameter, giving two functions: 

void zroots_f(fcomplex a[], int m, fcomplex roots[], int polish,  
             float eps);  
  
void zroots_d(dcomplex a[], int m, dcomplex roots[], int polish,  
             double eps);

Note that fcomplex and dcomplex are Numerical Recipes defined types that happen to have the same definition as the IDL types IDL_COMPLEX and IDL_DCOMPLEX, a convenient fact that eliminates some type conversion issues.

The definition of FZ_ROOTS2 from the IDL user perspective is:

Calling Sequence

Result = FZ_ROOTS2(C)

Arguments

C

A vector of length m+1 containing the coefficients of the polynomial, in ascending order.

Keywords

DOUBLE

FZ_ROOTS2 normally uses the type of C to determine the type of the computation. If DOUBLE is specified, it overrides this default. Setting DOUBLE to a non-zero value causes the computation type and the result to be double precision complex. Setting it to zero forces single precision complex.

EPS

The desired fractional accuracy. The default value is 2.0 ´ 10-6.

NO_POLISH

Set this keyword to suppress the usual polishing of the roots by Laguerre's method.

TC_INPUT

If present, TC_INPUT specifies a named variable that will be assigned the input value C, with its type converted to the type of the result.

Example

The following figure gives the code for fzroots2.c,. This is ANSI C code that implements FZ_ROOTS2. The line numbers are not part of the code and are present to make the discussion easier to follow. Each line of this routine is discussed below.

Table 15-2: fzroots2.c 

Table 15-2: fzroots2.c 
C 
1  
2  
3  
4  
5  
6  
7  
8  
9  
10  
11  
12  
13  
14  
15  
16  
17  
18  
19  
20  
21  
22  
23  
24  
25  
26  
27  
28  
29  
30  
31  
32  
33  
34  
35  
36  
37  
38  
39  
40  
41  
42  
43  
44  
45  
46  

#include <stdio.h>  
#include <stdarg.h>  
#include "idl_export.h"  
#include <nr/nr.h>  
   
IDL_VPTR fzroots2(int argc, IDL_VPTR *argv, char *argk)  
{  
  typedef struct {  
    IDL_KW_RESULT_FIRST_FIELD;	 /* Must be first entry in this 
structure */  
    int force_type;  
    IDL_LONG do_double;  
    double eps;  
    IDL_LONG no_polish;  
    IDL_VPTR tc_input;  
  } KW_RESULT;  
  static IDL_KW_PAR kw_pars[] = {  
    {"DOUBLE", IDL_TYP_LONG, 1, 0,  
     IDL_KW_OFFSETOF(force_type), IDL_KW_OFFSETOF(do_double) },  
    { "EPS", IDL_TYP_DOUBLE, 1, 0, 0, IDL_KW_OFFSETOF(eps) },  
    { "NO_POLISH", IDL_TYP_LONG, 1, IDL_KW_ZERO,  
      0, IDL_KW_OFFSETOF(no_polish) },  
    { "TC_INPUT", 0, 1, IDL_KW_OUT|IDL_KW_ZERO,  
      0, IDL_KW_OFFSETOF(tc_input) },   
    { NULL }  
  };  
   
  KW_RESULT kw;  
  IDL_VPTR result;  
  IDL_VPTR c_raw;  
  IDL_VPTR c_tc;  
  IDL_MEMINT m;  
  void *outdata;  
  IDL_ARRAY_DIM dim;  
  int rtype;  
  static IDL_ALLTYPES zero;  
   
  kw.eps = 2.0e-6;  
  (void) IDL_KWProcessByOffset(argc, argv, argk, 
kw_pars,&c_raw,1,&kw);  
   
  IDL_ENSURE_ARRAY(c_raw);  
  IDL_ENSURE_SIMPLE(c_raw);  
  if (c_raw->value.arr->n_dim != 1)   
  IDL_Message(IDL_M_NAMED_GENERIC, IDL_MSG_LONGJMP,  
              "Input argument must be a column vector.");
C 
47  
48  
49  
50  
51  
52  
53  
54  
55  
56  
57  
58  
59  
60  
61  
62  
63  
64  
65  
66  
67  
68  
69  
70  
71  
72  
73  
74  
75  
76  
77  
78  
79  
80  
81  
82  
83  

 m = c_raw->value.arr->dim[0];   
if (--m <= 0)  
    IDL_Message(IDL_M_NAMED_GENERIC, IDL_MSG_LONGJMP,   
                "Input array does not have enough elements");  
if (kw.tc_input)  
    IDL_StoreScalar(kw.tc_input, IDL_TYP_LONG, &zero);  
   
  if (kw.force_type) {   
    rtype = kw.do_double ? IDL_TYP_DCOMPLEX : IDL_TYP_COMPLEX;  
  } else {  
    rtype = ((c_raw->type == IDL_TYP_DOUBLE)  
             || (c_raw->type == IDL_TYP_DCOMPLEX))  
      ? IDL_TYP_DCOMPLEX : IDL_TYP_COMPLEX;  
  }   
  dim[0] = m;   
  outdata = (void *)  
    IDL_MakeTempArray(rtype,1,dim,IDL_ARR_INI_NOP,&result);   
   
  if (c_raw->type == rtype) {   
    c_tc = c_raw;   
  } else {   
    c_tc = IDL_BasicTypeConversion(1, &c_raw, rtype);   
  }   
   
  if (rtype == IDL_TYP_COMPLEX) {   
    zroots_f((fcomplex *) c_tc->value.arr->data, m,   
             ((fcomplex *)outdata)-1,!kw.no_polish,(float) kw.eps);   
  } else {   
    zroots_d((dcomplex *) c_tc->value.arr->data, m,   
             ((dcomplex *) outdata) - 1, !kw.no_polish, kw.eps);   
  }   
   
  if (kw.tc_input) IDL_VarCopy(c_tc, kw.tc_input);   
  else if (c_raw != c_tc) IDL_Deltmp(c_tc);   
   
  IDL_KW_FREE;  
  return result;   
}

4

nr.h is the header file provided with Numerical Recipes in C code.

6

FZROOTS2 has the usual three standard arguments.

10

kw.force_type will be TRUE if the user specifies the DOUBLE keyword. In this case, the value of the DOUBLE keyword will determine the result type without regard for the type of the input argument.

If the user specifies DOUBLE, a zero value forces a single precision complex result and non-zero forces double precision complex.

12

The value of the EPS keyword.

13

The value of the NO_POLISH keyword.

15

The value of the TC_INPUT keyword.

16

This array defines the keywords accepted by FZ_ROOTS2.

17

Since setting DOUBLE to 0 has a different meaning than not specifying the keyword at all, kw.force_type is used to detect the fact that the keyword is set independent of its value.

20

The EPS keyword allows the user to specify the kw.eps tolerance parameter. kw.eps is specified as double precision to avoid losing accuracy for double precision computations—it will be converted to single precision if necessary. The default value for this keyword is non-zero, so no zeroing is specified here. If the user includes the EPS keyword, the value will be placed in kw.eps, otherwise kw.eps will not be changed.

21

This keyword lets the user suppress the usual polishing performed by zroots(). Since specifying a value of 0 is equivalent to not specifying the keyword at all, IDL_KW_ZERO is used to initialize the variable.

23

If present, TC_INPUT is an output keyword that will have the type converted value of the input argument stored in it. By specifying IDL_KW_OUT and IDL_KW_ZERO, we ensure that TC_INPUT is either zero or a pointer to a valid IDL variable.

28

The results of keyword processing will all be written to this variable by IDL_KWProcessByOffset().

29

This variable will receive the function result.

30

The input argument prior to any type conversion.

31

The type converted input variable. If the input variable is already of the correct type, this will be the same as c_raw, otherwise it will be different.

32

The value of m to be passed to zroots().

33

Pointer to the data area of the result variable. We declare it as (void *) so that it can point to data of any type.

34

Used to specify dimensions of the result. This will always be a vector of m elements.

35

IDL type code for result variable.

36

Used by IDL_StoreScalar() to type check the TC_INPUT keyword. It is declared as static to ensure it is initialized to zero.

38

Set the default EPS value before doing keyword processing. If the user specifies EPS, the supplied value will override this. Otherwise, this value will still be in kw.eps and will be passed to zroots() unaltered.

39

Perform keyword processing.

42-43

Ensure that the input argument is an array, and is one of the basic types (not a file variable or structure).

44-46

The input variable must be a vector, and therefore should have only a single dimension.

47-50

Ensure that the input variable is long enough for m to be non-zero. m is one less than the number of elements in the input vector, so this is equivalent to saying that the input must have at least 2 elements.

51

If the TC_INPUT keyword was present, use IDL_StoreScalar() to make sure the named variable specified can receive the converted input value. A nice side effect of this operation is that any dynamic memory currently being used by this variable will be freed now instead of later after we have allocated other dynamic memory. This freed memory might be immediately reusable if it is large enough, which would reduce memory fragmentation and lower overall memory requirements.

54

If the user specified the DOUBLE keyword, it is used to control the resulting type, otherwise the input argument type is used to decide.

55

The DOUBLE keyword was specified. If it is non-zero, use IDL_TYP_DCOMPLEX, otherwise IDL_TYP_COMPLEX.

56-60

Use the input type to decide the result type. If the input is IDL_TYP_DOUBLE or IDL_TYP_DCOMPLEX, use IDL_TYP_DCOMPLEX, otherwise IDL_TYP_COMPLEX.

61-63

Create the output variable that will be passed back as the result of FZ_ROOTS2.

65-69

If necessary, convert the input argument to the result type. This is done after creation of the output variable because it is likely to have a short lifetime. If it does get freed at the end of this routine, it won't cause memory fragmentation by leaving a hole in the process virtual memory.

71

The version of zroots() to call depends on the data type of the result.

72-73

Single precision complex. Note that the outdata pointer is decremented by one element. This compensates for the fact that the Numerical Recipe routine will index it from [1..m] rather than [0..m-1] as is the usual C convention. Also, kw.eps is cast to single precision.

74-76

Double precision complex case.

79

If the user specified the TC_INPUT keyword, copy the type converted input into the keyword variable. Since VarCopy() frees its source variable if it is a temporary variable, we are relieved of the usual responsibility to call IDL_Deltmp() if c_tc contains a temporary variable created on line 66.

80

The user didn't specify the TC_INPUT keyword. In this case, if we allocated c_tc on line 66, we must free it before returning.

82

Free any resources allocated by keyword processing.

83

Return the result.

  IDL Online Help (March 06, 2007)