**[M-5] schurd()** -- Schur decomposition

__Syntax__

*void* **schurd(***X***,** *T***,** *Q***)**

*void* **_schurd(***X***,** *Q***)**

*void* ** schurdgroupby(***X***,** *f***,** *T***,** *Q***,** *w***,** *m***)**

*void* **_schurdgroupby(***X***,** *f***,** *Q***,** *w***,** *m***)**

where inputs are

*X*: *numeric matrix*
*f*: *pointer scalar* (points to a function used to group eigenvalues)

and outputs are

*T*: *numeric matrix* (Schur-form matrix)
*Q*: *numeric matrix* (orthogonal or unitary)
*w*: *numeric vector* of eigenvalues
*m*: *real scalar* (the number of eigenvalues satisfy the
grouping condition)

__Description__
**schurd(***X***,** *T***,** *Q***)** computes the Schur decomposition of a square, numeric
matrix, *X*, returning the Schur-form matrix, *T*, and the matrix of Schur
vectors, *Q*. *Q* is orthogonal if *X* is real and unitary if *X* is complex.

**_schurd(***X***,** *Q***)** does the same thing as **schurd()**, except that it returns *T*
in *X*.

**schurdgroupby(***X***,** *f***,** *T***,** *Q***,** *w***,** *m***)** computes the Schur decomposition and the
eigenvalues of a square, numeric matrix, *X*, and groups the results
according to whether a condition on each eigenvalue is satisfied.
**schurdgroupby()** returns the Schur-form matrix in *T*, the matrix of Schur
vectors in *Q*, the eigenvalues in *w*, and the number of eigenvalues for
which the condition is true in *m*. *f* is a pointer of the function that
implements the condition on each eigenvalue, as discussed below.

**_schurdgroupby(***X***,** *f***,** *Q***,** *w***,** *m***)** does the same thing as **schurdgroupby()**
except that it returns *T* in *X*.

**_schurd_la()** and **_schurdgroupby_la()** are the interfaces into the LAPACK
routines used to implement the above functions; see **[M-1] LAPACK**. Their
direct use is not recommended.

__Remarks__

Remarks are presented under the following headings:

Schur decomposition
Grouping the results

__Schur decomposition__

Many algorithms begin by obtaining the Schur decomposition of a square
matrix.

The Schur decomposition of matrix **X** can be written as

**Q**' * **X** * **Q** = **T**
where **T** is in Schur form, **Q**, the matrix of Schur vectors, is orthogonal
if **X** is real or unitary if **X** is complex.

A real, square matrix is in Schur form if it is block upper triangular
with 1 *x* 1 and 2 *x* 2 diagonal blocks. Each 2 *x* 2 diagonal block has equal
diagonal elements and opposite sign off-diagonal elements. A complex,
square matrix is in Schur form if it is upper triangular. The
eigenvalues of **X** are obtained from the Schur form by a few quick
computations.

In the example below, we define **X**, obtain the Schur decomposition, and
list **T**.

**: X=(.31,.69,.13,.56\.31,.5,.72,.42 \.68,.37,.71,.8 \.09,.16,.83,.9)**
**: schurd(X, T=., Q=.)**
**: T**
1 2 3 4
+-------------------------------------------------------------+
1 | 2.10742167 .1266712792 .0549744934 .3329112999 |
2 | 0 -.0766307549 .3470959084 .1042286546 |
3 | 0 -.4453774705 -.0766307549 .3000409803 |
4 | 0 0 0 .4658398402 |
+-------------------------------------------------------------+

__Grouping the results__

In many applications, there is a stable solution if the modulus of an
eigenvalue is less than one and an explosive solution if the modulus is
greater than or equal to one. One frequently handles these cases
differently and would group the Schur decomposition results into a block
corresponding to stable solutions and a block corresponding to explosive
solutions.

In the following example, we use **schurdgroupby()** to put the stable
solutions first. One of the arguments to **schurdgroupby()** is a pointer to
a function that accepts a complex scalar argument, an eigenvalue, and
returns 1 to select the eigenvalue and 0 otherwise. Here **isstable()**
returns 1 if the eigenvalue is less than 1:

**: real scalar isstable(scalar p)**
**> {**
**> return((abs(p)<1))**
**> }**
Using this function to group the results, we see that the Schur-form
matrix has been reordered.

**: schurdgroupby(X, &isstable(), T=., Q=., w=., m=.)**

**: T**
1 2 3 4
+-------------------------------------------------------------+
1 | -.0766307549 .445046622 .3029641608 -.0341867415 |
2 | -.3473539401 -.0766307549 -.1036266286 .0799058566 |
3 | 0 0 .4658398402 -.3475944606 |
4 | 0 0 0 2.10742167 |
+-------------------------------------------------------------+
Listing the moduli of the eigenvalues reveals that they are grouped into
stable and explosive groups.

**: abs(w)**
1 2 3 4
+---------------------------------------------------------+
1 | .4005757984 .4005757984 .4658398402 2.10742167 |
+---------------------------------------------------------+

**m** contains the number of stable solutions

**: m**
3

__Conformability__

**schurd(***X***,** *T***,** *Q* **)**:
*input:*
*X*: *n x n*
*output:*
*T*: *n x n*
*Q*: *n x n*
**_schurd(***X***,** *Q***)**:
*input:*
*X*: *n x n*
*output:*
*X*: *n x n*
*Q*: *n x n*
**schurdgroupby(***X***,** *f***,** *T***,** *Q***,** *w***,** *m***)**:
*input:*
*X*: *n x n*
*f*: 1 *x* 1
*output:*
*T*: *n x n*
*Q*: *n x n*
*w*: 1 *x n*
*m*: 1 *x* 1

**_schurdgroupby(***X***,** *f***,** *Q***,** *w***,** *m***)**:
*input:*
*X*: *n x n*
*f*: 1 *x* 1
*output:*
*X*: *n x n*
*Q*: *n x n*
*w*: 1 *x n*
*m*: 1 *x* 1

__Diagnostics__

**_schurd()** and **_schurdgroupby()** abort with error if *X* is a view.

**schurd()**, **_schurd()**, **schurdgroupby()**, and** _schurdgroupby()** return missing
results if *X* contains missing values.

**schurdgroupby()** groups the results via a matrix transform. If the
problem is very ill conditioned, applying this matrix transform can
cause changes in the eigenvalues. In extreme cases, the grouped
eigenvalues may no longer satisfy the condition used to perform the
grouping.

__Source code__

schurd.mata, _schurd.mata, schurdgroupby.mata, _schurdgroupby.mata