GALAHAD IR package

purpose

Given a sparse symmetric \(n \times n\) matrix \(A = a_{ij}\) and the factorization of \(A\) found by the GALAHAD package SLS, the ir package solves the system of linear equations \(A x = b\) using iterative refinement.

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/ir.pdf for additional details.

method

Iterative refinement proceeds as follows. First obtain the floating-point solution to \(A x = b\) using the factors of \(A\). Then iterate until either the desired residual accuracy (or the iteration limit is reached) as follows: evaluate the residual \(r = b - A x\), find the floating-point solution \(\delta x\) to \(A \delta x = r\), and replace \(x\) by \(x + \delta x\).

callable functions

overview of functions provided

// namespaces

namespace conf;

// typedefs

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

// structs

struct ir_control_type;
struct ir_inform_type;

// global functions

void ir_initialize(void **data, struct ir_control_type* control, ipc_ *status);
void ir_information(void **data, struct ir_inform_type* inform, ipc_ *status);

void ir_terminate(
    void **data,
    struct ir_control_type* control,
    struct ir_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 ir_initialize below will instead be

void ir_initialize_s_64(void **data, struct ir_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 ir_initialize(void **data, struct ir_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 ir_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 ir_information(void **data, struct ir_inform_type* inform, ipc_ *status)

Provides output information

Parameters:

data	holds private internal data
inform	is a struct containing output information (see ir_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 ir_terminate(
    void **data,
    struct ir_control_type* control,
    struct ir_inform_type* inform
)

Deallocate all internal private storage

Parameters:

data	holds private internal data
control	is a struct containing control information (see ir_control_type)
inform	is a struct containing output information (see ir_inform_type)

available structures

ir_control_type structure

#include <galahad_ir.h>

struct ir_control_type {
    // fields

    bool f_indexing;
    ipc_ error;
    ipc_ out;
    ipc_ print_level;
    ipc_ itref_max;
    rpc_ acceptable_residual_relative;
    rpc_ acceptable_residual_absolute;
    rpc_ required_residual_relative;
    bool record_residuals;
    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

unit for error messages

ipc_ out

unit for monitor output

ipc_ print_level

controls level of diagnostic output

ipc_ itref_max

maximum number of iterative refinements allowed

rpc_ acceptable_residual_relative

refinement will cease as soon as the residual \(\|Ax-b\|\) falls below max( acceptable_residual_relative \* \(\|b\|\), acceptable_residual_absolute )

rpc_ acceptable_residual_absolute

see acceptable_residual_relative

rpc_ required_residual_relative

refinement will be judged to have failed if the residual \(\|Ax-b\| \geq\) required_residual_relative \* \(\|b\|\). No checking if required_residual_relative < 0

bool record_residuals

record the initial and final residual

bool space_critical

if space is critical, ensure allocated arrays are no bigger than needed

bool deallocate_error_fatal

exit if any deallocation fails

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’

ir_inform_type structure

#include <galahad_ir.h>

struct ir_inform_type {
    // fields

    ipc_ status;
    ipc_ alloc_status;
    char bad_alloc[81];
    rpc_ norm_initial_residual;
    rpc_ norm_final_residual;
};

detailed documentation

inform derived type as a C struct

components

ipc_ status

the return status. Possible values are:

• 0

  the solution has been found.

• -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.

• -11

  Iterative refinement has not reduced the relative residual by more than control.required_relative_residual.

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.

rpc_ norm_initial_residual

the infinity norm of the initial residual

rpc_ norm_final_residual

the infinity norm of the final residual

Index