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