Free cookie consent management tool by TermsFeed Policy Generator

source: trunk/sources/ALGLIB/lu.cs @ 2643

Last change on this file since 2643 was 2563, checked in by gkronber, 15 years ago

Updated ALGLIB to latest version. #751 (Plugin for for data-modeling with ANN (integrated into CEDMA))

File size: 17.0 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 lu
32    {
33        public const int lunb = 8;
34
35
36        /*************************************************************************
37        LU decomposition of a general matrix of size MxN
38
39        The subroutine calculates the LU decomposition of a rectangular general
40        matrix with partial pivoting (with row permutations).
41
42        Input parameters:
43            A   -   matrix A whose indexes range within [0..M-1, 0..N-1].
44            M   -   number of rows in matrix A.
45            N   -   number of columns in matrix A.
46
47        Output parameters:
48            A   -   matrices L and U in compact form (see below).
49                    Array whose indexes range within [0..M-1, 0..N-1].
50            Pivots - permutation matrix in compact form (see below).
51                    Array whose index ranges within [0..Min(M-1,N-1)].
52
53        Matrix A is represented as A = P * L * U, where P is a permutation matrix,
54        matrix L - lower triangular (or lower trapezoid, if M>N) matrix,
55        U - upper triangular (or upper trapezoid, if M<N) matrix.
56
57        Let M be equal to 4 and N be equal to 3:
58
59                           (  1          )    ( U11 U12 U13  )
60        A = P1 * P2 * P3 * ( L21  1      )  * (     U22 U23  )
61                           ( L31 L32  1  )    (         U33  )
62                           ( L41 L42 L43 )
63
64        Matrix L has size MxMin(M,N), matrix U has size Min(M,N)xN, matrix P(i) is
65        a permutation of the identity matrix of size MxM with numbers I and Pivots[I].
66
67        The algorithm returns array Pivots and the following matrix which replaces
68        matrix A and contains matrices L and U in compact form (the example applies
69        to M=4, N=3).
70
71         ( U11 U12 U13 )
72         ( L21 U22 U23 )
73         ( L31 L32 U33 )
74         ( L41 L42 L43 )
75
76        As we can see, the unit diagonal isn't stored.
77
78          -- LAPACK routine (version 3.0) --
79             Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
80             Courant Institute, Argonne National Lab, and Rice University
81             June 30, 1992
82        *************************************************************************/
83        public static void rmatrixlu(ref double[,] a,
84            int m,
85            int n,
86            ref int[] pivots)
87        {
88            double[,] b = new double[0,0];
89            double[] t = new double[0];
90            int[] bp = new int[0];
91            int minmn = 0;
92            int i = 0;
93            int ip = 0;
94            int j = 0;
95            int j1 = 0;
96            int j2 = 0;
97            int cb = 0;
98            int nb = 0;
99            double v = 0;
100            int i_ = 0;
101            int i1_ = 0;
102
103            System.Diagnostics.Debug.Assert(lunb>=1, "RMatrixLU internal error");
104            nb = lunb;
105           
106            //
107            // Decide what to use - blocked or unblocked code
108            //
109            if( n<=1 | Math.Min(m, n)<=nb | nb==1 )
110            {
111               
112                //
113                // Unblocked code
114                //
115                rmatrixlu2(ref a, m, n, ref pivots);
116            }
117            else
118            {
119               
120                //
121                // Blocked code.
122                // First, prepare temporary matrix and indices
123                //
124                b = new double[m-1+1, nb-1+1];
125                t = new double[n-1+1];
126                pivots = new int[Math.Min(m, n)-1+1];
127                minmn = Math.Min(m, n);
128                j1 = 0;
129                j2 = Math.Min(minmn, nb)-1;
130               
131                //
132                // Main cycle
133                //
134                while( j1<minmn )
135                {
136                    cb = j2-j1+1;
137                   
138                    //
139                    // LU factorization of diagonal and subdiagonal blocks:
140                    // 1. Copy columns J1..J2 of A to B
141                    // 2. LU(B)
142                    // 3. Copy result back to A
143                    // 4. Copy pivots, apply pivots
144                    //
145                    for(i=j1; i<=m-1; i++)
146                    {
147                        i1_ = (j1) - (0);
148                        for(i_=0; i_<=cb-1;i_++)
149                        {
150                            b[i-j1,i_] = a[i,i_+i1_];
151                        }
152                    }
153                    rmatrixlu2(ref b, m-j1, cb, ref bp);
154                    for(i=j1; i<=m-1; i++)
155                    {
156                        i1_ = (0) - (j1);
157                        for(i_=j1; i_<=j2;i_++)
158                        {
159                            a[i,i_] = b[i-j1,i_+i1_];
160                        }
161                    }
162                    for(i=0; i<=cb-1; i++)
163                    {
164                        ip = bp[i];
165                        pivots[j1+i] = j1+ip;
166                        if( bp[i]!=i )
167                        {
168                            if( j1!=0 )
169                            {
170                               
171                                //
172                                // Interchange columns 0:J1-1
173                                //
174                                for(i_=0; i_<=j1-1;i_++)
175                                {
176                                    t[i_] = a[j1+i,i_];
177                                }
178                                for(i_=0; i_<=j1-1;i_++)
179                                {
180                                    a[j1+i,i_] = a[j1+ip,i_];
181                                }
182                                for(i_=0; i_<=j1-1;i_++)
183                                {
184                                    a[j1+ip,i_] = t[i_];
185                                }
186                            }
187                            if( j2<n-1 )
188                            {
189                               
190                                //
191                                // Interchange the rest of the matrix, if needed
192                                //
193                                for(i_=j2+1; i_<=n-1;i_++)
194                                {
195                                    t[i_] = a[j1+i,i_];
196                                }
197                                for(i_=j2+1; i_<=n-1;i_++)
198                                {
199                                    a[j1+i,i_] = a[j1+ip,i_];
200                                }
201                                for(i_=j2+1; i_<=n-1;i_++)
202                                {
203                                    a[j1+ip,i_] = t[i_];
204                                }
205                            }
206                        }
207                    }
208                   
209                    //
210                    // Compute block row of U
211                    //
212                    if( j2<n-1 )
213                    {
214                        for(i=j1+1; i<=j2; i++)
215                        {
216                            for(j=j1; j<=i-1; j++)
217                            {
218                                v = a[i,j];
219                                for(i_=j2+1; i_<=n-1;i_++)
220                                {
221                                    a[i,i_] = a[i,i_] - v*a[j,i_];
222                                }
223                            }
224                        }
225                    }
226                   
227                    //
228                    // Update trailing submatrix
229                    //
230                    if( j2<n-1 )
231                    {
232                        for(i=j2+1; i<=m-1; i++)
233                        {
234                            for(j=j1; j<=j2; j++)
235                            {
236                                v = a[i,j];
237                                for(i_=j2+1; i_<=n-1;i_++)
238                                {
239                                    a[i,i_] = a[i,i_] - v*a[j,i_];
240                                }
241                            }
242                        }
243                    }
244                   
245                    //
246                    // Next step
247                    //
248                    j1 = j2+1;
249                    j2 = Math.Min(minmn, j1+nb)-1;
250                }
251            }
252        }
253
254
255        public static void ludecomposition(ref double[,] a,
256            int m,
257            int n,
258            ref int[] pivots)
259        {
260            int i = 0;
261            int j = 0;
262            int jp = 0;
263            double[] t1 = new double[0];
264            double s = 0;
265            int i_ = 0;
266
267            pivots = new int[Math.Min(m, n)+1];
268            t1 = new double[Math.Max(m, n)+1];
269            System.Diagnostics.Debug.Assert(m>=0 & n>=0, "Error in LUDecomposition: incorrect function arguments");
270           
271            //
272            // Quick return if possible
273            //
274            if( m==0 | n==0 )
275            {
276                return;
277            }
278            for(j=1; j<=Math.Min(m, n); j++)
279            {
280               
281                //
282                // Find pivot and test for singularity.
283                //
284                jp = j;
285                for(i=j+1; i<=m; i++)
286                {
287                    if( (double)(Math.Abs(a[i,j]))>(double)(Math.Abs(a[jp,j])) )
288                    {
289                        jp = i;
290                    }
291                }
292                pivots[j] = jp;
293                if( (double)(a[jp,j])!=(double)(0) )
294                {
295                   
296                    //
297                    //Apply the interchange to rows
298                    //
299                    if( jp!=j )
300                    {
301                        for(i_=1; i_<=n;i_++)
302                        {
303                            t1[i_] = a[j,i_];
304                        }
305                        for(i_=1; i_<=n;i_++)
306                        {
307                            a[j,i_] = a[jp,i_];
308                        }
309                        for(i_=1; i_<=n;i_++)
310                        {
311                            a[jp,i_] = t1[i_];
312                        }
313                    }
314                   
315                    //
316                    //Compute elements J+1:M of J-th column.
317                    //
318                    if( j<m )
319                    {
320                       
321                        //
322                        // CALL DSCAL( M-J, ONE / A( J, J ), A( J+1, J ), 1 )
323                        //
324                        jp = j+1;
325                        s = 1/a[j,j];
326                        for(i_=jp; i_<=m;i_++)
327                        {
328                            a[i_,j] = s*a[i_,j];
329                        }
330                    }
331                }
332                if( j<Math.Min(m, n) )
333                {
334                   
335                    //
336                    //Update trailing submatrix.
337                    //CALL DGER( M-J, N-J, -ONE, A( J+1, J ), 1, A( J, J+1 ), LDA,A( J+1, J+1 ), LDA )
338                    //
339                    jp = j+1;
340                    for(i=j+1; i<=m; i++)
341                    {
342                        s = a[i,j];
343                        for(i_=jp; i_<=n;i_++)
344                        {
345                            a[i,i_] = a[i,i_] - s*a[j,i_];
346                        }
347                    }
348                }
349            }
350        }
351
352
353        public static void ludecompositionunpacked(double[,] a,
354            int m,
355            int n,
356            ref double[,] l,
357            ref double[,] u,
358            ref int[] pivots)
359        {
360            int i = 0;
361            int j = 0;
362            int minmn = 0;
363
364            a = (double[,])a.Clone();
365
366            if( m==0 | n==0 )
367            {
368                return;
369            }
370            minmn = Math.Min(m, n);
371            l = new double[m+1, minmn+1];
372            u = new double[minmn+1, n+1];
373            ludecomposition(ref a, m, n, ref pivots);
374            for(i=1; i<=m; i++)
375            {
376                for(j=1; j<=minmn; j++)
377                {
378                    if( j>i )
379                    {
380                        l[i,j] = 0;
381                    }
382                    if( j==i )
383                    {
384                        l[i,j] = 1;
385                    }
386                    if( j<i )
387                    {
388                        l[i,j] = a[i,j];
389                    }
390                }
391            }
392            for(i=1; i<=minmn; i++)
393            {
394                for(j=1; j<=n; j++)
395                {
396                    if( j<i )
397                    {
398                        u[i,j] = 0;
399                    }
400                    if( j>=i )
401                    {
402                        u[i,j] = a[i,j];
403                    }
404                }
405            }
406        }
407
408
409        /*************************************************************************
410        Level 2 BLAS version of RMatrixLU
411
412          -- LAPACK routine (version 3.0) --
413             Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
414             Courant Institute, Argonne National Lab, and Rice University
415             June 30, 1992
416        *************************************************************************/
417        private static void rmatrixlu2(ref double[,] a,
418            int m,
419            int n,
420            ref int[] pivots)
421        {
422            int i = 0;
423            int j = 0;
424            int jp = 0;
425            double[] t1 = new double[0];
426            double s = 0;
427            int i_ = 0;
428
429            pivots = new int[Math.Min(m-1, n-1)+1];
430            t1 = new double[Math.Max(m-1, n-1)+1];
431            System.Diagnostics.Debug.Assert(m>=0 & n>=0, "Error in LUDecomposition: incorrect function arguments");
432           
433            //
434            // Quick return if possible
435            //
436            if( m==0 | n==0 )
437            {
438                return;
439            }
440            for(j=0; j<=Math.Min(m-1, n-1); j++)
441            {
442               
443                //
444                // Find pivot and test for singularity.
445                //
446                jp = j;
447                for(i=j+1; i<=m-1; i++)
448                {
449                    if( (double)(Math.Abs(a[i,j]))>(double)(Math.Abs(a[jp,j])) )
450                    {
451                        jp = i;
452                    }
453                }
454                pivots[j] = jp;
455                if( (double)(a[jp,j])!=(double)(0) )
456                {
457                   
458                    //
459                    //Apply the interchange to rows
460                    //
461                    if( jp!=j )
462                    {
463                        for(i_=0; i_<=n-1;i_++)
464                        {
465                            t1[i_] = a[j,i_];
466                        }
467                        for(i_=0; i_<=n-1;i_++)
468                        {
469                            a[j,i_] = a[jp,i_];
470                        }
471                        for(i_=0; i_<=n-1;i_++)
472                        {
473                            a[jp,i_] = t1[i_];
474                        }
475                    }
476                   
477                    //
478                    //Compute elements J+1:M of J-th column.
479                    //
480                    if( j<m )
481                    {
482                        jp = j+1;
483                        s = 1/a[j,j];
484                        for(i_=jp; i_<=m-1;i_++)
485                        {
486                            a[i_,j] = s*a[i_,j];
487                        }
488                    }
489                }
490                if( j<Math.Min(m, n)-1 )
491                {
492                   
493                    //
494                    //Update trailing submatrix.
495                    //
496                    jp = j+1;
497                    for(i=j+1; i<=m-1; i++)
498                    {
499                        s = a[i,j];
500                        for(i_=jp; i_<=n-1;i_++)
501                        {
502                            a[i,i_] = a[i,i_] - s*a[j,i_];
503                        }
504                    }
505                }
506            }
507        }
508    }
509}
Note: See TracBrowser for help on using the repository browser.