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/cinverse.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: 8.6 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 cinverse
32    {
33        /*************************************************************************
34        Inversion of a complex matrix given by its LU decomposition.
35
36        Input parameters:
37            A       -   LU decomposition of the matrix (output of CMatrixLU subroutine).
38            Pivots  -   table of permutations which were made during the LU decomposition
39                        (the output of CMatrixLU subroutine).
40            N       -   size of matrix A.
41
42        Output parameters:
43            A       -   inverse of matrix A.
44                        Array whose indexes range within [0..N-1, 0..N-1].
45
46        Result:
47            True, if the matrix is not singular.
48            False, if the matrix is singular.
49
50          -- LAPACK routine (version 3.0) --
51             Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
52             Courant Institute, Argonne National Lab, and Rice University
53             February 29, 1992
54        *************************************************************************/
55        public static bool cmatrixluinverse(ref AP.Complex[,] a,
56            ref int[] pivots,
57            int n)
58        {
59            bool result = new bool();
60            AP.Complex[] work = new AP.Complex[0];
61            int i = 0;
62            int iws = 0;
63            int j = 0;
64            int jb = 0;
65            int jj = 0;
66            int jp = 0;
67            AP.Complex v = 0;
68            int i_ = 0;
69
70            result = true;
71           
72            //
73            // Quick return if possible
74            //
75            if( n==0 )
76            {
77                return result;
78            }
79            work = new AP.Complex[n-1+1];
80           
81            //
82            // Form inv(U)
83            //
84            if( !ctrinverse.cmatrixtrinverse(ref a, n, true, false) )
85            {
86                result = false;
87                return result;
88            }
89           
90            //
91            // Solve the equation inv(A)*L = inv(U) for inv(A).
92            //
93            for(j=n-1; j>=0; j--)
94            {
95               
96                //
97                // Copy current column of L to WORK and replace with zeros.
98                //
99                for(i=j+1; i<=n-1; i++)
100                {
101                    work[i] = a[i,j];
102                    a[i,j] = 0;
103                }
104               
105                //
106                // Compute current column of inv(A).
107                //
108                if( j<n-1 )
109                {
110                    for(i=0; i<=n-1; i++)
111                    {
112                        v = 0.0;
113                        for(i_=j+1; i_<=n-1;i_++)
114                        {
115                            v += a[i,i_]*work[i_];
116                        }
117                        a[i,j] = a[i,j]-v;
118                    }
119                }
120            }
121           
122            //
123            // Apply column interchanges.
124            //
125            for(j=n-2; j>=0; j--)
126            {
127                jp = pivots[j];
128                if( jp!=j )
129                {
130                    for(i_=0; i_<=n-1;i_++)
131                    {
132                        work[i_] = a[i_,j];
133                    }
134                    for(i_=0; i_<=n-1;i_++)
135                    {
136                        a[i_,j] = a[i_,jp];
137                    }
138                    for(i_=0; i_<=n-1;i_++)
139                    {
140                        a[i_,jp] = work[i_];
141                    }
142                }
143            }
144            return result;
145        }
146
147
148        /*************************************************************************
149        Inversion of a general complex matrix.
150
151        Input parameters:
152            A   -   matrix. Array whose indexes range within [0..N-1, 0..N-1].
153            N   -   size of matrix A.
154
155        Output parameters:
156            A   -   inverse of matrix A.
157                    Array whose indexes range within [0..N-1, 0..N-1].
158
159        Result:
160            True, if the matrix is not singular.
161            False, if the matrix is singular.
162
163          -- ALGLIB --
164             Copyright 2005 by Bochkanov Sergey
165        *************************************************************************/
166        public static bool cmatrixinverse(ref AP.Complex[,] a,
167            int n)
168        {
169            bool result = new bool();
170            int[] pivots = new int[0];
171
172            clu.cmatrixlu(ref a, n, n, ref pivots);
173            result = cmatrixluinverse(ref a, ref pivots, n);
174            return result;
175        }
176
177
178        public static bool complexinverselu(ref AP.Complex[,] a,
179            ref int[] pivots,
180            int n)
181        {
182            bool result = new bool();
183            AP.Complex[] work = new AP.Complex[0];
184            int i = 0;
185            int iws = 0;
186            int j = 0;
187            int jb = 0;
188            int jj = 0;
189            int jp = 0;
190            int jp1 = 0;
191            AP.Complex v = 0;
192            int i_ = 0;
193
194            result = true;
195           
196            //
197            // Quick return if possible
198            //
199            if( n==0 )
200            {
201                return result;
202            }
203            work = new AP.Complex[n+1];
204           
205            //
206            // Form inv(U)
207            //
208            if( !ctrinverse.complexinvtriangular(ref a, n, true, false) )
209            {
210                result = false;
211                return result;
212            }
213           
214            //
215            // Solve the equation inv(A)*L = inv(U) for inv(A).
216            //
217            for(j=n; j>=1; j--)
218            {
219               
220                //
221                // Copy current column of L to WORK and replace with zeros.
222                //
223                for(i=j+1; i<=n; i++)
224                {
225                    work[i] = a[i,j];
226                    a[i,j] = 0;
227                }
228               
229                //
230                // Compute current column of inv(A).
231                //
232                if( j<n )
233                {
234                    jp1 = j+1;
235                    for(i=1; i<=n; i++)
236                    {
237                        v = 0.0;
238                        for(i_=jp1; i_<=n;i_++)
239                        {
240                            v += a[i,i_]*work[i_];
241                        }
242                        a[i,j] = a[i,j]-v;
243                    }
244                }
245            }
246           
247            //
248            // Apply column interchanges.
249            //
250            for(j=n-1; j>=1; j--)
251            {
252                jp = pivots[j];
253                if( jp!=j )
254                {
255                    for(i_=1; i_<=n;i_++)
256                    {
257                        work[i_] = a[i_,j];
258                    }
259                    for(i_=1; i_<=n;i_++)
260                    {
261                        a[i_,j] = a[i_,jp];
262                    }
263                    for(i_=1; i_<=n;i_++)
264                    {
265                        a[i_,jp] = work[i_];
266                    }
267                }
268            }
269            return result;
270        }
271
272
273        public static bool complexinverse(ref AP.Complex[,] a,
274            int n)
275        {
276            bool result = new bool();
277            int[] pivots = new int[0];
278
279            clu.complexludecomposition(ref a, n, n, ref pivots);
280            result = complexinverselu(ref a, ref pivots, n);
281            return result;
282        }
283    }
284}
Note: See TracBrowser for help on using the repository browser.