# Interval Linear Solvers Code Examples

This section presents code examples of using the routines described in “Interval Linear Solvers”. These routines are intended for computing enclosures and estimates of the solution sets to interval linear systems of equations as well as for checking properties of interval matrices and their inversion.

## Example. Interval Gauss-Seidel method

Given an interval system of linear algebraic equations, interval Gauss-Seidel method (implemented as ?gegss routine) is often applied for enclosing a desired portion of the solution set that is bounded by a prescribed interval box.

Consider the following interval linear system of equations

proposed first by E. Hansen (see Hansen92). Does its solution set intersect the interval box

The following sample program answers the above question.

`PROGRAM DIGEGSS_EXAMPLE`
`!`
`! Example program enclosing the solution set to a square interval`
`! linear system by interval Gauss-Seidel iterative method`
`!`
`!----------------------------------------------------------------------!`
`USE INTERVAL_ARITHMETIC`
`IMPLICIT NONE`
`!----------------------------------------------------------------------!`
`INTEGER, PARAMETER :: DIM = 2`
`INTEGER :: NRHS, LDA, LDB, NITS,`
`INFO, I, J`
`REAL(8) :: EPSILON`
`TYPE(D_INTERVAL) :: A(DIM,DIM), B(DIM,1),`
`ENCL(DIM,1)`
`CHARACTER(1) :: TRANS`
`!----------------------------------------------------------------------!`
`PRINT 300`
`!----------------------------------------------------------------------!`
`! !`
`! Initializing the input data - !`
`! !`
`TRANS = 'N'`
`NRHS = 2`
`A(1,1) = DINTERVAL(2.,3.); A(1,2) = DINTERVAL(0.,1.);`
`A(2,1) = DINTERVAL(0.,1.); A(2,2) = DINTERVAL(2.,3.);`
`LDA = 2`
`B(1,1) = DINTERVAL(0.,120.); B(2,1) = DINTERVAL(60.,240.);`
`LDB = 2`
`EPSILON = 1.D-6`
`NITS = 20`
`!----------------------------------------------------------------------!`
`!`
`! Assigning the bounding box for the solution set -`
`DO I = 1, DIM`
`ENCL(I,1) = DINTERVAL(0.,200.)`
`END DO`
`!----------------------------------------------------------------------!`
`CALL DIGEGSS(TRANS,`
`DIM, NRHS, A, LDA, B, LDB, ENCL, EPSILON, NITS, INFO )`
`!----------------------------------------------------------------------!`
`!`
`! Outputting the solution`
`IF( INFO /= 0 ) THEN`
`PRINT 400`
`ELSE`
`PRINT 600`
`DO I = 1, DIM`
`PRINT *, '[', B(I,1), ']'`
`END DO`
`END IF`
`!----------------------------------------------------------------------!`
`300 FORMAT (/,' **** SOLVING INTERVAL LINEAR SYSTEM **** ',/, &`
`' by interval Gauss-Seidel method ')`
`400 FORMAT (/,' The interval Gauss-Seidel method fails. ')`
`600 FORMAT (/,' Outer interval estimate of the solution set:',/)`
`!----------------------------------------------------------------------!`
`END PROGRAM DIGEGSS_EXAMPLE`

Assigning double-precision intervals to the entries of the matrix A and right-hand side vector B is carried out by DINTERVAL function that turns two real numbers into the interval having these reals as endpoints. Running the above code produces the answer

`**** SOLVING INTERVAL LINEAR SYSTEM****`
`by interval Gauss-Seidel method`
`Outer interval estimate of the solution set:`
`[ 0.000000000000000E+000 60.0000000000000 ]`
`[ 0.000000000000000E+000 120.000000000000 ]`

One can make sure that the resulting box really encloses the required portion of the solution set after having a look at the corresponding graph from the paper Hansen92. Moreover, it is even the tightest possible enclosure.

## Example. Hansen-Bliek-Rohn procedure

The following Fortran 90 program illustrates the use of digehbs routine implementing “semiinterval” Hansen-Bliek-Rohn procedure for outer interval estimation of the solution sets to interval linear systems.

