Free cookie consent management tool by TermsFeed Policy Generator

source: branches/Persistence Test/ALGLIB/lq.cs @ 3772

Last change on this file since 3772 was 2430, checked in by gkronber, 15 years ago

Moved ALGLIB code into a separate plugin. #783

File size: 13.5 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        /*************************************************************************
250        Obsolete 1-based subroutine
251        See RMatrixLQ for 0-based replacement.
252        *************************************************************************/
253        public static void lqdecomposition(ref double[,] a,
254            int m,
255            int n,
256            ref double[] tau)
257        {
258            double[] work = new double[0];
259            double[] t = new double[0];
260            int i = 0;
261            int k = 0;
262            int nmip1 = 0;
263            int minmn = 0;
264            int maxmn = 0;
265            double tmp = 0;
266            int i_ = 0;
267            int i1_ = 0;
268
269            minmn = Math.Min(m, n);
270            maxmn = Math.Max(m, n);
271            work = new double[m+1];
272            t = new double[n+1];
273            tau = new double[minmn+1];
274           
275            //
276            // Test the input arguments
277            //
278            k = Math.Min(m, n);
279            for(i=1; i<=k; i++)
280            {
281               
282                //
283                // Generate elementary reflector H(i) to annihilate A(i,i+1:n)
284                //
285                nmip1 = n-i+1;
286                i1_ = (i) - (1);
287                for(i_=1; i_<=nmip1;i_++)
288                {
289                    t[i_] = a[i,i_+i1_];
290                }
291                reflections.generatereflection(ref t, nmip1, ref tmp);
292                tau[i] = tmp;
293                i1_ = (1) - (i);
294                for(i_=i; i_<=n;i_++)
295                {
296                    a[i,i_] = t[i_+i1_];
297                }
298                t[1] = 1;
299                if( i<n )
300                {
301                   
302                    //
303                    // Apply H(i) to A(i+1:m,i:n) from the right
304                    //
305                    reflections.applyreflectionfromtheright(ref a, tau[i], ref t, i+1, m, i, n, ref work);
306                }
307            }
308        }
309
310
311        /*************************************************************************
312        Obsolete 1-based subroutine
313        See RMatrixLQUnpackQ for 0-based replacement.
314        *************************************************************************/
315        public static void unpackqfromlq(ref double[,] a,
316            int m,
317            int n,
318            ref double[] tau,
319            int qrows,
320            ref double[,] q)
321        {
322            int i = 0;
323            int j = 0;
324            int k = 0;
325            int minmn = 0;
326            double[] v = new double[0];
327            double[] work = new double[0];
328            int vm = 0;
329            int i_ = 0;
330            int i1_ = 0;
331
332            System.Diagnostics.Debug.Assert(qrows<=n, "UnpackQFromLQ: QRows>N!");
333            if( m==0 | n==0 | qrows==0 )
334            {
335                return;
336            }
337           
338            //
339            // init
340            //
341            minmn = Math.Min(m, n);
342            k = Math.Min(minmn, qrows);
343            q = new double[qrows+1, n+1];
344            v = new double[n+1];
345            work = new double[qrows+1];
346            for(i=1; i<=qrows; i++)
347            {
348                for(j=1; j<=n; j++)
349                {
350                    if( i==j )
351                    {
352                        q[i,j] = 1;
353                    }
354                    else
355                    {
356                        q[i,j] = 0;
357                    }
358                }
359            }
360           
361            //
362            // unpack Q
363            //
364            for(i=k; i>=1; i--)
365            {
366               
367                //
368                // Apply H(i)
369                //
370                vm = n-i+1;
371                i1_ = (i) - (1);
372                for(i_=1; i_<=vm;i_++)
373                {
374                    v[i_] = a[i,i_+i1_];
375                }
376                v[1] = 1;
377                reflections.applyreflectionfromtheright(ref q, tau[i], ref v, 1, qrows, i, n, ref work);
378            }
379        }
380
381
382        /*************************************************************************
383        Obsolete 1-based subroutine
384        *************************************************************************/
385        public static void lqdecompositionunpacked(double[,] a,
386            int m,
387            int n,
388            ref double[,] l,
389            ref double[,] q)
390        {
391            int i = 0;
392            int j = 0;
393            double[] tau = new double[0];
394
395            a = (double[,])a.Clone();
396
397            if( n<=0 )
398            {
399                return;
400            }
401            q = new double[n+1, n+1];
402            l = new double[m+1, n+1];
403           
404            //
405            // LQDecomposition
406            //
407            lqdecomposition(ref a, m, n, ref tau);
408           
409            //
410            // L
411            //
412            for(i=1; i<=m; i++)
413            {
414                for(j=1; j<=n; j++)
415                {
416                    if( j>i )
417                    {
418                        l[i,j] = 0;
419                    }
420                    else
421                    {
422                        l[i,j] = a[i,j];
423                    }
424                }
425            }
426           
427            //
428            // Q
429            //
430            unpackqfromlq(ref a, m, n, ref tau, n, ref q);
431        }
432    }
433}
Note: See TracBrowser for help on using the repository browser.