Free cookie consent management tool by TermsFeed Policy Generator

source: trunk/sources/ALGLIB/lq.cs @ 2575

Last change on this file since 2575 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: 12.8 KB
Line 
1/*************************************************************************
2Copyright (c) 2005-2007, Sergey Bochkanov (ALGLIB project).
3
4>>> SOURCE LICENSE >>>
5This program is free software; you can redistribute it and/or modify
6it under the terms of the GNU General Public License as published by
7the Free Software Foundation (www.fsf.org); either version 2 of the
8License, or (at your option) any later version.
9
10This program is distributed in the hope that it will be useful,
11but WITHOUT ANY WARRANTY; without even the implied warranty of
12MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
13GNU General Public License for more details.
14
15A copy of the GNU General Public License is available at
16http://www.fsf.org/licensing/licenses
17
18>>> END OF LICENSE >>>
19*************************************************************************/
20
21using System;
22
23namespace alglib
24{
25    public class lq
26    {
27        /*************************************************************************
28        LQ decomposition of a rectangular matrix of size MxN
29
30        Input parameters:
31            A   -   matrix A whose indexes range within [0..M-1, 0..N-1].
32            M   -   number of rows in matrix A.
33            N   -   number of columns in matrix A.
34
35        Output parameters:
36            A   -   matrices L and Q in compact form (see below)
37            Tau -   array of scalar factors which are used to form
38                    matrix Q. Array whose index ranges within [0..Min(M,N)-1].
39
40        Matrix A is represented as A = LQ, where Q is an orthogonal matrix of size
41        MxM, L - lower triangular (or lower trapezoid) matrix of size M x N.
42
43        The elements of matrix L are located on and below  the  main  diagonal  of
44        matrix A. The elements which are located in Tau array and above  the  main
45        diagonal of matrix A are used to form matrix Q as follows:
46
47        Matrix Q is represented as a product of elementary reflections
48
49        Q = H(k-1)*H(k-2)*...*H(1)*H(0),
50
51        where k = min(m,n), and each H(i) is of the form
52
53        H(i) = 1 - tau * v * (v^T)
54
55        where tau is a scalar stored in Tau[I]; v - real vector, so that v(0:i-1)=0,
56        v(i) = 1, v(i+1:n-1) stored in A(i,i+1:n-1).
57
58          -- ALGLIB --
59             Copyright 2005-2007 by Bochkanov Sergey
60        *************************************************************************/
61        public static void rmatrixlq(ref double[,] a,
62            int m,
63            int n,
64            ref double[] tau)
65        {
66            double[] work = new double[0];
67            double[] t = new double[0];
68            int i = 0;
69            int k = 0;
70            int minmn = 0;
71            int maxmn = 0;
72            double tmp = 0;
73            int i_ = 0;
74            int i1_ = 0;
75
76            minmn = Math.Min(m, n);
77            maxmn = Math.Max(m, n);
78            work = new double[m+1];
79            t = new double[n+1];
80            tau = new double[minmn-1+1];
81            k = Math.Min(m, n);
82            for(i=0; i<=k-1; i++)
83            {
84               
85                //
86                // Generate elementary reflector H(i) to annihilate A(i,i+1:n-1)
87                //
88                i1_ = (i) - (1);
89                for(i_=1; i_<=n-i;i_++)
90                {
91                    t[i_] = a[i,i_+i1_];
92                }
93                reflections.generatereflection(ref t, n-i, ref tmp);
94                tau[i] = tmp;
95                i1_ = (1) - (i);
96                for(i_=i; i_<=n-1;i_++)
97                {
98                    a[i,i_] = t[i_+i1_];
99                }
100                t[1] = 1;
101                if( i<n )
102                {
103                   
104                    //
105                    // Apply H(i) to A(i+1:m,i:n) from the right
106                    //
107                    reflections.applyreflectionfromtheright(ref a, tau[i], ref t, i+1, m-1, i, n-1, ref work);
108                }
109            }
110        }
111
112
113        /*************************************************************************
114        Partial unpacking of matrix Q from the LQ decomposition of a matrix A
115
116        Input parameters:
117            A       -   matrices L and Q in compact form.
118                        Output of RMatrixLQ subroutine.
119            M       -   number of rows in given matrix A. M>=0.
120            N       -   number of columns in given matrix A. N>=0.
121            Tau     -   scalar factors which are used to form Q.
122                        Output of the RMatrixLQ subroutine.
123            QRows   -   required number of rows in matrix Q. N>=QRows>=0.
124
125        Output parameters:
126            Q       -   first QRows rows of matrix Q. Array whose indexes range
127                        within [0..QRows-1, 0..N-1]. If QRows=0, the array remains
128                        unchanged.
129
130          -- ALGLIB --
131             Copyright 2005 by Bochkanov Sergey
132        *************************************************************************/
133        public static void rmatrixlqunpackq(ref double[,] a,
134            int m,
135            int n,
136            ref double[] tau,
137            int qrows,
138            ref double[,] q)
139        {
140            int i = 0;
141            int j = 0;
142            int k = 0;
143            int minmn = 0;
144            double[] v = new double[0];
145            double[] work = new double[0];
146            int i_ = 0;
147            int i1_ = 0;
148
149            System.Diagnostics.Debug.Assert(qrows<=n, "RMatrixLQUnpackQ: QRows>N!");
150            if( m<=0 | n<=0 | qrows<=0 )
151            {
152                return;
153            }
154           
155            //
156            // init
157            //
158            minmn = Math.Min(m, n);
159            k = Math.Min(minmn, qrows);
160            q = new double[qrows-1+1, n-1+1];
161            v = new double[n+1];
162            work = new double[qrows+1];
163            for(i=0; i<=qrows-1; i++)
164            {
165                for(j=0; j<=n-1; j++)
166                {
167                    if( i==j )
168                    {
169                        q[i,j] = 1;
170                    }
171                    else
172                    {
173                        q[i,j] = 0;
174                    }
175                }
176            }
177           
178            //
179            // unpack Q
180            //
181            for(i=k-1; i>=0; i--)
182            {
183               
184                //
185                // Apply H(i)
186                //
187                i1_ = (i) - (1);
188                for(i_=1; i_<=n-i;i_++)
189                {
190                    v[i_] = a[i,i_+i1_];
191                }
192                v[1] = 1;
193                reflections.applyreflectionfromtheright(ref q, tau[i], ref v, 0, qrows-1, i, n-1, ref work);
194            }
195        }
196
197
198        /*************************************************************************
199        Unpacking of matrix L from the LQ decomposition of a matrix A
200
201        Input parameters:
202            A       -   matrices Q and L in compact form.
203                        Output of RMatrixLQ subroutine.
204            M       -   number of rows in given matrix A. M>=0.
205            N       -   number of columns in given matrix A. N>=0.
206
207        Output parameters:
208            L       -   matrix L, array[0..M-1, 0..N-1].
209
210          -- ALGLIB --
211             Copyright 2005 by Bochkanov Sergey
212        *************************************************************************/
213        public static void rmatrixlqunpackl(ref double[,] a,
214            int m,
215            int n,
216            ref double[,] l)
217        {
218            int i = 0;
219            int k = 0;
220            int i_ = 0;
221
222            if( m<=0 | n<=0 )
223            {
224                return;
225            }
226            l = new double[m-1+1, n-1+1];
227            for(i=0; i<=n-1; i++)
228            {
229                l[0,i] = 0;
230            }
231            for(i=1; i<=m-1; i++)
232            {
233                for(i_=0; i_<=n-1;i_++)
234                {
235                    l[i,i_] = l[0,i_];
236                }
237            }
238            for(i=0; i<=m-1; i++)
239            {
240                k = Math.Min(i, n-1);
241                for(i_=0; i_<=k;i_++)
242                {
243                    l[i,i_] = a[i,i_];
244                }
245            }
246        }
247
248
249        public static void lqdecomposition(ref double[,] a,
250            int m,
251            int n,
252            ref double[] tau)
253        {
254            double[] work = new double[0];
255            double[] t = new double[0];
256            int i = 0;
257            int k = 0;
258            int nmip1 = 0;
259            int minmn = 0;
260            int maxmn = 0;
261            double tmp = 0;
262            int i_ = 0;
263            int i1_ = 0;
264
265            minmn = Math.Min(m, n);
266            maxmn = Math.Max(m, n);
267            work = new double[m+1];
268            t = new double[n+1];
269            tau = new double[minmn+1];
270           
271            //
272            // Test the input arguments
273            //
274            k = Math.Min(m, n);
275            for(i=1; i<=k; i++)
276            {
277               
278                //
279                // Generate elementary reflector H(i) to annihilate A(i,i+1:n)
280                //
281                nmip1 = n-i+1;
282                i1_ = (i) - (1);
283                for(i_=1; i_<=nmip1;i_++)
284                {
285                    t[i_] = a[i,i_+i1_];
286                }
287                reflections.generatereflection(ref t, nmip1, ref tmp);
288                tau[i] = tmp;
289                i1_ = (1) - (i);
290                for(i_=i; i_<=n;i_++)
291                {
292                    a[i,i_] = t[i_+i1_];
293                }
294                t[1] = 1;
295                if( i<n )
296                {
297                   
298                    //
299                    // Apply H(i) to A(i+1:m,i:n) from the right
300                    //
301                    reflections.applyreflectionfromtheright(ref a, tau[i], ref t, i+1, m, i, n, ref work);
302                }
303            }
304        }
305
306
307        public static void unpackqfromlq(ref double[,] a,
308            int m,
309            int n,
310            ref double[] tau,
311            int qrows,
312            ref double[,] q)
313        {
314            int i = 0;
315            int j = 0;
316            int k = 0;
317            int minmn = 0;
318            double[] v = new double[0];
319            double[] work = new double[0];
320            int vm = 0;
321            int i_ = 0;
322            int i1_ = 0;
323
324            System.Diagnostics.Debug.Assert(qrows<=n, "UnpackQFromLQ: QRows>N!");
325            if( m==0 | n==0 | qrows==0 )
326            {
327                return;
328            }
329           
330            //
331            // init
332            //
333            minmn = Math.Min(m, n);
334            k = Math.Min(minmn, qrows);
335            q = new double[qrows+1, n+1];
336            v = new double[n+1];
337            work = new double[qrows+1];
338            for(i=1; i<=qrows; i++)
339            {
340                for(j=1; j<=n; j++)
341                {
342                    if( i==j )
343                    {
344                        q[i,j] = 1;
345                    }
346                    else
347                    {
348                        q[i,j] = 0;
349                    }
350                }
351            }
352           
353            //
354            // unpack Q
355            //
356            for(i=k; i>=1; i--)
357            {
358               
359                //
360                // Apply H(i)
361                //
362                vm = n-i+1;
363                i1_ = (i) - (1);
364                for(i_=1; i_<=vm;i_++)
365                {
366                    v[i_] = a[i,i_+i1_];
367                }
368                v[1] = 1;
369                reflections.applyreflectionfromtheright(ref q, tau[i], ref v, 1, qrows, i, n, ref work);
370            }
371        }
372
373
374        public static void lqdecompositionunpacked(double[,] a,
375            int m,
376            int n,
377            ref double[,] l,
378            ref double[,] q)
379        {
380            int i = 0;
381            int j = 0;
382            double[] tau = new double[0];
383
384            a = (double[,])a.Clone();
385
386            if( n<=0 )
387            {
388                return;
389            }
390            q = new double[n+1, n+1];
391            l = new double[m+1, n+1];
392           
393            //
394            // LQDecomposition
395            //
396            lqdecomposition(ref a, m, n, ref tau);
397           
398            //
399            // L
400            //
401            for(i=1; i<=m; i++)
402            {
403                for(j=1; j<=n; j++)
404                {
405                    if( j>i )
406                    {
407                        l[i,j] = 0;
408                    }
409                    else
410                    {
411                        l[i,j] = a[i,j];
412                    }
413                }
414            }
415           
416            //
417            // Q
418            //
419            unpackqfromlq(ref a, m, n, ref tau, n, ref q);
420        }
421    }
422}
Note: See TracBrowser for help on using the repository browser.