Free cookie consent management tool by TermsFeed Policy Generator

source: trunk/sources/ALGLIB/qr.cs @ 2504

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

Moved ALGLIB code into a separate plugin. #783

File size: 14.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 qr
32    {
33        /*************************************************************************
34        QR decomposition of a rectangular matrix of size MxN
35
36        Input parameters:
37            A   -   matrix A whose indexes range within [0..M-1, 0..N-1].
38            M   -   number of rows in matrix A.
39            N   -   number of columns in matrix A.
40
41        Output parameters:
42            A   -   matrices Q and R in compact form (see below).
43            Tau -   array of scalar factors which are used to form
44                    matrix Q. Array whose index ranges within [0.. Min(M-1,N-1)].
45
46        Matrix A is represented as A = QR, where Q is an orthogonal matrix of size
47        MxM, R - upper triangular (or upper trapezoid) matrix of size M x N.
48
49        The elements of matrix R are located on and above the main diagonal of
50        matrix A. The elements which are located in Tau array and below the main
51        diagonal of matrix A are used to form matrix Q as follows:
52
53        Matrix Q is represented as a product of elementary reflections
54
55        Q = H(0)*H(2)*...*H(k-1),
56
57        where k = min(m,n), and each H(i) is in the form
58
59        H(i) = 1 - tau * v * (v^T)
60
61        where tau is a scalar stored in Tau[I]; v - real vector,
62        so that v(0:i-1) = 0, v(i) = 1, v(i+1:m-1) stored in A(i+1:m-1,i).
63
64          -- LAPACK routine (version 3.0) --
65             Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
66             Courant Institute, Argonne National Lab, and Rice University
67             February 29, 1992.
68             Translation from FORTRAN to pseudocode (AlgoPascal)
69             by Sergey Bochkanov, ALGLIB project, 2005-2007.
70        *************************************************************************/
71        public static void rmatrixqr(ref double[,] a,
72            int m,
73            int n,
74            ref double[] tau)
75        {
76            double[] work = new double[0];
77            double[] t = new double[0];
78            int i = 0;
79            int k = 0;
80            int minmn = 0;
81            double tmp = 0;
82            int i_ = 0;
83            int i1_ = 0;
84
85            if( m<=0 | n<=0 )
86            {
87                return;
88            }
89            minmn = Math.Min(m, n);
90            work = new double[n-1+1];
91            t = new double[m+1];
92            tau = new double[minmn-1+1];
93           
94            //
95            // Test the input arguments
96            //
97            k = minmn;
98            for(i=0; i<=k-1; i++)
99            {
100               
101                //
102                // Generate elementary reflector H(i) to annihilate A(i+1:m,i)
103                //
104                i1_ = (i) - (1);
105                for(i_=1; i_<=m-i;i_++)
106                {
107                    t[i_] = a[i_+i1_,i];
108                }
109                reflections.generatereflection(ref t, m-i, ref tmp);
110                tau[i] = tmp;
111                i1_ = (1) - (i);
112                for(i_=i; i_<=m-1;i_++)
113                {
114                    a[i_,i] = t[i_+i1_];
115                }
116                t[1] = 1;
117                if( i<n )
118                {
119                   
120                    //
121                    // Apply H(i) to A(i:m-1,i+1:n-1) from the left
122                    //
123                    reflections.applyreflectionfromtheleft(ref a, tau[i], ref t, i, m-1, i+1, n-1, ref work);
124                }
125            }
126        }
127
128
129        /*************************************************************************
130        Partial unpacking of matrix Q from the QR decomposition of a matrix A
131
132        Input parameters:
133            A       -   matrices Q and R in compact form.
134                        Output of RMatrixQR subroutine.
135            M       -   number of rows in given matrix A. M>=0.
136            N       -   number of columns in given matrix A. N>=0.
137            Tau     -   scalar factors which are used to form Q.
138                        Output of the RMatrixQR subroutine.
139            QColumns -  required number of columns of matrix Q. M>=QColumns>=0.
140
141        Output parameters:
142            Q       -   first QColumns columns of matrix Q.
143                        Array whose indexes range within [0..M-1, 0..QColumns-1].
144                        If QColumns=0, the array remains unchanged.
145
146          -- ALGLIB --
147             Copyright 2005 by Bochkanov Sergey
148        *************************************************************************/
149        public static void rmatrixqrunpackq(ref double[,] a,
150            int m,
151            int n,
152            ref double[] tau,
153            int qcolumns,
154            ref double[,] q)
155        {
156            int i = 0;
157            int j = 0;
158            int k = 0;
159            int minmn = 0;
160            double[] v = new double[0];
161            double[] work = new double[0];
162            int i_ = 0;
163            int i1_ = 0;
164
165            System.Diagnostics.Debug.Assert(qcolumns<=m, "UnpackQFromQR: QColumns>M!");
166            if( m<=0 | n<=0 | qcolumns<=0 )
167            {
168                return;
169            }
170           
171            //
172            // init
173            //
174            minmn = Math.Min(m, n);
175            k = Math.Min(minmn, qcolumns);
176            q = new double[m-1+1, qcolumns-1+1];
177            v = new double[m+1];
178            work = new double[qcolumns-1+1];
179            for(i=0; i<=m-1; i++)
180            {
181                for(j=0; j<=qcolumns-1; j++)
182                {
183                    if( i==j )
184                    {
185                        q[i,j] = 1;
186                    }
187                    else
188                    {
189                        q[i,j] = 0;
190                    }
191                }
192            }
193           
194            //
195            // unpack Q
196            //
197            for(i=k-1; i>=0; i--)
198            {
199               
200                //
201                // Apply H(i)
202                //
203                i1_ = (i) - (1);
204                for(i_=1; i_<=m-i;i_++)
205                {
206                    v[i_] = a[i_+i1_,i];
207                }
208                v[1] = 1;
209                reflections.applyreflectionfromtheleft(ref q, tau[i], ref v, i, m-1, 0, qcolumns-1, ref work);
210            }
211        }
212
213
214        /*************************************************************************
215        Unpacking of matrix R from the QR decomposition of a matrix A
216
217        Input parameters:
218            A       -   matrices Q and R in compact form.
219                        Output of RMatrixQR subroutine.
220            M       -   number of rows in given matrix A. M>=0.
221            N       -   number of columns in given matrix A. N>=0.
222
223        Output parameters:
224            R       -   matrix R, array[0..M-1, 0..N-1].
225
226          -- ALGLIB --
227             Copyright 2005 by Bochkanov Sergey
228        *************************************************************************/
229        public static void rmatrixqrunpackr(ref double[,] a,
230            int m,
231            int n,
232            ref double[,] r)
233        {
234            int i = 0;
235            int k = 0;
236            int i_ = 0;
237
238            if( m<=0 | n<=0 )
239            {
240                return;
241            }
242            k = Math.Min(m, n);
243            r = new double[m-1+1, n-1+1];
244            for(i=0; i<=n-1; i++)
245            {
246                r[0,i] = 0;
247            }
248            for(i=1; i<=m-1; i++)
249            {
250                for(i_=0; i_<=n-1;i_++)
251                {
252                    r[i,i_] = r[0,i_];
253                }
254            }
255            for(i=0; i<=k-1; i++)
256            {
257                for(i_=i; i_<=n-1;i_++)
258                {
259                    r[i,i_] = a[i,i_];
260                }
261            }
262        }
263
264
265        /*************************************************************************
266        Obsolete 1-based subroutine. See RMatrixQR for 0-based replacement.
267        *************************************************************************/
268        public static void qrdecomposition(ref double[,] a,
269            int m,
270            int n,
271            ref double[] tau)
272        {
273            double[] work = new double[0];
274            double[] t = new double[0];
275            int i = 0;
276            int k = 0;
277            int mmip1 = 0;
278            int minmn = 0;
279            double tmp = 0;
280            int i_ = 0;
281            int i1_ = 0;
282
283            minmn = Math.Min(m, n);
284            work = new double[n+1];
285            t = new double[m+1];
286            tau = new double[minmn+1];
287           
288            //
289            // Test the input arguments
290            //
291            k = Math.Min(m, n);
292            for(i=1; i<=k; i++)
293            {
294               
295                //
296                // Generate elementary reflector H(i) to annihilate A(i+1:m,i)
297                //
298                mmip1 = m-i+1;
299                i1_ = (i) - (1);
300                for(i_=1; i_<=mmip1;i_++)
301                {
302                    t[i_] = a[i_+i1_,i];
303                }
304                reflections.generatereflection(ref t, mmip1, ref tmp);
305                tau[i] = tmp;
306                i1_ = (1) - (i);
307                for(i_=i; i_<=m;i_++)
308                {
309                    a[i_,i] = t[i_+i1_];
310                }
311                t[1] = 1;
312                if( i<n )
313                {
314                   
315                    //
316                    // Apply H(i) to A(i:m,i+1:n) from the left
317                    //
318                    reflections.applyreflectionfromtheleft(ref a, tau[i], ref t, i, m, i+1, n, ref work);
319                }
320            }
321        }
322
323
324        /*************************************************************************
325        Obsolete 1-based subroutine. See RMatrixQRUnpackQ for 0-based replacement.
326        *************************************************************************/
327        public static void unpackqfromqr(ref double[,] a,
328            int m,
329            int n,
330            ref double[] tau,
331            int qcolumns,
332            ref double[,] q)
333        {
334            int i = 0;
335            int j = 0;
336            int k = 0;
337            int minmn = 0;
338            double[] v = new double[0];
339            double[] work = new double[0];
340            int vm = 0;
341            int i_ = 0;
342            int i1_ = 0;
343
344            System.Diagnostics.Debug.Assert(qcolumns<=m, "UnpackQFromQR: QColumns>M!");
345            if( m==0 | n==0 | qcolumns==0 )
346            {
347                return;
348            }
349           
350            //
351            // init
352            //
353            minmn = Math.Min(m, n);
354            k = Math.Min(minmn, qcolumns);
355            q = new double[m+1, qcolumns+1];
356            v = new double[m+1];
357            work = new double[qcolumns+1];
358            for(i=1; i<=m; i++)
359            {
360                for(j=1; j<=qcolumns; j++)
361                {
362                    if( i==j )
363                    {
364                        q[i,j] = 1;
365                    }
366                    else
367                    {
368                        q[i,j] = 0;
369                    }
370                }
371            }
372           
373            //
374            // unpack Q
375            //
376            for(i=k; i>=1; i--)
377            {
378               
379                //
380                // Apply H(i)
381                //
382                vm = m-i+1;
383                i1_ = (i) - (1);
384                for(i_=1; i_<=vm;i_++)
385                {
386                    v[i_] = a[i_+i1_,i];
387                }
388                v[1] = 1;
389                reflections.applyreflectionfromtheleft(ref q, tau[i], ref v, i, m, 1, qcolumns, ref work);
390            }
391        }
392
393
394        /*************************************************************************
395        Obsolete 1-based subroutine. See RMatrixQR for 0-based replacement.
396        *************************************************************************/
397        public static void qrdecompositionunpacked(double[,] a,
398            int m,
399            int n,
400            ref double[,] q,
401            ref double[,] r)
402        {
403            int i = 0;
404            int k = 0;
405            double[] tau = new double[0];
406            double[] work = new double[0];
407            double[] v = new double[0];
408            int i_ = 0;
409
410            a = (double[,])a.Clone();
411
412            k = Math.Min(m, n);
413            if( n<=0 )
414            {
415                return;
416            }
417            work = new double[m+1];
418            v = new double[m+1];
419            q = new double[m+1, m+1];
420            r = new double[m+1, n+1];
421           
422            //
423            // QRDecomposition
424            //
425            qrdecomposition(ref a, m, n, ref tau);
426           
427            //
428            // R
429            //
430            for(i=1; i<=n; i++)
431            {
432                r[1,i] = 0;
433            }
434            for(i=2; i<=m; i++)
435            {
436                for(i_=1; i_<=n;i_++)
437                {
438                    r[i,i_] = r[1,i_];
439                }
440            }
441            for(i=1; i<=k; i++)
442            {
443                for(i_=i; i_<=n;i_++)
444                {
445                    r[i,i_] = a[i,i_];
446                }
447            }
448           
449            //
450            // Q
451            //
452            unpackqfromqr(ref a, m, n, ref tau, m, ref q);
453        }
454    }
455}
Note: See TracBrowser for help on using the repository browser.