Integer type:  int32  int64  nag_int  show int32  show int32  show int64  show int64  show nag_int  show nag_int

PDF version (NAG web site, 64-bit version, 64-bit version)
Chapter Contents
Chapter Introduction
NAG Toolbox

# NAG Toolbox: nag_linsys_real_toeplitz_update (f04mf)

## Purpose

nag_linsys_real_toeplitz_update (f04mf) updates the solution of the equations $Tx=b$, where $T$ is a real symmetric positive definite Toeplitz matrix.

## Syntax

[x, p, work, ifail] = f04mf(t, b, x, work, 'n', n)
[x, p, work, ifail] = nag_linsys_real_toeplitz_update(t, b, x, work, 'n', n)

## Description

nag_linsys_real_toeplitz_update (f04mf) solves the equations
 $Tnxn=bn,$
where ${T}_{n}$ is the $n$ by $n$ symmetric positive definite Toeplitz matrix
 $Tn= τ0 τ1 τ2 … τn-1 τ1 τ0 τ1 … τn-2 τ2 τ1 τ0 … τn-3 . . . . τn-1 τn-2 τn-3 … τ0$
and ${b}_{n}$ is the $n$-element vector ${b}_{n}={\left({\beta }_{1}{\beta }_{2}\dots {\beta }_{n}\right)}^{\mathrm{T}}$, given the solution of the equations
 $Tn-1xn-1=bn-1.$
This function will normally be used to successively solve the equations
 $Tkxk=bk, k= 1,2,…,n.$
If it is desired to solve the equations for a single value of $n$, then function nag_linsys_real_toeplitz_solve (f04ff) may be called. This function uses the method of Levinson (see Levinson (1947) and Golub and Van Loan (1996)).

## References

Bunch J R (1985) Stability of methods for solving Toeplitz systems of equations SIAM J. Sci. Statist. Comput. 6 349–364
Bunch J R (1987) The weak and strong stability of algorithms in numerical linear algebra Linear Algebra Appl. 88/89 49–66
Cybenko G (1980) The numerical stability of the Levinson–Durbin algorithm for Toeplitz systems of equations SIAM J. Sci. Statist. Comput. 1 303–319
Golub G H and Van Loan C F (1996) Matrix Computations (3rd Edition) Johns Hopkins University Press, Baltimore
Levinson N (1947) The Weiner RMS error criterion in filter design and prediction J. Math. Phys. 25 261–278

## Parameters

### Compulsory Input Parameters

1:     $\mathrm{t}\left(:\right)$ – double array
The dimension of the array t must be at least $\mathrm{max}\phantom{\rule{0.125em}{0ex}}\left(1,{\mathbf{n}}\right)$
${\mathbf{t}}\left(\mathit{i}+1\right)$ must contain the value ${\tau }_{\mathit{i}}$, for $\mathit{i}=0,1,\dots ,{\mathbf{n}}-1$.
Constraint: ${\mathbf{t}}\left(1\right)>0.0$. Note that if this is not true, then the Toeplitz matrix cannot be positive definite.
2:     $\mathrm{b}\left(:\right)$ – double array
The dimension of the array b must be at least $\mathrm{max}\phantom{\rule{0.125em}{0ex}}\left(1,{\mathbf{n}}\right)$
The right-hand side vector ${b}_{n}$.
3:     $\mathrm{x}\left(:\right)$ – double array
The dimension of the array x must be at least $\mathrm{max}\phantom{\rule{0.125em}{0ex}}\left(1,{\mathbf{n}}\right)$
With ${\mathbf{n}}>1$ the ($n-1$) elements of the solution vector ${x}_{n-1}$ as returned by a previous call to nag_linsys_real_toeplitz_update (f04mf). The element ${\mathbf{x}}\left({\mathbf{n}}\right)$ need not be specified.
4:     $\mathrm{work}\left(:\right)$ – double array
The dimension of the array work must be at least $\mathrm{max}\phantom{\rule{0.125em}{0ex}}\left(1,2×{\mathbf{n}}-1\right)$
With ${\mathbf{n}}>2$ the elements of work should be as returned from a previous call to nag_linsys_real_toeplitz_update (f04mf) with (${\mathbf{n}}-1$) as the argument n.

### Optional Input Parameters

1:     $\mathrm{n}$int64int32nag_int scalar
Default: the dimension of the arrays t, b, x.
The order of the Toeplitz matrix $T$.
Constraint: ${\mathbf{n}}\ge 0$. When ${\mathbf{n}}=0$, then an immediate return is effected.

### Output Parameters

