1 | /*************************************************************************
|
---|
2 | Copyright (c) 1992-2007 The University of Tennessee. All rights reserved.
|
---|
3 |
|
---|
4 | Contributors:
|
---|
5 | * Sergey Bochkanov (ALGLIB project). Translation from FORTRAN to
|
---|
6 | pseudocode.
|
---|
7 |
|
---|
8 | See subroutines comments for additional copyrights.
|
---|
9 |
|
---|
10 | Redistribution and use in source and binary forms, with or without
|
---|
11 | modification, are permitted provided that the following conditions are
|
---|
12 | met:
|
---|
13 |
|
---|
14 | - Redistributions of source code must retain the above copyright
|
---|
15 | notice, this list of conditions and the following disclaimer.
|
---|
16 |
|
---|
17 | - Redistributions in binary form must reproduce the above copyright
|
---|
18 | notice, this list of conditions and the following disclaimer listed
|
---|
19 | in this license in the documentation and/or other materials
|
---|
20 | provided with the distribution.
|
---|
21 |
|
---|
22 | - Neither the name of the copyright holders nor the names of its
|
---|
23 | contributors may be used to endorse or promote products derived from
|
---|
24 | this software without specific prior written permission.
|
---|
25 |
|
---|
26 | THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
|
---|
27 | "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
|
---|
28 | LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
|
---|
29 | A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT
|
---|
30 | OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
|
---|
31 | SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
|
---|
32 | LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
|
---|
33 | DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
|
---|
34 | THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
|
---|
35 | (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
|
---|
36 | OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
|
---|
37 | *************************************************************************/
|
---|
38 |
|
---|
39 | using System;
|
---|
40 |
|
---|
41 | class qr
|
---|
42 | {
|
---|
43 | /*************************************************************************
|
---|
44 | QR decomposition of a rectangular matrix of size MxN
|
---|
45 |
|
---|
46 | Input parameters:
|
---|
47 | A - matrix A whose indexes range within [0..M-1, 0..N-1].
|
---|
48 | M - number of rows in matrix A.
|
---|
49 | N - number of columns in matrix A.
|
---|
50 |
|
---|
51 | Output parameters:
|
---|
52 | A - matrices Q and R in compact form (see below).
|
---|
53 | Tau - array of scalar factors which are used to form
|
---|
54 | matrix Q. Array whose index ranges within [0.. Min(M-1,N-1)].
|
---|
55 |
|
---|
56 | Matrix A is represented as A = QR, where Q is an orthogonal matrix of size
|
---|
57 | MxM, R - upper triangular (or upper trapezoid) matrix of size M x N.
|
---|
58 |
|
---|
59 | The elements of matrix R are located on and above the main diagonal of
|
---|
60 | matrix A. The elements which are located in Tau array and below the main
|
---|
61 | diagonal of matrix A are used to form matrix Q as follows:
|
---|
62 |
|
---|
63 | Matrix Q is represented as a product of elementary reflections
|
---|
64 |
|
---|
65 | Q = H(0)*H(2)*...*H(k-1),
|
---|
66 |
|
---|
67 | where k = min(m,n), and each H(i) is in the form
|
---|
68 |
|
---|
69 | H(i) = 1 - tau * v * (v^T)
|
---|
70 |
|
---|
71 | where tau is a scalar stored in Tau[I]; v - real vector,
|
---|
72 | so that v(0:i-1) = 0, v(i) = 1, v(i+1:m-1) stored in A(i+1:m-1,i).
|
---|
73 |
|
---|
74 | -- LAPACK routine (version 3.0) --
|
---|
75 | Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
|
---|
76 | Courant Institute, Argonne National Lab, and Rice University
|
---|
77 | February 29, 1992.
|
---|
78 | Translation from FORTRAN to pseudocode (AlgoPascal)
|
---|
79 | by Sergey Bochkanov, ALGLIB project, 2005-2007.
|
---|
80 | *************************************************************************/
|
---|
81 | public static void rmatrixqr(ref double[,] a,
|
---|
82 | int m,
|
---|
83 | int n,
|
---|
84 | ref double[] tau)
|
---|
85 | {
|
---|
86 | double[] work = new double[0];
|
---|
87 | double[] t = new double[0];
|
---|
88 | int i = 0;
|
---|
89 | int k = 0;
|
---|
90 | int minmn = 0;
|
---|
91 | double tmp = 0;
|
---|
92 | int i_ = 0;
|
---|
93 | int i1_ = 0;
|
---|
94 |
|
---|
95 | if( m<=0 | n<=0 )
|
---|
96 | {
|
---|
97 | return;
|
---|
98 | }
|
---|
99 | minmn = Math.Min(m, n);
|
---|
100 | work = new double[n-1+1];
|
---|
101 | t = new double[m+1];
|
---|
102 | tau = new double[minmn-1+1];
|
---|
103 |
|
---|
104 | //
|
---|
105 | // Test the input arguments
|
---|
106 | //
|
---|
107 | k = minmn;
|
---|
108 | for(i=0; i<=k-1; i++)
|
---|
109 | {
|
---|
110 |
|
---|
111 | //
|
---|
112 | // Generate elementary reflector H(i) to annihilate A(i+1:m,i)
|
---|
113 | //
|
---|
114 | i1_ = (i) - (1);
|
---|
115 | for(i_=1; i_<=m-i;i_++)
|
---|
116 | {
|
---|
117 | t[i_] = a[i_+i1_,i];
|
---|
118 | }
|
---|
119 | reflections.generatereflection(ref t, m-i, ref tmp);
|
---|
120 | tau[i] = tmp;
|
---|
121 | i1_ = (1) - (i);
|
---|
122 | for(i_=i; i_<=m-1;i_++)
|
---|
123 | {
|
---|
124 | a[i_,i] = t[i_+i1_];
|
---|
125 | }
|
---|
126 | t[1] = 1;
|
---|
127 | if( i<n )
|
---|
128 | {
|
---|
129 |
|
---|
130 | //
|
---|
131 | // Apply H(i) to A(i:m-1,i+1:n-1) from the left
|
---|
132 | //
|
---|
133 | reflections.applyreflectionfromtheleft(ref a, tau[i], ref t, i, m-1, i+1, n-1, ref work);
|
---|
134 | }
|
---|
135 | }
|
---|
136 | }
|
---|
137 |
|
---|
138 |
|
---|
139 | /*************************************************************************
|
---|
140 | Partial unpacking of matrix Q from the QR decomposition of a matrix A
|
---|
141 |
|
---|
142 | Input parameters:
|
---|
143 | A - matrices Q and R in compact form.
|
---|
144 | Output of RMatrixQR subroutine.
|
---|
145 | M - number of rows in given matrix A. M>=0.
|
---|
146 | N - number of columns in given matrix A. N>=0.
|
---|
147 | Tau - scalar factors which are used to form Q.
|
---|
148 | Output of the RMatrixQR subroutine.
|
---|
149 | QColumns - required number of columns of matrix Q. M>=QColumns>=0.
|
---|
150 |
|
---|
151 | Output parameters:
|
---|
152 | Q - first QColumns columns of matrix Q.
|
---|
153 | Array whose indexes range within [0..M-1, 0..QColumns-1].
|
---|
154 | If QColumns=0, the array remains unchanged.
|
---|
155 |
|
---|
156 | -- ALGLIB --
|
---|
157 | Copyright 2005 by Bochkanov Sergey
|
---|
158 | *************************************************************************/
|
---|
159 | public static void rmatrixqrunpackq(ref double[,] a,
|
---|
160 | int m,
|
---|
161 | int n,
|
---|
162 | ref double[] tau,
|
---|
163 | int qcolumns,
|
---|
164 | ref double[,] q)
|
---|
165 | {
|
---|
166 | int i = 0;
|
---|
167 | int j = 0;
|
---|
168 | int k = 0;
|
---|
169 | int minmn = 0;
|
---|
170 | double[] v = new double[0];
|
---|
171 | double[] work = new double[0];
|
---|
172 | int i_ = 0;
|
---|
173 | int i1_ = 0;
|
---|
174 |
|
---|
175 | System.Diagnostics.Debug.Assert(qcolumns<=m, "UnpackQFromQR: QColumns>M!");
|
---|
176 | if( m<=0 | n<=0 | qcolumns<=0 )
|
---|
177 | {
|
---|
178 | return;
|
---|
179 | }
|
---|
180 |
|
---|
181 | //
|
---|
182 | // init
|
---|
183 | //
|
---|
184 | minmn = Math.Min(m, n);
|
---|
185 | k = Math.Min(minmn, qcolumns);
|
---|
186 | q = new double[m-1+1, qcolumns-1+1];
|
---|
187 | v = new double[m+1];
|
---|
188 | work = new double[qcolumns-1+1];
|
---|
189 | for(i=0; i<=m-1; i++)
|
---|
190 | {
|
---|
191 | for(j=0; j<=qcolumns-1; j++)
|
---|
192 | {
|
---|
193 | if( i==j )
|
---|
194 | {
|
---|
195 | q[i,j] = 1;
|
---|
196 | }
|
---|
197 | else
|
---|
198 | {
|
---|
199 | q[i,j] = 0;
|
---|
200 | }
|
---|
201 | }
|
---|
202 | }
|
---|
203 |
|
---|
204 | //
|
---|
205 | // unpack Q
|
---|
206 | //
|
---|
207 | for(i=k-1; i>=0; i--)
|
---|
208 | {
|
---|
209 |
|
---|
210 | //
|
---|
211 | // Apply H(i)
|
---|
212 | //
|
---|
213 | i1_ = (i) - (1);
|
---|
214 | for(i_=1; i_<=m-i;i_++)
|
---|
215 | {
|
---|
216 | v[i_] = a[i_+i1_,i];
|
---|
217 | }
|
---|
218 | v[1] = 1;
|
---|
219 | reflections.applyreflectionfromtheleft(ref q, tau[i], ref v, i, m-1, 0, qcolumns-1, ref work);
|
---|
220 | }
|
---|
221 | }
|
---|
222 |
|
---|
223 |
|
---|
224 | /*************************************************************************
|
---|
225 | Unpacking of matrix R from the QR decomposition of a matrix A
|
---|
226 |
|
---|
227 | Input parameters:
|
---|
228 | A - matrices Q and R in compact form.
|
---|
229 | Output of RMatrixQR subroutine.
|
---|
230 | M - number of rows in given matrix A. M>=0.
|
---|
231 | N - number of columns in given matrix A. N>=0.
|
---|
232 |
|
---|
233 | Output parameters:
|
---|
234 | R - matrix R, array[0..M-1, 0..N-1].
|
---|
235 |
|
---|
236 | -- ALGLIB --
|
---|
237 | Copyright 2005 by Bochkanov Sergey
|
---|
238 | *************************************************************************/
|
---|
239 | public static void rmatrixqrunpackr(ref double[,] a,
|
---|
240 | int m,
|
---|
241 | int n,
|
---|
242 | ref double[,] r)
|
---|
243 | {
|
---|
244 | int i = 0;
|
---|
245 | int k = 0;
|
---|
246 | int i_ = 0;
|
---|
247 |
|
---|
248 | if( m<=0 | n<=0 )
|
---|
249 | {
|
---|
250 | return;
|
---|
251 | }
|
---|
252 | k = Math.Min(m, n);
|
---|
253 | r = new double[m-1+1, n-1+1];
|
---|
254 | for(i=0; i<=n-1; i++)
|
---|
255 | {
|
---|
256 | r[0,i] = 0;
|
---|
257 | }
|
---|
258 | for(i=1; i<=m-1; i++)
|
---|
259 | {
|
---|
260 | for(i_=0; i_<=n-1;i_++)
|
---|
261 | {
|
---|
262 | r[i,i_] = r[0,i_];
|
---|
263 | }
|
---|
264 | }
|
---|
265 | for(i=0; i<=k-1; i++)
|
---|
266 | {
|
---|
267 | for(i_=i; i_<=n-1;i_++)
|
---|
268 | {
|
---|
269 | r[i,i_] = a[i,i_];
|
---|
270 | }
|
---|
271 | }
|
---|
272 | }
|
---|
273 |
|
---|
274 |
|
---|
275 | /*************************************************************************
|
---|
276 | Obsolete 1-based subroutine. See RMatrixQR for 0-based replacement.
|
---|
277 | *************************************************************************/
|
---|
278 | public static void qrdecomposition(ref double[,] a,
|
---|
279 | int m,
|
---|
280 | int n,
|
---|
281 | ref double[] tau)
|
---|
282 | {
|
---|
283 | double[] work = new double[0];
|
---|
284 | double[] t = new double[0];
|
---|
285 | int i = 0;
|
---|
286 | int k = 0;
|
---|
287 | int mmip1 = 0;
|
---|
288 | int minmn = 0;
|
---|
289 | double tmp = 0;
|
---|
290 | int i_ = 0;
|
---|
291 | int i1_ = 0;
|
---|
292 |
|
---|
293 | minmn = Math.Min(m, n);
|
---|
294 | work = new double[n+1];
|
---|
295 | t = new double[m+1];
|
---|
296 | tau = new double[minmn+1];
|
---|
297 |
|
---|
298 | //
|
---|
299 | // Test the input arguments
|
---|
300 | //
|
---|
301 | k = Math.Min(m, n);
|
---|
302 | for(i=1; i<=k; i++)
|
---|
303 | {
|
---|
304 |
|
---|
305 | //
|
---|
306 | // Generate elementary reflector H(i) to annihilate A(i+1:m,i)
|
---|
307 | //
|
---|
308 | mmip1 = m-i+1;
|
---|
309 | i1_ = (i) - (1);
|
---|
310 | for(i_=1; i_<=mmip1;i_++)
|
---|
311 | {
|
---|
312 | t[i_] = a[i_+i1_,i];
|
---|
313 | }
|
---|
314 | reflections.generatereflection(ref t, mmip1, ref tmp);
|
---|
315 | tau[i] = tmp;
|
---|
316 | i1_ = (1) - (i);
|
---|
317 | for(i_=i; i_<=m;i_++)
|
---|
318 | {
|
---|
319 | a[i_,i] = t[i_+i1_];
|
---|
320 | }
|
---|
321 | t[1] = 1;
|
---|
322 | if( i<n )
|
---|
323 | {
|
---|
324 |
|
---|
325 | //
|
---|
326 | // Apply H(i) to A(i:m,i+1:n) from the left
|
---|
327 | //
|
---|
328 | reflections.applyreflectionfromtheleft(ref a, tau[i], ref t, i, m, i+1, n, ref work);
|
---|
329 | }
|
---|
330 | }
|
---|
331 | }
|
---|
332 |
|
---|
333 |
|
---|
334 | /*************************************************************************
|
---|
335 | Obsolete 1-based subroutine. See RMatrixQRUnpackQ for 0-based replacement.
|
---|
336 | *************************************************************************/
|
---|
337 | public static void unpackqfromqr(ref double[,] a,
|
---|
338 | int m,
|
---|
339 | int n,
|
---|
340 | ref double[] tau,
|
---|
341 | int qcolumns,
|
---|
342 | ref double[,] q)
|
---|
343 | {
|
---|
344 | int i = 0;
|
---|
345 | int j = 0;
|
---|
346 | int k = 0;
|
---|
347 | int minmn = 0;
|
---|
348 | double[] v = new double[0];
|
---|
349 | double[] work = new double[0];
|
---|
350 | int vm = 0;
|
---|
351 | int i_ = 0;
|
---|
352 | int i1_ = 0;
|
---|
353 |
|
---|
354 | System.Diagnostics.Debug.Assert(qcolumns<=m, "UnpackQFromQR: QColumns>M!");
|
---|
355 | if( m==0 | n==0 | qcolumns==0 )
|
---|
356 | {
|
---|
357 | return;
|
---|
358 | }
|
---|
359 |
|
---|
360 | //
|
---|
361 | // init
|
---|
362 | //
|
---|
363 | minmn = Math.Min(m, n);
|
---|
364 | k = Math.Min(minmn, qcolumns);
|
---|
365 | q = new double[m+1, qcolumns+1];
|
---|
366 | v = new double[m+1];
|
---|
367 | work = new double[qcolumns+1];
|
---|
368 | for(i=1; i<=m; i++)
|
---|
369 | {
|
---|
370 | for(j=1; j<=qcolumns; j++)
|
---|
371 | {
|
---|
372 | if( i==j )
|
---|
373 | {
|
---|
374 | q[i,j] = 1;
|
---|
375 | }
|
---|
376 | else
|
---|
377 | {
|
---|
378 | q[i,j] = 0;
|
---|
379 | }
|
---|
380 | }
|
---|
381 | }
|
---|
382 |
|
---|
383 | //
|
---|
384 | // unpack Q
|
---|
385 | //
|
---|
386 | for(i=k; i>=1; i--)
|
---|
387 | {
|
---|
388 |
|
---|
389 | //
|
---|
390 | // Apply H(i)
|
---|
391 | //
|
---|
392 | vm = m-i+1;
|
---|
393 | i1_ = (i) - (1);
|
---|
394 | for(i_=1; i_<=vm;i_++)
|
---|
395 | {
|
---|
396 | v[i_] = a[i_+i1_,i];
|
---|
397 | }
|
---|
398 | v[1] = 1;
|
---|
399 | reflections.applyreflectionfromtheleft(ref q, tau[i], ref v, i, m, 1, qcolumns, ref work);
|
---|
400 | }
|
---|
401 | }
|
---|
402 |
|
---|
403 |
|
---|
404 | /*************************************************************************
|
---|
405 | Obsolete 1-based subroutine. See RMatrixQR for 0-based replacement.
|
---|
406 | *************************************************************************/
|
---|
407 | public static void qrdecompositionunpacked(double[,] a,
|
---|
408 | int m,
|
---|
409 | int n,
|
---|
410 | ref double[,] q,
|
---|
411 | ref double[,] r)
|
---|
412 | {
|
---|
413 | int i = 0;
|
---|
414 | int k = 0;
|
---|
415 | double[] tau = new double[0];
|
---|
416 | double[] work = new double[0];
|
---|
417 | double[] v = new double[0];
|
---|
418 | int i_ = 0;
|
---|
419 |
|
---|
420 | a = (double[,])a.Clone();
|
---|
421 |
|
---|
422 | k = Math.Min(m, n);
|
---|
423 | if( n<=0 )
|
---|
424 | {
|
---|
425 | return;
|
---|
426 | }
|
---|
427 | work = new double[m+1];
|
---|
428 | v = new double[m+1];
|
---|
429 | q = new double[m+1, m+1];
|
---|
430 | r = new double[m+1, n+1];
|
---|
431 |
|
---|
432 | //
|
---|
433 | // QRDecomposition
|
---|
434 | //
|
---|
435 | qrdecomposition(ref a, m, n, ref tau);
|
---|
436 |
|
---|
437 | //
|
---|
438 | // R
|
---|
439 | //
|
---|
440 | for(i=1; i<=n; i++)
|
---|
441 | {
|
---|
442 | r[1,i] = 0;
|
---|
443 | }
|
---|
444 | for(i=2; i<=m; i++)
|
---|
445 | {
|
---|
446 | for(i_=1; i_<=n;i_++)
|
---|
447 | {
|
---|
448 | r[i,i_] = r[1,i_];
|
---|
449 | }
|
---|
450 | }
|
---|
451 | for(i=1; i<=k; i++)
|
---|
452 | {
|
---|
453 | for(i_=i; i_<=n;i_++)
|
---|
454 | {
|
---|
455 | r[i,i_] = a[i,i_];
|
---|
456 | }
|
---|
457 | }
|
---|
458 |
|
---|
459 | //
|
---|
460 | // Q
|
---|
461 | //
|
---|
462 | unpackqfromqr(ref a, m, n, ref tau, m, ref q);
|
---|
463 | }
|
---|
464 | }
|
---|