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\).

parametric real type T#

Below, the symbol T refers to a parametric real type that may be Float32 (single precision) or Float64 (double precision). Calable functions as described are with T as Float64, but variants (with the additional suffix _s, e.g., ir_initialize_s) are available with T as Float32.

callable functions#

    function ir_initialize(data, control, status)

Set default control values and initialize private data

Parameters:

data

holds private internal data

control

is a structure containing control information (see ir_control_type)

status

is a scalar variable of type Int32 that gives the exit status from the package. Possible values are (currently):

  • 0

    The initialization was successful.

    function ir_information(data, inform, status)

Provides output information

Parameters:

data

holds private internal data

inform

is a structure containing output information (see ir_inform_type)

status

is a scalar variable of type Int32 that gives the exit status from the package. Possible values are (currently):

  • 0

    The values were recorded successfully

    function ir_terminate(data, control, inform)

Deallocate all internal private storage

Parameters:

data

holds private internal data

control

is a structure containing control information (see ir_control_type)

inform

is a structure containing output information (see ir_inform_type)

available structures#

ir_control_type structure#

    struct ir_control_type{T}
      f_indexing::Bool
      error::Int32
      out::Int32
      print_level::Int32
      itref_max::Int32
      acceptable_residual_relative::T
      acceptable_residual_absolute::T
      required_residual_relative::T
      record_residuals::Bool
      space_critical::Bool
      deallocate_error_fatal::Bool
      prefix::NTuple{31,Cchar}

detailed documentation#

control derived type as a Julia structure

components#

Bool f_indexing

use C or Fortran sparse matrix indexing

Int32 error

unit for error messages

Int32 out

unit for monitor output

Int32 print_level

controls level of diagnostic output

Int32 itref_max

maximum number of iterative refinements allowed

T acceptable_residual_relative

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

T acceptable_residual_absolute

see acceptable_residual_relative

T 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

NTuple{31,Cchar} prefix

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#

    struct ir_inform_type{T}
      status::Int32
      alloc_status::Int32
      bad_alloc::NTuple{81,Cchar}
      norm_initial_residual::T
      norm_final_residual::T

detailed documentation#

inform derived type as a Julia structure

components#

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

Int32 alloc_status

the status of the last attempted allocation/deallocation.

NTuple{81,Cchar} bad_alloc

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

T norm_initial_residual

the infinity norm of the initial residual

T norm_final_residual

the infinity norm of the final residual