[2806] | 1 | /*************************************************************************
|
---|
| 2 | Copyright (c) 1992-2007 The University of Tennessee. All rights reserved.
|
---|
| 3 |
|
---|
| 4 | Contributors:
|
---|
| 5 | * Sergey Bochkanov (ALGLIB project). Translation from FORTRAN to
|
---|
| 6 | pseudocode.
|
---|
| 7 |
|
---|
| 8 | See subroutines comments for additional copyrights.
|
---|
| 9 |
|
---|
| 10 | >>> SOURCE LICENSE >>>
|
---|
| 11 | This program is free software; you can redistribute it and/or modify
|
---|
| 12 | it under the terms of the GNU General Public License as published by
|
---|
| 13 | the Free Software Foundation (www.fsf.org); either version 2 of the
|
---|
| 14 | License, or (at your option) any later version.
|
---|
| 15 |
|
---|
| 16 | This program is distributed in the hope that it will be useful,
|
---|
| 17 | but WITHOUT ANY WARRANTY; without even the implied warranty of
|
---|
| 18 | MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
|
---|
| 19 | GNU General Public License for more details.
|
---|
| 20 |
|
---|
| 21 | A copy of the GNU General Public License is available at
|
---|
| 22 | http://www.fsf.org/licensing/licenses
|
---|
| 23 |
|
---|
| 24 | >>> END OF LICENSE >>>
|
---|
| 25 | *************************************************************************/
|
---|
| 26 |
|
---|
| 27 | using System;
|
---|
| 28 |
|
---|
| 29 | namespace alglib
|
---|
| 30 | {
|
---|
| 31 | public class schur
|
---|
| 32 | {
|
---|
| 33 | /*************************************************************************
|
---|
| 34 | Subroutine performing the Schur decomposition of a general matrix by using
|
---|
| 35 | the QR algorithm with multiple shifts.
|
---|
| 36 |
|
---|
| 37 | The source matrix A is represented as S'*A*S = T, where S is an orthogonal
|
---|
| 38 | matrix (Schur vectors), T - upper quasi-triangular matrix (with blocks of
|
---|
| 39 | sizes 1x1 and 2x2 on the main diagonal).
|
---|
| 40 |
|
---|
| 41 | Input parameters:
|
---|
| 42 | A - matrix to be decomposed.
|
---|
| 43 | Array whose indexes range within [0..N-1, 0..N-1].
|
---|
| 44 | N - size of A, N>=0.
|
---|
| 45 |
|
---|
| 46 |
|
---|
| 47 | Output parameters:
|
---|
| 48 | A - contains matrix T.
|
---|
| 49 | Array whose indexes range within [0..N-1, 0..N-1].
|
---|
| 50 | S - contains Schur vectors.
|
---|
| 51 | Array whose indexes range within [0..N-1, 0..N-1].
|
---|
| 52 |
|
---|
| 53 | Note 1:
|
---|
| 54 | The block structure of matrix T can be easily recognized: since all
|
---|
| 55 | the elements below the blocks are zeros, the elements a[i+1,i] which
|
---|
| 56 | are equal to 0 show the block border.
|
---|
| 57 |
|
---|
| 58 | Note 2:
|
---|
| 59 | The algorithm performance depends on the value of the internal parameter
|
---|
| 60 | NS of the InternalSchurDecomposition subroutine which defines the number
|
---|
| 61 | of shifts in the QR algorithm (similarly to the block width in block-matrix
|
---|
| 62 | algorithms in linear algebra). If you require maximum performance on
|
---|
| 63 | your machine, it is recommended to adjust this parameter manually.
|
---|
| 64 |
|
---|
| 65 | Result:
|
---|
| 66 | True,
|
---|
| 67 | if the algorithm has converged and parameters A and S contain the result.
|
---|
| 68 | False,
|
---|
| 69 | if the algorithm has not converged.
|
---|
| 70 |
|
---|
| 71 | Algorithm implemented on the basis of the DHSEQR subroutine (LAPACK 3.0 library).
|
---|
| 72 | *************************************************************************/
|
---|
| 73 | public static bool rmatrixschur(ref double[,] a,
|
---|
| 74 | int n,
|
---|
| 75 | ref double[,] s)
|
---|
| 76 | {
|
---|
| 77 | bool result = new bool();
|
---|
| 78 | double[] tau = new double[0];
|
---|
| 79 | double[] wi = new double[0];
|
---|
| 80 | double[] wr = new double[0];
|
---|
| 81 | double[,] a1 = new double[0,0];
|
---|
| 82 | double[,] s1 = new double[0,0];
|
---|
| 83 | int info = 0;
|
---|
| 84 | int i = 0;
|
---|
| 85 | int j = 0;
|
---|
| 86 |
|
---|
| 87 |
|
---|
| 88 | //
|
---|
| 89 | // Upper Hessenberg form of the 0-based matrix
|
---|
| 90 | //
|
---|
| 91 | hessenberg.rmatrixhessenberg(ref a, n, ref tau);
|
---|
| 92 | hessenberg.rmatrixhessenbergunpackq(ref a, n, ref tau, ref s);
|
---|
| 93 |
|
---|
| 94 | //
|
---|
| 95 | // Convert from 0-based arrays to 1-based,
|
---|
| 96 | // then call InternalSchurDecomposition
|
---|
| 97 | // Awkward, of course, but Schur decompisiton subroutine
|
---|
| 98 | // is too complex to fix it.
|
---|
| 99 | //
|
---|
| 100 | //
|
---|
| 101 | a1 = new double[n+1, n+1];
|
---|
| 102 | s1 = new double[n+1, n+1];
|
---|
| 103 | for(i=1; i<=n; i++)
|
---|
| 104 | {
|
---|
| 105 | for(j=1; j<=n; j++)
|
---|
| 106 | {
|
---|
| 107 | a1[i,j] = a[i-1,j-1];
|
---|
| 108 | s1[i,j] = s[i-1,j-1];
|
---|
| 109 | }
|
---|
| 110 | }
|
---|
| 111 | hsschur.internalschurdecomposition(ref a1, n, 1, 1, ref wr, ref wi, ref s1, ref info);
|
---|
| 112 | result = info==0;
|
---|
| 113 |
|
---|
| 114 | //
|
---|
| 115 | // convert from 1-based arrays to -based
|
---|
| 116 | //
|
---|
| 117 | for(i=1; i<=n; i++)
|
---|
| 118 | {
|
---|
| 119 | for(j=1; j<=n; j++)
|
---|
| 120 | {
|
---|
| 121 | a[i-1,j-1] = a1[i,j];
|
---|
| 122 | s[i-1,j-1] = s1[i,j];
|
---|
| 123 | }
|
---|
| 124 | }
|
---|
| 125 | return result;
|
---|
| 126 | }
|
---|
| 127 |
|
---|
| 128 |
|
---|
| 129 | public static bool schurdecomposition(ref double[,] a,
|
---|
| 130 | int n,
|
---|
| 131 | ref double[,] s)
|
---|
| 132 | {
|
---|
| 133 | bool result = new bool();
|
---|
| 134 | double[] tau = new double[0];
|
---|
| 135 | double[] wi = new double[0];
|
---|
| 136 | double[] wr = new double[0];
|
---|
| 137 | int info = 0;
|
---|
| 138 |
|
---|
| 139 | hessenberg.toupperhessenberg(ref a, n, ref tau);
|
---|
| 140 | hessenberg.unpackqfromupperhessenberg(ref a, n, ref tau, ref s);
|
---|
| 141 | hsschur.internalschurdecomposition(ref a, n, 1, 1, ref wr, ref wi, ref s, ref info);
|
---|
| 142 | result = info==0;
|
---|
| 143 | return result;
|
---|
| 144 | }
|
---|
| 145 | }
|
---|
| 146 | }
|
---|