`PROGRAM DIGEHBS_EXAMPLE`
`!`
`! Example program for enclosing the solution set to square interval`
`! interval system of equations by Hansen-Bliek-Rohn procedure`
`!`
`!----------------------------------------------------------------------!`
`USE INTERVAL_ARITHMETIC`
`IMPLICIT NONE`
`!----------------------------------------------------------------------!`
`INTEGER, PARAMETER :: DIM = 2`
`INTEGER :: LDA, LDB, INFO, I, J`
`TYPE(D_INTERVAL), ALLOCATABLE ::`
`A(:,:), B(:)`
`CHARACTER(1) :: TRANS`
`!----------------------------------------------------------------------!`
`PRINT 300`
`!----------------------------------------------------------------------!`
`!`
`! Initializing the input data -`
`!`
`TRANS = 'N'`
`ALLOCATE( A(DIM,DIM), B(DIM) )`
`A(1,1) = DINTERVAL(2.,4.); A(1,2) = DINTERVAL(-2.,1.)`
`A(2,1) = DINTERVAL(-1.,2.); A(2,2) = DINTERVAL(2.,2.)`
`LDA = 2`
`B(1) = DINTERVAL(0.,2.); B(2) = DINTERVAL(0.,2.)`
`LDB = 2`
`!----------------------------------------------------------------------!`
`CALL DIGEHBS( TRANS, DIM, A, LDA,`
`B, LDB, INFO )`
`!----------------------------------------------------------------------!`
`IF( INFO /= 0 ) THEN`
`PRINT 400`
`ELSE`
`PRINT 600`
`DO I = 1, DIM`
`PRINT *, I, ') [', B(I), ']'`
`END DO`
`END IF`
`!----------------------------------------------------------------------!`
`DEALLOCATE( A, B )`
`!----------------------------------------------------------------------!`
`300 FORMAT (/,' **** SOLVING INTERVAL LINEAR SYSTEM ****',/, &`
`' by Hansen-Bliek-Rohn procedure ',/)`
`400 FORMAT (/,' The matrix of the system is not an H-matrix,',/, &`
`' Hansen-Bliek-Rohn procedure fails. ',/)`
`600 FORMAT (/,' Enclosure of the solution set: ',/)`
`!----------------------------------------------------------------------!`
`END PROGRAM DIGEHBS_EXAMPLE`

However, the output of the program looks like

`**** SOLVING INTERVAL LINEAR SYSTEM****`
`by Hansen-Bliek-Rohn procedure`
`The matrix of the system is not an H-matrix,`
`Hansen-Bliek-Rohn procedure fails.`
``

This result is obtained because the program is applied to the interval linear system

where the interval matrix is not an H-matrix (that is, it does not have diagonal dominance).

However, preconditioning by digemip routine helps to resolve the problem. The following modified program, which incorporates preliminary preconditioning of the interval linear system under solution, makes the matrix diagonally dominant and produces an acceptable answer to the problem.

`PROGRAM DIGEMIP_DIGEHBS_EXAMPLE`
`!`
`! Example program for enclosing the solution set to square interval`
`! interval system of equations by Hansen-Bliek-Rohn procedure`
`!`
`!----------------------------------------------------------------------!`
`USE INTERVAL_ARITHMETIC`
`IMPLICIT NONE`
`!----------------------------------------------------------------------!`
`INTEGER, PARAMETER :: DIM = 2, NRHS = 1`
`INTEGER :: LDA, LDB, INFO, I, J`
`TYPE(D_INTERVAL), ALLOCATABLE ::`
`A(:,:), B(:)`
`CHARACTER(1) :: TRANS`
`!----------------------------------------------------------------------!`
`PRINT 300`
`!----------------------------------------------------------------------!`
`!`
`! Initializing the input data -`
`!`
`TRANS = 'N'`
`ALLOCATE( A(DIM,DIM), B(DIM) )`
`A(1,1) = DINTERVAL(2.,4.); A(1,2) = DINTERVAL(-2.,1.)`
`A(2,1) = DINTERVAL(-1.,2.); A(2,2) = DINTERVAL(2.,2.)`
`LDA = 2`
`B(1) = DINTERVAL(0.,2.); B(2) = DINTERVAL(0.,2.)`
`LDB = 2`
`!----------------------------------------------------------------------!`
`CALL DIGEMIP( DIM, NRHS, A, LDA,`
`B, LDB, INFO )`
`CALL DIGEHBS( TRANS, DIM, A, LDA,`
`B, LDB, INFO )`
`!----------------------------------------------------------------------!`
`IF( INFO /= 0 ) THEN`
`PRINT 400`
`ELSE`
`PRINT 600`
`DO I = 1, DIM`
`PRINT *, I, ') [', B(I), ']'`
`END DO`
`END IF`
`DEALLOCATE( A, B )`
`!----------------------------------------------------------------------!`
`300 FORMAT (/,' **** SOLVING INTERVAL LINEAR SYSTEM ****',/, &`
`' by Hansen-Bliek-Rohn procedure ',/)`
`400 FORMAT (/,' The matrix of the system is not an H-matrix,',/, &`
`' Hansen-Bliek-Rohn procedure fails. ',/)`
`600 FORMAT (/,' Enclosure of the solution set: ',/)`
`!----------------------------------------------------------------------!`
`END PROGRAM DIGEMIP_DIGEHBS_EXAMPLE`

