Free cookie consent management tool by TermsFeed Policy Generator

source: trunk/sources/ALGLIB/qr.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: 13.7 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        public static void qrdecomposition(ref double[,] a,
266            int m,
267            int n,
268            ref double[] tau)
269        {
270            double[] work = new double[0];
271            double[] t = new double[0];
272            int i = 0;
273            int k = 0;
274            int mmip1 = 0;
275            int minmn = 0;
276            double tmp = 0;
277            int i_ = 0;
278            int i1_ = 0;
279
280            minmn = Math.Min(m, n);
281            work = new double[n+1];
282            t = new double[m+1];
283            tau = new double[minmn+1];
284           
285            //
286            // Test the input arguments
287            //
288            k = Math.Min(m, n);
289            for(i=1; i<=k; i++)
290            {
291               
292                //
293                // Generate elementary reflector H(i) to annihilate A(i+1:m,i)
294                //
295                mmip1 = m-i+1;
296                i1_ = (i) - (1);
297                for(i_=1; i_<=mmip1;i_++)
298                {
299                    t[i_] = a[i_+i1_,i];
300                }
301                reflections.generatereflection(ref t, mmip1, ref tmp);
302                tau[i] = tmp;
303                i1_ = (1) - (i);
304                for(i_=i; i_<=m;i_++)
305                {
306                    a[i_,i] = t[i_+i1_];
307                }
308                t[1] = 1;
309                if( i<n )
310                {
311                   
312                    //
313                    // Apply H(i) to A(i:m,i+1:n) from the left
314                    //
315                    reflections.applyreflectionfromtheleft(ref a, tau[i], ref t, i, m, i+1, n, ref work);
316                }
317            }
318        }
319
320
321        public static void unpackqfromqr(ref double[,] a,
322            int m,
323            int n,
324            ref double[] tau,
325            int qcolumns,
326            ref double[,] q)
327        {
328            int i = 0;
329            int j = 0;
330            int k = 0;
331            int minmn = 0;
332            double[] v = new double[0];
333            double[] work = new double[0];
334            int vm = 0;
335            int i_ = 0;
336            int i1_ = 0;
337
338            System.Diagnostics.Debug.Assert(qcolumns<=m, "UnpackQFromQR: QColumns>M!");
339            if( m==0 | n==0 | qcolumns==0 )
340            {
341                return;
342            }
343           
344            //
345            // init
346            //
347            minmn = Math.Min(m, n);
348            k = Math.Min(minmn, qcolumns);
349            q = new double[m+1, qcolumns+1];
350            v = new double[m+1];
351            work = new double[qcolumns+1];
352            for(i=1; i<=m; i++)
353            {
354                for(j=1; j<=qcolumns; j++)
355                {
356                    if( i==j )
357                    {
358                        q[i,j] = 1;
359                    }
360                    else
361                    {
362                        q[i,j] = 0;
363                    }
364                }
365            }
366           
367            //
368            // unpack Q
369            //
370            for(i=k; i>=1; i--)
371            {
372               
373                //
374                // Apply H(i)
375                //
376                vm = m-i+1;
377                i1_ = (i) - (1);
378                for(i_=1; i_<=vm;i_++)
379                {
380                    v[i_] = a[i_+i1_,i];
381                }
382                v[1] = 1;
383                reflections.applyreflectionfromtheleft(ref q, tau[i], ref v, i, m, 1, qcolumns, ref work);
384            }
385        }
386
387
388        public static void qrdecompositionunpacked(double[,] a,
389            int m,
390            int n,
391            ref double[,] q,
392            ref double[,] r)
393        {
394            int i = 0;
395            int k = 0;
396            double[] tau = new double[0];
397            double[] work = new double[0];
398            double[] v = new double[0];
399            int i_ = 0;
400
401            a = (double[,])a.Clone();
402
403            k = Math.Min(m, n);
404            if( n<=0 )
405            {
406                return;
407            }
408            work = new double[m+1];
409            v = new double[m+1];
410            q = new double[m+1, m+1];
411            r = new double[m+1, n+1];
412           
413            //
414            // QRDecomposition
415            //
416            qrdecomposition(ref a, m, n, ref tau);
417           
418            //
419            // R
420            //
421            for(i=1; i<=n; i++)
422            {
423                r[1,i] = 0;
424            }
425            for(i=2; i<=m; i++)
426            {
427                for(i_=1; i_<=n;i_++)
428                {
429                    r[i,i_] = r[1,i_];
430                }
431            }
432            for(i=1; i<=k; i++)
433            {
434                for(i_=i; i_<=n;i_++)
435                {
436                    r[i,i_] = a[i,i_];
437                }
438            }
439           
440            //
441            // Q
442            //
443            unpackqfromqr(ref a, m, n, ref tau, m, ref q);
444        }
445    }
446}
Note: See TracBrowser for help on using the repository browser.