1:     $\mathrm{x}\left(:\right)$ – double array
The dimension of the array x will be $\mathrm{max}\phantom{\rule{0.125em}{0ex}}\left(1,{\mathbf{n}}\right)$
The solution vector ${x}_{n}$.
2:     $\mathrm{p}$ – double scalar
The reflection coefficient ${p}_{n-1}$. (See Further Comments.)
3:     $\mathrm{work}\left(:\right)$ – double array
The dimension of the array work will be $\mathrm{max}\phantom{\rule{0.125em}{0ex}}\left(1,2×{\mathbf{n}}-1\right)$
The first (${\mathbf{n}}-1$) elements of work contain the solution to the Yule–Walker equations
 $Tn-1yn-1=-tn-1,$
where ${t}_{n-1}={\left({\tau }_{1}{\tau }_{2}\dots {\tau }_{n-1}\right)}^{{\mathbf{t}}}$.
4:     $\mathrm{ifail}$int64int32nag_int scalar
${\mathbf{ifail}}={\mathbf{0}}$ unless the function detects an error (see Error Indicators and Warnings).

## Error Indicators and Warnings

Errors or warnings detected by the function:

Cases prefixed with W are classified as warnings and do not generate an error of type NAG:error_n. See nag_issue_warnings.

${\mathbf{ifail}}=-1$
 On entry, ${\mathbf{n}}<0$, or ${\mathbf{t}}\left(0\right)\le 0.0$.
W  ${\mathbf{ifail}}=1$
The Toeplitz matrix ${T}_{n}$ is not positive definite to working accuracy. If, on exit, p is close to unity, then ${T}_{n}$ was probably close to being singular.
${\mathbf{ifail}}=-99$
An unexpected error has been triggered by this routine. Please contact NAG.
${\mathbf{ifail}}=-399$
Your licence key may have expired or may not have been installed correctly.
${\mathbf{ifail}}=-999$
Dynamic memory allocation failed.

## Accuracy

The computed solution of the equations certainly satisfies
 $r=Tnxn-bn,$
where ${‖r‖}_{1}$ is approximately bounded by
 $r1≤cεCTn,$
$c$ being a modest function of $n$, $\epsilon$ being the machine precision and $C\left(T\right)$ being the condition number of $T$ with respect to inversion. This bound is almost certainly pessimistic, but it seems unlikely that the method of Levinson is backward stable, so caution should be exercised when ${T}_{n}$ is ill-conditioned. The following bound on ${T}_{n}^{-1}$ holds:
 $max1∏i=1 n-11-pi2 , 1∏i=1 n-11-pi ≤Tn-11≤∏i=1 n-1 1+pi 1-pi .$
(See Golub and Van Loan (1996).) The norm of ${T}_{n}^{-1}$ may also be estimated using function nag_linsys_real_gen_norm_rcomm (f04yd). For further information on stability issues see Bunch (1985), Bunch (1987), Cybenko (1980) and Golub and Van Loan (1996).

The number of floating-point operations used by this function is approximately $8n$.
If ${y}_{i}$ is the solution of the equations
 $Tiyi=-τ1τ2…τiT,$
then the reflection coefficient ${p}_{i}$ is defined as the $i$th element of ${y}_{i}$.

## Example

This example finds the solution of the equations ${T}_{k}{x}_{k}={b}_{k}$, $k=1,2,3,4$, where
 $T4= 4 3 2 1 3 4 3 2 2 3 4 3 1 2 3 4 and b4= 1 1 1 1 .$
```function f04mf_example

fprintf('f04mf example results\n\n');

t = [4; 3; 2; 1];
b = [1; 1; 1; 1];
x = [0];
work = zeros(9,1);
fprintf('  order  refl. coeff   solution\n');
for k=1:4
[x, p, work, ifail] = f04mf( ...
t(1:k), b(1:k), x, work);
if k > 1
fprintf('%6d%11.4f    ', k, p);
else
fprintf('%6d%11.4s    ', k, ' ');
end
fprintf('%8.4f',transpose(x));
fprintf('\n');
if k < 4
x = [x; 0]; % Extend x by one element
end
end

```
```f04mf example results

order  refl. coeff   solution
1                 0.2500
2    -0.7500      0.1429  0.1429
3     0.1429      0.1667  0.0000  0.1667
4     0.1667      0.2000  0.0000 -0.0000  0.2000
```

PDF version (NAG web site, 64-bit version, 64-bit version)
Chapter Contents
Chapter Introduction
NAG Toolbox

© The Numerical Algorithms Group Ltd, Oxford, UK. 2009–2015