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/clu.cs

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

extracted external libraries and adapted dependent plugins (ticket #837)
File size: 10.4 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 clu
32    {
33        /*************************************************************************
34        LU decomposition of a complex general matrix of size MxN
35
36        The subroutine calculates the LU decomposition of a rectangular general
37        matrix with partial pivoting (with row permutations).
38
39        Input parameters:
40            A   -   matrix A whose indexes range within [0..M-1, 0..N-1].
41            M   -   number of rows in matrix A.
42            N   -   number of columns in matrix A.
43
44        Output parameters:
45            A   -   matrices L and U in compact form (see below).
46                    Array whose indexes range within [0..M-1, 0..N-1].
47            Pivots - permutation matrix in compact form (see below).
48                    Array whose index ranges within [0..Min(M-1,N-1)].
49
50        Matrix A is represented as A = P * L * U, where P is a permutation matrix,
51        matrix L - lower triangular (or lower trapezoid, if M>N) matrix,
52        U - upper triangular (or upper trapezoid, if M<N) matrix.
53
54        Let M be equal to 4 and N be equal to 3:
55
56                           (  1          )    ( U11 U12 U13  )
57        A = P1 * P2 * P3 * ( L21  1      )  * (     U22 U23  )
58                           ( L31 L32  1  )    (         U33  )
59                           ( L41 L42 L43 )
60
61        Matrix L has size MxMin(M,N), matrix U has size Min(M,N)xN, matrix P(i) is
62        a permutation of the identity matrix of size MxM with numbers I and Pivots[I].
63
64        The algorithm returns array Pivots and the following matrix which replaces
65        matrix A and contains matrices L and U in compact form (the example applies
66        to M=4, N=3).
67
68         ( U11 U12 U13 )
69         ( L21 U22 U23 )
70         ( L31 L32 U33 )
71         ( L41 L42 L43 )
72
73        As we can see, the unit diagonal isn't stored.
74
75          -- LAPACK routine (version 3.0) --
76             Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
77             Courant Institute, Argonne National Lab, and Rice University
78             June 30, 1992
79        *************************************************************************/
80        public static void cmatrixlu(ref AP.Complex[,] a,
81            int m,
82            int n,
83            ref int[] pivots)
84        {
85            int i = 0;
86            int j = 0;
87            int jp = 0;
88            AP.Complex[] t1 = new AP.Complex[0];
89            AP.Complex s = 0;
90            int i_ = 0;
91
92            pivots = new int[Math.Min(m-1, n-1)+1];
93            t1 = new AP.Complex[Math.Max(m-1, n-1)+1];
94            System.Diagnostics.Debug.Assert(m>=0 & n>=0, "Error in LUDecomposition: incorrect function arguments");
95           
96            //
97            // Quick return if possible
98            //
99            if( m==0 | n==0 )
100            {
101                return;
102            }
103            for(j=0; j<=Math.Min(m-1, n-1); j++)
104            {
105               
106                //
107                // Find pivot and test for singularity.
108                //
109                jp = j;
110                for(i=j+1; i<=m-1; i++)
111                {
112                    if( (double)(AP.Math.AbsComplex(a[i,j]))>(double)(AP.Math.AbsComplex(a[jp,j])) )
113                    {
114                        jp = i;
115                    }
116                }
117                pivots[j] = jp;
118                if( a[jp,j]!=0 )
119                {
120                   
121                    //
122                    //Apply the interchange to rows
123                    //
124                    if( jp!=j )
125                    {
126                        for(i_=0; i_<=n-1;i_++)
127                        {
128                            t1[i_] = a[j,i_];
129                        }
130                        for(i_=0; i_<=n-1;i_++)
131                        {
132                            a[j,i_] = a[jp,i_];
133                        }
134                        for(i_=0; i_<=n-1;i_++)
135                        {
136                            a[jp,i_] = t1[i_];
137                        }
138                    }
139                   
140                    //
141                    //Compute elements J+1:M of J-th column.
142                    //
143                    if( j<m )
144                    {
145                        jp = j+1;
146                        s = 1/a[j,j];
147                        for(i_=jp; i_<=m-1;i_++)
148                        {
149                            a[i_,j] = s*a[i_,j];
150                        }
151                    }
152                }
153                if( j<Math.Min(m, n)-1 )
154                {
155                   
156                    //
157                    //Update trailing submatrix.
158                    //
159                    jp = j+1;
160                    for(i=j+1; i<=m-1; i++)
161                    {
162                        s = a[i,j];
163                        for(i_=jp; i_<=n-1;i_++)
164                        {
165                            a[i,i_] = a[i,i_] - s*a[j,i_];
166                        }
167                    }
168                }
169            }
170        }
171
172
173        public static void complexludecomposition(ref AP.Complex[,] a,
174            int m,
175            int n,
176            ref int[] pivots)
177        {
178            int i = 0;
179            int j = 0;
180            int jp = 0;
181            AP.Complex[] t1 = new AP.Complex[0];
182            AP.Complex s = 0;
183            int i_ = 0;
184
185            pivots = new int[Math.Min(m, n)+1];
186            t1 = new AP.Complex[Math.Max(m, n)+1];
187            System.Diagnostics.Debug.Assert(m>=0 & n>=0, "Error in ComplexLUDecomposition: incorrect function arguments");
188           
189            //
190            // Quick return if possible
191            //
192            if( m==0 | n==0 )
193            {
194                return;
195            }
196            for(j=1; j<=Math.Min(m, n); j++)
197            {
198               
199                //
200                // Find pivot and test for singularity.
201                //
202                jp = j;
203                for(i=j+1; i<=m; i++)
204                {
205                    if( (double)(AP.Math.AbsComplex(a[i,j]))>(double)(AP.Math.AbsComplex(a[jp,j])) )
206                    {
207                        jp = i;
208                    }
209                }
210                pivots[j] = jp;
211                if( a[jp,j]!=0 )
212                {
213                   
214                    //
215                    //Apply the interchange to rows
216                    //
217                    if( jp!=j )
218                    {
219                        for(i_=1; i_<=n;i_++)
220                        {
221                            t1[i_] = a[j,i_];
222                        }
223                        for(i_=1; i_<=n;i_++)
224                        {
225                            a[j,i_] = a[jp,i_];
226                        }
227                        for(i_=1; i_<=n;i_++)
228                        {
229                            a[jp,i_] = t1[i_];
230                        }
231                    }
232                   
233                    //
234                    //Compute elements J+1:M of J-th column.
235                    //
236                    if( j<m )
237                    {
238                        jp = j+1;
239                        s = 1/a[j,j];
240                        for(i_=jp; i_<=m;i_++)
241                        {
242                            a[i_,j] = s*a[i_,j];
243                        }
244                    }
245                }
246                if( j<Math.Min(m, n) )
247                {
248                   
249                    //
250                    //Update trailing submatrix.
251                    //
252                    jp = j+1;
253                    for(i=j+1; i<=m; i++)
254                    {
255                        s = a[i,j];
256                        for(i_=jp; i_<=n;i_++)
257                        {
258                            a[i,i_] = a[i,i_] - s*a[j,i_];
259                        }
260                    }
261                }
262            }
263        }
264
265
266        public static void complexludecompositionunpacked(AP.Complex[,] a,
267            int m,
268            int n,
269            ref AP.Complex[,] l,
270            ref AP.Complex[,] u,
271            ref int[] pivots)
272        {
273            int i = 0;
274            int j = 0;
275            int minmn = 0;
276
277            a = (AP.Complex[,])a.Clone();
278
279            if( m==0 | n==0 )
280            {
281                return;
282            }
283            minmn = Math.Min(m, n);
284            l = new AP.Complex[m+1, minmn+1];
285            u = new AP.Complex[minmn+1, n+1];
286            complexludecomposition(ref a, m, n, ref pivots);
287            for(i=1; i<=m; i++)
288            {
289                for(j=1; j<=minmn; j++)
290                {
291                    if( j>i )
292                    {
293                        l[i,j] = 0;
294                    }
295                    if( j==i )
296                    {
297                        l[i,j] = 1;
298                    }
299                    if( j<i )
300                    {
301                        l[i,j] = a[i,j];
302                    }
303                }
304            }
305            for(i=1; i<=minmn; i++)
306            {
307                for(j=1; j<=n; j++)
308                {
309                    if( j<i )
310                    {
311                        u[i,j] = 0;
312                    }
313                    if( j>=i )
314                    {
315                        u[i,j] = a[i,j];
316                    }
317                }
318            }
319        }
320    }
321}
Note: See TracBrowser for help on using the repository browser.