# NAG FL Interfacef08rnf (zuncsd)

## 1Purpose

f08rnf computes the CS decomposition of a complex $m$ by $m$ unitary matrix $X$, partitioned into a $2$ by $2$ array of submatrices.

## 2Specification

Fortran Interface
 Subroutine f08rnf ( m, p, q, x11, x12, x21, x22, u1, ldu1, u2, ldu2, v1t, v2t, work, info)
 Integer, Intent (In) :: m, p, q, ldx11, ldx12, ldx21, ldx22, ldu1, ldu2, ldv1t, ldv2t, lwork, lrwork Integer, Intent (Out) :: iwork(m-min(p,m-p,q,m-q)), info Real (Kind=nag_wp), Intent (Out) :: theta(min(p,m-p,q,m-q)), rwork(max(1,lrwork)) Complex (Kind=nag_wp), Intent (Inout) :: x11(ldx11,*), x12(ldx12,*), x21(ldx21,*), x22(ldx22,*), u1(ldu1,*), u2(ldu2,*), v1t(ldv1t,*), v2t(ldv2t,*) Complex (Kind=nag_wp), Intent (Out) :: work(max(1,lwork)) Character (1), Intent (In) :: jobu1, jobu2, jobv1t, jobv2t, trans, signs
#include <nag.h>
 void f08rnf_ (const char *jobu1, const char *jobu2, const char *jobv1t, const char *jobv2t, const char *trans, const char *signs, const Integer *m, const Integer *p, const Integer *q, Complex x11[], const Integer *ldx11, Complex x12[], const Integer *ldx12, Complex x21[], const Integer *ldx21, Complex x22[], const Integer *ldx22, double theta[], Complex u1[], const Integer *ldu1, Complex u2[], const Integer *ldu2, Complex v1t[], const Integer *ldv1t, Complex v2t[], const Integer *ldv2t, Complex work[], const Integer *lwork, double rwork[], const Integer *lrwork, Integer iwork[], Integer *info, const Charlen length_jobu1, const Charlen length_jobu2, const Charlen length_jobv1t, const Charlen length_jobv2t, const Charlen length_trans, const Charlen length_signs)
The routine may be called by the names f08rnf, nagf_lapackeig_zuncsd or its LAPACK name zuncsd.

## 3Description

The $m$ by $m$ unitary matrix $X$ is partitioned as
 $X= X11 X12 X21 X22$
where ${X}_{11}$ is a $p$ by $q$ submatrix and the dimensions of the other submatrices ${X}_{12}$, ${X}_{21}$ and ${X}_{22}$ are such that $X$ remains $m$ by $m$.
The CS decomposition of $X$ is $X=U{\Sigma }_{p}{V}^{\mathrm{T}}$ where $U$, $V$ and ${\Sigma }_{p}$ are $m$ by $m$ matrices, such that
 $U= U1 0 0 U2$
is a unitary matrix containing the $p$ by $p$ unitary matrix ${U}_{1}$ and the $\left(m-p\right)$ by $\left(m-p\right)$ unitary matrix ${U}_{2}$;
 $V= V1 0 0 V2$
is a unitary matrix containing the $q$ by $q$ unitary matrix ${V}_{1}$ and the $\left(m-q\right)$ by $\left(m-q\right)$ unitary matrix ${V}_{2}$; and
 $Σp= I11 0 0 0 C 0 0 -S 0 0 0 -I12 0 0 I22 0 0 S C 0 0 I21 0 0$