This time, the output of the program is

`**** SOLVING INTERVAL LINEAR SYSTEM****`
`by Hansen-Bliek-Rohn procedure`
`Enclosure of the solution set:`
`1 ) [ -4.23529411764708 10.7058823529412 ]`
`2 ) [ -6.70588235294119 10.8235294117647 ]`

(the last digits may change for various computer architectures).

## Example. Computing enclosure for inverse interval matrix

Given an interval 2 × 2 matrix

the following Fortran 90 code computes an enclosure of its inverse interval matrix:

`PROGRAM SIGESZI_EXAMPLE`
`!`
`!Example program inverting an interval matrix by Sczulz iterative procedure`
`!`
`!----------------------------------------------------------------------!`
`USE INTERVAL_ARITHMETIC`
`IMPLICIT NONE`
`!----------------------------------------------------------------------!`
`INTEGER, PARAMETER :: DIM = 2, LDA = 2`
`INTEGER :: INFO, I, J`
`TYPE(S_INTERVAL), ALLOCATABLE ::`
`A(:,:)`
`!----------------------------------------------------------------------!`
`PRINT 300`
`!----------------------------------------------------------------------!`
`!`
`! Initislizing the input data -`
`!`
`ALLOCATE( A(LDA,DIM) )`
`A(1,1) = SINTERVAL(3.,3.); A(1,2) = SINTERVAL(0.,1.)`
`A(2,1) = SINTERVAL(1.,2.); A(2,2) = SINTERVAL(2.,3.)`
`!----------------------------------------------------------------------!`
`CALL SIGESZI ( DIM, A, LDA, INFO )`
`!----------------------------------------------------------------------!`
`PRINT 600`
`DO I = 1, DIM`
`PRINT *, ('[', A(I,J), ']', J = 1,`
`DIM )`
`END DO`
`DEALLOCATE( A )`
`!----------------------------------------------------------------------!`
`300 FORMAT (/,' **** INVERTING INTERVAL MATRIX ****',/, &`
`' by interval Schulz method ' )`
`400 FORMAT (/,' Schulz inversion procedure failed. ',/)`
`600 FORMAT (/,' Enclosure of the inverse matrix ',/)`
`!----------------------------------------------------------------------!`
`END PROGRAM SIGESZI_EXAMPLE`

The output listing (with small variations depending on the architecture) looks as follows:

`**** INVERTING INTERVAL MATRIX ****`
`by interval Schulz method`
`Enclosure of the inverse matrix`
`[ 0.2407409 0.5000001 ][ -0.2500000 0.1018518 ]`
`[ -0.5000000 5.5555239E-02 ][ 0.1388889 0.7500001 ]`

At the same time, if we widen the (1,1) entry of the matrix to the interval [2, 3], the sigeszi procedure fails to compute a finite enclosure of the inverse to the new interval matrix

Nevertheless, the interval linear system with such matrix can be successfully solved by specialized routines, for example, by interval Gauss method or interval Gauss-Seidel method (see Example “Interval Gauss-Seidel method”).

Parent topic: Code Examples