source: git/Singular/svd/qr.h @ aad4ca4

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