contains the $r$ by $r$ non-negative diagonal submatrices $C$ and $S$ satisfying ${C}^{2}+{S}^{2}=I$, where $r=\mathrm{min}\phantom{\rule{0.125em}{0ex}}\left(p,m-p,q,m-q\right)$ and the top left partition is $p$ by $q$.
The identity matrix ${I}_{11}$ is of order $\mathrm{min}\phantom{\rule{0.125em}{0ex}}\left(p,q\right)-r$ and vanishes if $\mathrm{min}\phantom{\rule{0.125em}{0ex}}\left(p,q\right)=r$.
The identity matrix ${I}_{12}$ is of order $\mathrm{min}\phantom{\rule{0.125em}{0ex}}\left(p,m-q\right)-r$ and vanishes if $\mathrm{min}\phantom{\rule{0.125em}{0ex}}\left(p,m-q\right)=r$.
The identity matrix ${I}_{21}$ is of order $\mathrm{min}\phantom{\rule{0.125em}{0ex}}\left(m-p,q\right)-r$ and vanishes if $\mathrm{min}\phantom{\rule{0.125em}{0ex}}\left(m-p,q\right)=r$.
The identity matrix ${I}_{22}$ is of order $\mathrm{min}\phantom{\rule{0.125em}{0ex}}\left(m-p,m-q\right)-r$ and vanishes if $\mathrm{min}\phantom{\rule{0.125em}{0ex}}\left(m-p,m-q\right)=r$.
In each of the four cases $r=p,q,m-p,m-q$ at least two of the identity matrices vanish.
The indicated zeros represent augmentations by additional rows or columns (but not both) to the square diagonal matrices formed by ${I}_{ij}$ and $C$ or $S$.
${\Sigma }_{p}$ does not need to be stored in full; it is sufficient to return only the values ${\theta }_{i}$ for $i=1,2,\dots ,r$ where ${C}_{ii}=\mathrm{cos}\left({\theta }_{i}\right)$ and ${S}_{ii}=\mathrm{sin}\left({\theta }_{i}\right)$.
The algorithm used to perform the complete $CS$ decomposition is described fully in Sutton (2009) including discussions of the stability and accuracy of the algorithm.
Anderson E, Bai Z, Bischof C, Blackford S, Demmel J, Dongarra J J, Du Croz J J, Greenbaum A, Hammarling S, McKenney A and Sorensen D (1999) LAPACK Users' Guide (3rd Edition) SIAM, Philadelphia https://www.netlib.org/lapack/lug
Golub G H and Van Loan C F (2012) Matrix Computations (4th Edition) Johns Hopkins University Press, Baltimore
Sutton B D (2009) Computing the complete $CS$ decomposition Numerical Algorithms (Volume 50) 1017–1398 Springer US 33–65 https://dx.doi.org/10.1007/s11075-008-9215-6

## 5Arguments

