source: git/Singular/svd/reflections.h @ 7fb70f

spielwiese
Last change on this file since 7fb70f was 7fb70f, checked in by Hans Schoenemann <hannes@…>, 6 years ago
add (from 3-1): system("svd",m)
  • Property mode set to 100644
File size: 10.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 _reflections_h
40#define _reflections_h
41
42#include "ap.h"
43#include "amp.h"
44namespace reflections
45{
46    template<unsigned int Precision>
47    void generatereflection(ap::template_1d_array< amp::ampf<Precision> >& x,
48        int n,
49        amp::ampf<Precision>& tau);
50    template<unsigned int Precision>
51    void applyreflectionfromtheleft(ap::template_2d_array< amp::ampf<Precision> >& c,
52        amp::ampf<Precision> tau,
53        const ap::template_1d_array< amp::ampf<Precision> >& v,
54        int m1,
55        int m2,
56        int n1,
57        int n2,
58        ap::template_1d_array< amp::ampf<Precision> >& work);
59    template<unsigned int Precision>
60    void applyreflectionfromtheright(ap::template_2d_array< amp::ampf<Precision> >& c,
61        amp::ampf<Precision> tau,
62        const ap::template_1d_array< amp::ampf<Precision> >& v,
63        int m1,
64        int m2,
65        int n1,
66        int n2,
67        ap::template_1d_array< amp::ampf<Precision> >& work);
68
69
70    /*************************************************************************
71    Generation of an elementary reflection transformation
72
73    The subroutine generates elementary reflection H of order N, so that, for
74    a given X, the following equality holds true:
75
76        ( X(1) )   ( Beta )
77    H * (  ..  ) = (  0   )
78        ( X(n) )   (  0   )
79
80    where
81                  ( V(1) )
82    H = 1 - Tau * (  ..  ) * ( V(1), ..., V(n) )
83                  ( V(n) )
84
85    where the first component of vector V equals 1.
86
87    Input parameters:
88        X   -   vector. Array whose index ranges within [1..N].
89        N   -   reflection order.
90
91    Output parameters:
92        X   -   components from 2 to N are replaced with vector V.
93                The first component is replaced with parameter Beta.
94        Tau -   scalar value Tau. If X is a null vector, Tau equals 0,
95                otherwise 1 <= Tau <= 2.
96
97    This subroutine is the modification of the DLARFG subroutines from
98    the LAPACK library. It has a similar functionality except for the
99    fact that it doesn’t handle errors when the intermediate results
100    cause an overflow.
101
102
103    MODIFICATIONS:
104        24.12.2005 sign(Alpha) was replaced with an analogous to the Fortran SIGN code.
105
106      -- LAPACK auxiliary routine (version 3.0) --
107         Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
108         Courant Institute, Argonne National Lab, and Rice University
109         September 30, 1994
110    *************************************************************************/
111    template<unsigned int Precision>
112    void generatereflection(ap::template_1d_array< amp::ampf<Precision> >& x,
113        int n,
114        amp::ampf<Precision>& tau)
115    {
116        int j;
117        amp::ampf<Precision> alpha;
118        amp::ampf<Precision> xnorm;
119        amp::ampf<Precision> v;
120        amp::ampf<Precision> beta;
121        amp::ampf<Precision> mx;
122
123
124       
125        //
126        // Executable Statements ..
127        //
128        if( n<=1 )
129        {
130            tau = 0;
131            return;
132        }
133       
134        //
135        // XNORM = DNRM2( N-1, X, INCX )
136        //
137        alpha = x(1);
138        mx = 0;
139        for(j=2; j<=n; j++)
140        {
141            mx = amp::maximum<Precision>(amp::abs<Precision>(x(j)), mx);
142        }
143        xnorm = 0;
144        if( mx!=0 )
145        {
146            for(j=2; j<=n; j++)
147            {
148                xnorm = xnorm+amp::sqr<Precision>(x(j)/mx);
149            }
150            xnorm = amp::sqrt<Precision>(xnorm)*mx;
151        }
152        if( xnorm==0 )
153        {
154           
155            //
156            // H  =  I
157            //
158            tau = 0;
159            return;
160        }
161       
162        //
163        // general case
164        //
165        mx = amp::maximum<Precision>(amp::abs<Precision>(alpha), amp::abs<Precision>(xnorm));
166        beta = -mx*amp::sqrt<Precision>(amp::sqr<Precision>(alpha/mx)+amp::sqr<Precision>(xnorm/mx));
167        if( alpha<0 )
168        {
169            beta = -beta;
170        }
171        tau = (beta-alpha)/beta;
172        v = 1/(alpha-beta);
173        ap::vmul(x.getvector(2, n), v);
174        x(1) = beta;
175    }
176
177
178    /*************************************************************************
179    Application of an elementary reflection to a rectangular matrix of size MxN
180
181    The algorithm pre-multiplies the matrix by an elementary reflection transformation
182    which is given by column V and scalar Tau (see the description of the
183    GenerateReflection procedure). Not the whole matrix but only a part of it
184    is transformed (rows from M1 to M2, columns from N1 to N2). Only the elements
185    of this submatrix are changed.
186
187    Input parameters:
188        C       -   matrix to be transformed.
189        Tau     -   scalar defining the transformation.
190        V       -   column defining the transformation.
191                    Array whose index ranges within [1..M2-M1+1].
192        M1, M2  -   range of rows to be transformed.
193        N1, N2  -   range of columns to be transformed.
194        WORK    -   working array whose indexes goes from N1 to N2.
195
196    Output parameters:
197        C       -   the result of multiplying the input matrix C by the
198                    transformation matrix which is given by Tau and V.
199                    If N1>N2 or M1>M2, C is not modified.
200
201      -- LAPACK auxiliary routine (version 3.0) --
202         Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
203         Courant Institute, Argonne National Lab, and Rice University
204         September 30, 1994
205    *************************************************************************/
206    template<unsigned int Precision>
207    void applyreflectionfromtheleft(ap::template_2d_array< amp::ampf<Precision> >& c,
208        amp::ampf<Precision> tau,
209        const ap::template_1d_array< amp::ampf<Precision> >& v,
210        int m1,
211        int m2,
212        int n1,
213        int n2,
214        ap::template_1d_array< amp::ampf<Precision> >& work)
215    {
216        amp::ampf<Precision> t;
217        int i;
218        int vm;
219
220
221        if( tau==0 || n1>n2 || m1>m2 )
222        {
223            return;
224        }
225       
226        //
227        // w := C' * v
228        //
229        vm = m2-m1+1;
230        for(i=n1; i<=n2; i++)
231        {
232            work(i) = 0;
233        }
234        for(i=m1; i<=m2; i++)
235        {
236            t = v(i+1-m1);
237            ap::vadd(work.getvector(n1, n2), c.getrow(i, n1, n2), t);
238        }
239       
240        //
241        // C := C - tau * v * w'
242        //
243        for(i=m1; i<=m2; i++)
244        {
245            t = v(i-m1+1)*tau;
246            ap::vsub(c.getrow(i, n1, n2), work.getvector(n1, n2), t);
247        }
248    }
249
250
251    /*************************************************************************
252    Application of an elementary reflection to a rectangular matrix of size MxN
253
254    The algorithm post-multiplies the matrix by an elementary reflection transformation
255    which is given by column V and scalar Tau (see the description of the
256    GenerateReflection procedure). Not the whole matrix but only a part of it
257    is transformed (rows from M1 to M2, columns from N1 to N2). Only the
258    elements of this submatrix are changed.
259
260    Input parameters:
261        C       -   matrix to be transformed.
262        Tau     -   scalar defining the transformation.
263        V       -   column defining the transformation.
264                    Array whose index ranges within [1..N2-N1+1].
265        M1, M2  -   range of rows to be transformed.
266        N1, N2  -   range of columns to be transformed.
267        WORK    -   working array whose indexes goes from M1 to M2.
268
269    Output parameters:
270        C       -   the result of multiplying the input matrix C by the
271                    transformation matrix which is given by Tau and V.
272                    If N1>N2 or M1>M2, C is not modified.
273
274      -- LAPACK auxiliary routine (version 3.0) --
275         Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
276         Courant Institute, Argonne National Lab, and Rice University
277         September 30, 1994
278    *************************************************************************/
279    template<unsigned int Precision>
280    void applyreflectionfromtheright(ap::template_2d_array< amp::ampf<Precision> >& c,
281        amp::ampf<Precision> tau,
282        const ap::template_1d_array< amp::ampf<Precision> >& v,
283        int m1,
284        int m2,
285        int n1,
286        int n2,
287        ap::template_1d_array< amp::ampf<Precision> >& work)
288    {
289        amp::ampf<Precision> t;
290        int i;
291        int vm;
292
293
294        if( tau==0 || n1>n2 || m1>m2 )
295        {
296            return;
297        }
298       
299        //
300        // w := C * v
301        //
302        vm = n2-n1+1;
303        for(i=m1; i<=m2; i++)
304        {
305            t = ap::vdotproduct(c.getrow(i, n1, n2), v.getvector(1, vm));
306            work(i) = t;
307        }
308       
309        //
310        // C := C - w * v'
311        //
312        for(i=m1; i<=m2; i++)
313        {
314            t = work(i)*tau;
315            ap::vsub(c.getrow(i, n1, n2), v.getvector(1, vm), t);
316        }
317    }
318} // namespace
319
320#endif
Note: See TracBrowser for help on using the repository browser.