Free cookie consent management tool by TermsFeed Policy Generator

source: branches/3.2/sources/HeuristicLab.ExtLibs/HeuristicLab.ALGLIB/2.1.2.2591/ALGLIB-2.1.2.2591/schur.cs @ 4539

Last change on this file since 4539 was 2645, checked in by mkommend, 15 years ago

extracted external libraries and adapted dependent plugins (ticket #837)

File size: 5.3 KB
Line 
1/*************************************************************************
2Copyright (c) 1992-2007 The University of Tennessee.  All rights reserved.
3
4Contributors:
5    * Sergey Bochkanov (ALGLIB project). Translation from FORTRAN to
6      pseudocode.
7
8See subroutines comments for additional copyrights.
9
10>>> SOURCE LICENSE >>>
11This program is free software; you can redistribute it and/or modify
12it under the terms of the GNU General Public License as published by
13the Free Software Foundation (www.fsf.org); either version 2 of the
14License, or (at your option) any later version.
15
16This program is distributed in the hope that it will be useful,
17but WITHOUT ANY WARRANTY; without even the implied warranty of
18MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
19GNU General Public License for more details.
20
21A copy of the GNU General Public License is available at
22http://www.fsf.org/licensing/licenses
23
24>>> END OF LICENSE >>>
25*************************************************************************/
26
27using System;
28
29namespace 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}
Note: See TracBrowser for help on using the repository browser.