GALAHAD ROOTS package

purpose

The roots package uses classical formulae together with Newton’s method to find all the real roots of a real polynomial.

Currently only the options and inform dictionaries are exposed; these are provided and used by other GALAHAD packages with Python interfaces. Please contact us if you would like full functionality!

See Section 4 of $GALAHAD/doc/roots.pdf for additional details.

method

Littlewood and Ferrari’s algorithms are used to find estimates of the real roots of cubic and quartic polynomials, respectively; a stabilized version of the well-known formula is used in the quadratic case. Newton’s method and/or methods based on the companion matrix are used to further refine the computed roots if necessary. Madsen and Reid’s (1975) method is used for polynomials whose degree exceeds four.

reference

The basic method is that given by

K. Madsen and J. K. Reid, ``FORTRAN Subroutines for Finding Polynomial Zeros’’. Technical Report A.E.R.E. R.7986, Computer Science and System Division, A.E.R.E. Harwell, Oxfordshire (1975)

callable functions

overview of functions provided

// namespaces

namespace conf;

// typedefs

typedef float spc_;
typedef double rpc_;
typedef int ipc_;

// structs

struct roots_control_type;
struct roots_inform_type;

// global functions

void roots_initialize(
    void **data,
    struct roots_control_type* control,
    ipc_ *status
);

void roots_information(
    void **data,
    struct roots_inform_type* inform,
    ipc_ *status
);

void roots_terminate(
    void **data,
    struct roots_control_type* control,
    struct roots_inform_type* inform
);

typedefs

typedef float spc_

spc_ is real single precision

typedef double rpc_

rpc_ is the real working precision used, but may be changed to float by defining the preprocessor variable REAL_32 or (if supported) to __real128 using the variable REAL_128.

typedef int ipc_

ipc_ is the default integer word length used, but may be changed to int64_t by defining the preprocessor variable INTEGER_64.

function and structure names

The function and structure names described below are appropriate for the default real working precision (double) and integer word length (int32_t). To use the functions and structures with different precisions and integer word lengths, an additional suffix must be added to their names (and the arguments set accordingly). The appropriate suffices are:

_s for single precision (float) reals and standard 32-bit (int32_t) integers;

_q for quadruple precision (__real128) reals (if supported) and standard 32-bit (int32_t) integers;

_64 for standard precision (double) reals and 64-bit (int64_t) integers;

_s_64 for single precision (float) reals and 64-bit (int64_t) integers; and

_q_64 for quadruple precision (__real128) reals (if supported) and 64-bit (int64_t) integers.

Thus a call to roots_initialize below will instead be

void roots_initialize_s_64(void **data, struct roots_control_type_s_64* control,
                         int64_t *status)

if single precision (float) reals and 64-bit (int64_t) integers are required. Thus it is possible to call functions for this package with more that one precision and/or integer word length at same time. An example is provided for the package expo, and the obvious modifications apply equally here.

function calls

void roots_initialize(
    void **data,
    struct roots_control_type* control,
    ipc_ *status
)

Set default control values and initialize private data

Parameters:

data	holds private internal data
control	is a struct containing control information (see roots_control_type)
status	is a scalar variable of type ipc_, that gives the exit status from the package. Possible values are (currently): 0 The initialization was successful.

void roots_information(
    void **data,
    struct roots_inform_type* inform,
    ipc_ *status
)

Provides output information

Parameters:

data	holds private internal data
inform	is a struct containing output information (see roots_inform_type)
status	is a scalar variable of type ipc_, that gives the exit status from the package. Possible values are (currently): 0 The values were recorded successfully

void roots_terminate(
    void **data,
    struct roots_control_type* control,
    struct roots_inform_type* inform
)

Deallocate all internal private storage

Parameters:

data	holds private internal data
control	is a struct containing control information (see roots_control_type)
inform	is a struct containing output information (see roots_inform_type)

available structures

roots_control_type structure

#include <galahad_roots.h>

struct roots_control_type {
    // fields

    bool f_indexing;
    ipc_ error;
    ipc_ out;
    ipc_ print_level;
    rpc_ tol;
    rpc_ zero_coef;
    rpc_ zero_f;
    bool space_critical;
    bool deallocate_error_fatal;
    char prefix[31];
};

detailed documentation

control derived type as a C struct

components

bool f_indexing

use C or Fortran sparse matrix indexing

ipc_ error

error and warning diagnostics occur on stream error

ipc_ out

general output occurs on stream out

ipc_ print_level

the level of output required is specified by print_level

rpc_ tol

the required accuracy of the roots

rpc_ zero_coef

any coefficient smaller in absolute value than zero_coef will be regarde to be zero

rpc_ zero_f

any value of the polynomial smaller in absolute value than zero_f will be regarded as giving a root

bool space_critical

if .space_critical true, every effort will be made to use as little space as possible. This may result in longer computation time

bool deallocate_error_fatal

if .deallocate_error_fatal is true, any array/pointer deallocation error will terminate execution. Otherwise, computation will continue

char prefix[31]

all output lines will be prefixed by .prefix(2:LEN(TRIM(.prefix))-1) where .prefix contains the required string enclosed in quotes, e.g. “string” or ‘string’

roots_inform_type structure

#include <galahad_roots.h>

struct roots_inform_type {
    // fields

    ipc_ status;
    ipc_ alloc_status;
    char bad_alloc[81];
};

detailed documentation

inform derived type as a C struct

components

ipc_ status

return status. Possible values are:

• 0

  the call was successful.

• -1

  An allocation error occurred. A message indicating the offending array is written on unit control.error, and the returned allocation status and a string containing the name of the offending array are held in inform.alloc_status and inform.bad_alloc respectively.

• -2

  A deallocation error occurred. A message indicating the offending array is written on unit control.error and the returned allocation status and a string containing the name of the offending array are held in inform.alloc_status and inform.bad_alloc respectively.

• -3

  Either the specified degree of the polynomial in degree is less than 0, or the declared dimension of the array roots is smaller than the specified degree.

ipc_ alloc_status

the status of the last attempted allocation/deallocation

char bad_alloc[81]

the name of the array for which an allocation/deallocation error occurred

Index