1: $\mathbf{jobu1}$Character(1) Input
On entry:
• if ${\mathbf{jobu1}}=\text{'Y'}$, ${U}_{1}$ is computed;
• otherwise, ${U}_{1}$ is not computed.
2: $\mathbf{jobu2}$Character(1) Input
On entry:
• if ${\mathbf{jobu2}}=\text{'Y'}$, ${U}_{2}$ is computed;
• otherwise, ${U}_{2}$ is not computed.
3: $\mathbf{jobv1t}$Character(1) Input
On entry:
• if ${\mathbf{jobv1t}}=\text{'Y'}$, ${V}_{1}^{\mathrm{T}}$ is computed;
• otherwise, ${V}_{1}^{\mathrm{T}}$ is not computed.
4: $\mathbf{jobv2t}$Character(1) Input
On entry:
• if ${\mathbf{jobv2t}}=\text{'Y'}$, ${V}_{2}^{\mathrm{T}}$ is computed;
• otherwise, ${V}_{2}^{\mathrm{T}}$ is not computed.
5: $\mathbf{trans}$Character(1) Input
On entry:
• if ${\mathbf{trans}}=\text{'T'}$, $X$, ${U}_{1}$, ${U}_{2}$, ${V}_{1}^{\mathrm{T}}$ and ${V}_{2}^{\mathrm{T}}$ are stored in row-major order;
• otherwise, $X$, ${U}_{1}$, ${U}_{2}$, ${V}_{1}^{\mathrm{T}}$ and ${V}_{2}^{\mathrm{T}}$ are stored in column-major order.
6: $\mathbf{signs}$Character(1) Input
On entry:
• if ${\mathbf{signs}}=\text{'O'}$, the lower-left block is made nonpositive (the other convention);
• otherwise, the upper-right block is made nonpositive (the default convention).
7: $\mathbf{m}$Integer Input
On entry: $m$, the number of rows and columns in the unitary matrix $X$.
Constraint: ${\mathbf{m}}\ge 0$.
8: $\mathbf{p}$Integer Input
On entry: $p$, the number of rows in ${X}_{11}$ and ${X}_{12}$.
Constraint: $0\le {\mathbf{p}}\le {\mathbf{m}}$.
9: $\mathbf{q}$Integer Input
On entry: $q$, the number of columns in ${X}_{11}$ and ${X}_{21}$.
Constraint: $0\le {\mathbf{q}}\le {\mathbf{m}}$.
10: $\mathbf{x11}\left({\mathbf{ldx11}},*\right)$Complex (Kind=nag_wp) array Input/Output
Note: the second dimension of the array x11 must be at least $\mathrm{max}\phantom{\rule{0.125em}{0ex}}\left(1,{\mathbf{p}}\right)$ if ${\mathbf{trans}}=\text{'T'}$, and at least ${\mathbf{q}}$ otherwise.
On entry: the upper left partition of the unitary matrix $X$ whose CSD is desired.
On exit: contains details of the unitary matrix used in a simultaneous bidiagonalization process.
11: $\mathbf{ldx11}$Integer Input
On entry: the first dimension of the array x11 as declared in the (sub)program from which f08rnf is called.
Constraints:
• if ${\mathbf{trans}}=\text{'T'}$, ${\mathbf{ldx11}}\ge \mathrm{max}\phantom{\rule{0.125em}{0ex}}\left(1,{\mathbf{q}}\right)$;
• otherwise ${\mathbf{ldx11}}\ge \mathrm{max}\phantom{\rule{0.125em}{0ex}}\left(1,{\mathbf{p}}\right)$.
12: $\mathbf{x12}\left({\mathbf{ldx12}},*\right)$Complex (Kind=nag_wp) array Input/Output
Note: the second dimension of the array x12 must be at least $\mathrm{max}\phantom{\rule{0.125em}{0ex}}\left(1,{\mathbf{p}}\right)$ if ${\mathbf{trans}}=\text{'T'}$, and at least ${\mathbf{m}}-{\mathbf{q}}$ otherwise.
On entry: the upper right partition of the unitary matrix $X$ whose CSD is desired.
On exit: contains details of the unitary matrix used in a simultaneous bidiagonalization process.
13: $\mathbf{ldx12}$Integer Input
On entry: the first dimension of the array x12 as declared in the (sub)program from which f08rnf is called.
Constraints:
• if ${\mathbf{trans}}=\text{'T'}$, ${\mathbf{ldx12}}\ge \mathrm{max}\phantom{\rule{0.125em}{0ex}}\left(1,{\mathbf{m}}-{\mathbf{q}}\right)$;
• otherwise ${\mathbf{ldx12}}\ge \mathrm{max}\phantom{\rule{0.125em}{0ex}}\left(1,{\mathbf{p}}\right)$.
14: $\mathbf{x21}\left({\mathbf{ldx21}},*\right)$Complex (Kind=nag_wp) array Input/Output
Note: the second dimension of the array x21 must be at least $\mathrm{max}\phantom{\rule{0.125em}{0ex}}\left(1,{\mathbf{m}}-{\mathbf{p}}\right)$ if ${\mathbf{trans}}=\text{'T'}$, and at least ${\mathbf{q}}$ otherwise.
On entry: the lower left partition of the unitary matrix $X$ whose CSD is desired.
On exit: contains details of the unitary matrix used in a simultaneous bidiagonalization process.
15: $\mathbf{ldx21}$Integer Input
On entry: the first dimension of the array x21 as declared in the (sub)program from which f08rnf is called.
Constraints:
• if ${\mathbf{trans}}=\text{'T'}$, ${\mathbf{ldx21}}\ge \mathrm{max}\phantom{\rule{0.125em}{0ex}}\left(1,{\mathbf{q}}\right)$;
• otherwise ${\mathbf{ldx21}}\ge \mathrm{max}\phantom{\rule{0.125em}{0ex}}\left(1,{\mathbf{m}}-{\mathbf{p}}\right)$.
16: $\mathbf{x22}\left({\mathbf{ldx22}},*\right)$Complex (Kind=nag_wp) array Input/Output
Note: the second dimension of the array x22 must be at least $\mathrm{max}\phantom{\rule{0.125em}{0ex}}\left(1,{\mathbf{m}}-{\mathbf{p}}\right)$ if ${\mathbf{trans}}=\text{'T'}$, and at least ${\mathbf{m}}-{\mathbf{q}}$ otherwise.
On entry: the lower right partition of the unitary matrix $X$ CSD is desired.
On exit: contains details of the unitary matrix used in a simultaneous bidiagonalization process.
17: $\mathbf{ldx22}$Integer Input
On entry: the first dimension of the array x22 as declared in the (sub)program from which f08rnf is called.
Constraints:
• if ${\mathbf{trans}}=\text{'T'}$, ${\mathbf{ldx22}}\ge \mathrm{max}\phantom{\rule{0.125em}{0ex}}\left(1,{\mathbf{m}}-{\mathbf{q}}\right)$;
• otherwise ${\mathbf{ldx22}}\ge \mathrm{max}\phantom{\rule{0.125em}{0ex}}\left(1,{\mathbf{m}}-{\mathbf{p}}\right)$.
18: $\mathbf{theta}\left(\mathrm{min}\phantom{\rule{0.125em}{0ex}}\left({\mathbf{p}},{\mathbf{m}}-{\mathbf{p}},{\mathbf{q}},{\mathbf{m}}-{\mathbf{q}}\right)\right)$Real (Kind=nag_wp) array Output
On exit: the values ${\theta }_{i}$ for $i=1,2,\dots ,r$ where $r=\mathrm{min}\phantom{\rule{0.125em}{0ex}}\left(p,m-p,q,m-q\right)$. The diagonal submatrices $C$ and $S$ of ${\Sigma }_{p}$ are constructed from these values as
• $C=\mathrm{diag}\left(\mathrm{cos}\left({\mathbf{theta}}\left(1\right)\right),\dots ,\mathrm{cos}\left({\mathbf{theta}}\left(r\right)\right)\right)$ and
• $S=\mathrm{diag}\left(\mathrm{sin}\left({\mathbf{theta}}\left(1\right)\right),\dots ,\mathrm{sin}\left({\mathbf{theta}}\left(r\right)\right)\right)$.
19: $\mathbf{u1}\left({\mathbf{ldu1}},*\right)$Complex (Kind=nag_wp) array Output
Note: the second dimension of the array u1 must be at least $\mathrm{max}\phantom{\rule{0.125em}{0ex}}\left(1,{\mathbf{p}}\right)$ if ${\mathbf{jobu1}}=\text{'Y'}$, and at least $1$ otherwise.
On exit: if ${\mathbf{jobu1}}=\text{'Y'}$, u1 contains the $p$ by $p$ unitary matrix ${U}_{1}$.
20: $\mathbf{ldu1}$Integer Input
On entry: the first dimension of the array u1 as declared in the (sub)program from which f08rnf is called.
Constraint: if ${\mathbf{jobu1}}=\text{'Y'}$, ${\mathbf{ldu1}}\ge \mathrm{max}\phantom{\rule{0.125em}{0ex}}\left(1,{\mathbf{p}}\right)$.
21: $\mathbf{u2}\left({\mathbf{ldu2}},*\right)$Complex (Kind=nag_wp) array Output
Note: the second dimension of the array u2 must be at least $\mathrm{max}\phantom{\rule{0.125em}{0ex}}\left(1,{\mathbf{m}}-{\mathbf{p}}\right)$ if ${\mathbf{jobu2}}=\text{'Y'}$, and at least $1$ otherwise.
On exit: if ${\mathbf{jobu2}}=\text{'Y'}$, u2 contains the $m-p$ by $m-p$ unitary matrix ${U}_{2}$.
22: $\mathbf{ldu2}$Integer Input
On entry: the first dimension of the array u2 as declared in the (sub)program from which f08rnf is called.
Constraint: if ${\mathbf{jobu2}}=\text{'Y'}$, ${\mathbf{ldu2}}\ge \mathrm{max}\phantom{\rule{0.125em}{0ex}}\left(1,{\mathbf{m}}-{\mathbf{p}}\right)$.
23: $\mathbf{v1t}\left({\mathbf{ldv1t}},*\right)$Complex (Kind=nag_wp) array Output
Note: the second dimension of the array v1t must be at least $\mathrm{max}\phantom{\rule{0.125em}{0ex}}\left(1,{\mathbf{q}}\right)$ if ${\mathbf{jobv1t}}=\text{'Y'}$, and at least $1$ otherwise.
On exit: if ${\mathbf{jobv1t}}=\text{'Y'}$, v1t contains the $q$ by $q$ unitary matrix ${{V}_{1}}^{\mathrm{H}}$.
24: $\mathbf{ldv1t}$Integer Input
On entry: the first dimension of the array v1t as declared in the (sub)program from which f08rnf is called.
Constraint: if ${\mathbf{jobv1t}}=\text{'Y'}$, ${\mathbf{ldv1t}}\ge \mathrm{max}\phantom{\rule{0.125em}{0ex}}\left(1,{\mathbf{q}}\right)$.
25: $\mathbf{v2t}\left({\mathbf{ldv2t}},*\right)$Complex (Kind=nag_wp) array Output
Note: the second dimension of the array v2t must be at least $\mathrm{max}\phantom{\rule{0.125em}{0ex}}\left(1,{\mathbf{m}}-{\mathbf{q}}\right)$ if ${\mathbf{jobv2t}}=\text{'Y'}$, and at least $1$ otherwise.
On exit: if ${\mathbf{jobv2t}}=\text{'Y'}$, v2t contains the $m-q$ by $m-q$ unitary matrix ${{V}_{2}}^{\mathrm{H}}$.
26: $\mathbf{ldv2t}$Integer Input
On entry: the first dimension of the array v2t as declared in the (sub)program from which f08rnf is called.
Constraint: if ${\mathbf{jobv2t}}=\text{'Y'}$, ${\mathbf{ldv2t}}\ge \mathrm{max}\phantom{\rule{0.125em}{0ex}}\left(1,{\mathbf{m}}-{\mathbf{q}}\right)$.
27: $\mathbf{work}\left(\mathrm{max}\phantom{\rule{0.125em}{0ex}}\left(1,{\mathbf{lwork}}\right)\right)$Complex (Kind=nag_wp) array Workspace
On exit: if ${\mathbf{info}}={\mathbf{0}}$, ${\mathbf{work}}\left(1\right)$ returns the optimal lwork.
If ${\mathbf{info}}>{\mathbf{0}}$ on exit, ${\mathbf{work}}\left(2:\mathit{r}\right)$ contains the values $\mathrm{PHI}\left(1\right),\dots \mathrm{PHI}\left(\mathit{r}-1\right)$ that, together with ${\mathrm{xref}}^{arg}\left(1\right),\dots {\mathrm{xref}}^{arg}\left(\mathit{r}\right)$, define the matrix in intermediate bidiagonal-block form remaining after nonconvergence. info specifies the number of nonzero PHI's.
28: $\mathbf{lwork}$Integer Input
On entry: the dimension of the array work as declared in the (sub)program from which f08rnf is called.
If ${\mathbf{lwork}}=-1$, a workspace query is assumed; the routine only calculates the optimal size of the work array, returns this value as the first entry of the work array, and no error message related to lwork is issued.
The minimum workspace required is $\mathrm{max}\phantom{\rule{0.125em}{0ex}}\left(1,{\mathbf{p}}\right)+\mathrm{max}\phantom{\rule{0.125em}{0ex}}\left(1,{\mathbf{m}}-{\mathbf{p}}\right)+\mathrm{max}\phantom{\rule{0.125em}{0ex}}\left(1,{\mathbf{q}}\right)+\phantom{\rule{0ex}{0ex}}\mathrm{max}\phantom{\rule{0.125em}{0ex}}\left(1,{\mathbf{m}}-{\mathbf{q}}\right)+\mathrm{max}\phantom{\rule{0.125em}{0ex}}\left(1,{\mathbf{p}},{\mathbf{m}}-{\mathbf{p}},{\mathbf{q}},{\mathbf{m}}-{\mathbf{q}}\right)+1$; the optimal amount of workspace depends on internal block sizes and the relative problem dimensions.
Constraint:
${\mathbf{lwork}}=-1$ or ${\mathbf{lwork}}\ge \mathrm{max}\phantom{\rule{0.125em}{0ex}}\left(1,{\mathbf{p}}\right)+\mathrm{max}\phantom{\rule{0.125em}{0ex}}\left(1,{\mathbf{m}}-{\mathbf{p}}\right)+\mathrm{max}\phantom{\rule{0.125em}{0ex}}\left(1,{\mathbf{q}}\right)+\phantom{\rule{0ex}{0ex}}\mathrm{max}\phantom{\rule{0.125em}{0ex}}\left(1,{\mathbf{m}}-{\mathbf{q}}\right)+\mathrm{max}\phantom{\rule{0.125em}{0ex}}\left(1,{\mathbf{p}},{\mathbf{m}}-{\mathbf{p}},{\mathbf{q}},{\mathbf{m}}-{\mathbf{q}}\right)+1$.
29: $\mathbf{rwork}\left(\mathrm{max}\phantom{\rule{0.125em}{0ex}}\left(1,{\mathbf{lrwork}}\right)\right)$Real (Kind=nag_wp) array Workspace
30: $\mathbf{lrwork}$Integer Input
On entry: the dimension of the array rwork as declared in the (sub)program from which f08rnf is called.
If ${\mathbf{lrwork}}=-1$, a workspace query is assumed; the routine only calculates the optimal size of the work array, returns this value as the first entry of the work array, and no error message related to lrwork is issued. Otherwise the required workspace is $5×\mathrm{max}\phantom{\rule{0.125em}{0ex}}\left(1,{\mathbf{q}}-1\right)+4×\mathrm{max}\phantom{\rule{0.125em}{0ex}}\left(1,{\mathbf{q}}\right)+\mathrm{max}\phantom{\rule{0.125em}{0ex}}\left(1,8×{\mathbf{q}}\right)+1$ which equates to $11$ for ${\mathbf{q}}=0$, $18$ for ${\mathbf{q}}=1$ and $17×{\mathbf{q}}-4$ when ${\mathbf{q}}>1$.
Constraint:
${\mathbf{lrwork}}=-1$ or ${\mathbf{lrwork}}\ge 5×\mathrm{max}\phantom{\rule{0.125em}{0ex}}\left(1,{\mathbf{q}}-1\right)+4×\mathrm{max}\phantom{\rule{0.125em}{0ex}}\left(1,{\mathbf{q}}\right)+\mathrm{max}\phantom{\rule{0.125em}{0ex}}\left(1,8×{\mathbf{q}}\right)+1$.
31: $\mathbf{iwork}\left({\mathbf{m}}-\mathrm{min}\phantom{\rule{0.125em}{0ex}}\left({\mathbf{p}},{\mathbf{m}}-{\mathbf{p}},{\mathbf{q}},{\mathbf{m}}-{\mathbf{q}}\right)\right)$Integer array Workspace
32: $\mathbf{info}$Integer Output
On exit: ${\mathbf{info}}=0$ unless the routine detects an error (see Section 6).

