Free cookie consent management tool by TermsFeed Policy Generator

source: trunk/sources/HeuristicLab.ExtLibs/HeuristicLab.ALGLIB/2.3.0/ALGLIB-2.3.0/cinverse.cs @ 3021

Last change on this file since 3021 was 2806, checked in by gkronber, 15 years ago

Added plugin for new version of ALGLIB. #875 (Update ALGLIB sources)

File size: 5.7 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            trfac.cmatrixlu(ref a, n, n, ref pivots);
173            result = cmatrixluinverse(ref a, ref pivots, n);
174            return result;
175        }
176    }
177}
Note: See TracBrowser for help on using the repository browser.