## 6Error Indicators and Warnings

${\mathbf{info}}<0$
If ${\mathbf{info}}=-i$, argument $i$ had an illegal value. An explanatory message is output, and execution of the program is terminated.
${\mathbf{info}}>0$
The Jacobi-type procedure failed to converge during an internal reduction to bidiagonal-block form. The process requires convergence to $\mathrm{min}\phantom{\rule{0.125em}{0ex}}\left({\mathbf{p}},{\mathbf{m}}-{\mathbf{p}},{\mathbf{q}},{\mathbf{m}}-{\mathbf{q}}\right)$ values, the value of info gives the number of converged values.

## 7Accuracy

The computed $CS$ decomposition is nearly the exact $CS$ decomposition for the nearby matrix $\left(X+E\right)$, where
 $E2 = Oε ,$
and $\epsilon$ is the machine precision.

## 8Parallelism and Performance

f08rnf is threaded by NAG for parallel execution in multithreaded implementations of the NAG Library.
f08rnf makes calls to BLAS and/or LAPACK routines, which may be threaded within the vendor library used by this implementation. Consult the documentation for the vendor library for further information.
Please consult the X06 Chapter Introduction for information on how to control and interrogate the OpenMP environment used within this routine. Please also consult the Users' Note for your implementation for any additional implementation-specific information.

The total number of floating-point operations required to perform the full $CS$ decomposition is approximately $2{m}^{3}$.
The real analogue of this routine is f08raf.

## 10Example

This example finds the full CS decomposition of a unitary $6$ by $6$ matrix $X$ (see Section 10.2) partitioned so that the top left block is $2$ by $4$.
The decomposition is performed both on submatrices of the unitary matrix $X$ and on separated partition matrices. Code is also provided to perform a recombining check if required.

### 10.1Program Text

Program Text (f08rnfe.f90)

### 10.2Program Data

Program Data (f08rnfe.d)

### 10.3Program Results

Program Results (f08rnfe.r)