source: git/Singular/svd/rotations.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.3 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 _rotations_h
40#define _rotations_h
41
42#include "ap.h"
43#include "amp.h"
44namespace rotations
45{
46    template<unsigned int Precision>
47    void applyrotationsfromtheleft(bool isforward,
48        int m1,
49        int m2,
50        int n1,
51        int n2,
52        const ap::template_1d_array< amp::ampf<Precision> >& c,
53        const ap::template_1d_array< amp::ampf<Precision> >& s,
54        ap::template_2d_array< amp::ampf<Precision> >& a,
55        ap::template_1d_array< amp::ampf<Precision> >& work);
56    template<unsigned int Precision>
57    void applyrotationsfromtheright(bool isforward,
58        int m1,
59        int m2,
60        int n1,
61        int n2,
62        const ap::template_1d_array< amp::ampf<Precision> >& c,
63        const ap::template_1d_array< amp::ampf<Precision> >& s,
64        ap::template_2d_array< amp::ampf<Precision> >& a,
65        ap::template_1d_array< amp::ampf<Precision> >& work);
66    template<unsigned int Precision>
67    void generaterotation(amp::ampf<Precision> f,
68        amp::ampf<Precision> g,
69        amp::ampf<Precision>& cs,
70        amp::ampf<Precision>& sn,
71        amp::ampf<Precision>& r);
72
73
74    /*************************************************************************
75    Application of a sequence of  elementary rotations to a matrix
76
77    The algorithm pre-multiplies the matrix by a sequence of rotation
78    transformations which is given by arrays C and S. Depending on the value
79    of the IsForward parameter either 1 and 2, 3 and 4 and so on (if IsForward=true)
80    rows are rotated, or the rows N and N-1, N-2 and N-3 and so on, are rotated.
81
82    Not the whole matrix but only a part of it is transformed (rows from M1 to
83    M2, columns from N1 to N2). Only the elements of this submatrix are changed.
84
85    Input parameters:
86        IsForward   -   the sequence of the rotation application.
87        M1,M2       -   the range of rows to be transformed.
88        N1, N2      -   the range of columns to be transformed.
89        C,S         -   transformation coefficients.
90                        Array whose index ranges within [1..M2-M1].
91        A           -   processed matrix.
92        WORK        -   working array whose index ranges within [N1..N2].
93
94    Output parameters:
95        A           -   transformed matrix.
96
97    Utility subroutine.
98    *************************************************************************/
99    template<unsigned int Precision>
100    void applyrotationsfromtheleft(bool isforward,
101        int m1,
102        int m2,
103        int n1,
104        int n2,
105        const ap::template_1d_array< amp::ampf<Precision> >& c,
106        const ap::template_1d_array< amp::ampf<Precision> >& s,
107        ap::template_2d_array< amp::ampf<Precision> >& a,
108        ap::template_1d_array< amp::ampf<Precision> >& work)
109    {
110        int j;
111        int jp1;
112        amp::ampf<Precision> ctemp;
113        amp::ampf<Precision> stemp;
114        amp::ampf<Precision> temp;
115
116
117        if( m1>m2 || n1>n2 )
118        {
119            return;
120        }
121
122        //
123        // Form  P * A
124        //
125        if( isforward )
126        {
127            if( n1!=n2 )
128            {
129
130                //
131                // Common case: N1<>N2
132                //
133                for(j=m1; j<=m2-1; j++)
134                {
135                    ctemp = c(j-m1+1);
136                    stemp = s(j-m1+1);
137                    if( ctemp!=1 || stemp!=0 )
138                    {
139                        jp1 = j+1;
140                        ap::vmove(work.getvector(n1, n2), a.getrow(jp1, n1, n2), ctemp);
141                        ap::vsub(work.getvector(n1, n2), a.getrow(j, n1, n2), stemp);
142                        ap::vmul(a.getrow(j, n1, n2), ctemp);
143                        ap::vadd(a.getrow(j, n1, n2), a.getrow(jp1, n1, n2), stemp);
144                        ap::vmove(a.getrow(jp1, n1, n2), work.getvector(n1, n2));
145                    }
146                }
147            }
148            else
149            {
150
151                //
152                // Special case: N1=N2
153                //
154                for(j=m1; j<=m2-1; j++)
155                {
156                    ctemp = c(j-m1+1);
157                    stemp = s(j-m1+1);
158                    if( ctemp!=1 || stemp!=0 )
159                    {
160                        temp = a(j+1,n1);
161                        a(j+1,n1) = ctemp*temp-stemp*a(j,n1);
162                        a(j,n1) = stemp*temp+ctemp*a(j,n1);
163                    }
164                }
165            }
166        }
167        else
168        {
169            if( n1!=n2 )
170            {
171
172                //
173                // Common case: N1<>N2
174                //
175                for(j=m2-1; j>=m1; j--)
176                {
177                    ctemp = c(j-m1+1);
178                    stemp = s(j-m1+1);
179                    if( ctemp!=1 || stemp!=0 )
180                    {
181                        jp1 = j+1;
182                        ap::vmove(work.getvector(n1, n2), a.getrow(jp1, n1, n2), ctemp);
183                        ap::vsub(work.getvector(n1, n2), a.getrow(j, n1, n2), stemp);
184                        ap::vmul(a.getrow(j, n1, n2), ctemp);
185                        ap::vadd(a.getrow(j, n1, n2), a.getrow(jp1, n1, n2), stemp);
186                        ap::vmove(a.getrow(jp1, n1, n2), work.getvector(n1, n2));
187                    }
188                }
189            }
190            else
191            {
192
193                //
194                // Special case: N1=N2
195                //
196                for(j=m2-1; j>=m1; j--)
197                {
198                    ctemp = c(j-m1+1);
199                    stemp = s(j-m1+1);
200                    if( ctemp!=1 || stemp!=0 )
201                    {
202                        temp = a(j+1,n1);
203                        a(j+1,n1) = ctemp*temp-stemp*a(j,n1);
204                        a(j,n1) = stemp*temp+ctemp*a(j,n1);
205                    }
206                }
207            }
208        }
209    }
210
211
212    /*************************************************************************
213    Application of a sequence of  elementary rotations to a matrix
214
215    The algorithm post-multiplies the matrix by a sequence of rotation
216    transformations which is given by arrays C and S. Depending on the value
217    of the IsForward parameter either 1 and 2, 3 and 4 and so on (if IsForward=true)
218    rows are rotated, or the rows N and N-1, N-2 and N-3 and so on are rotated.
219
220    Not the whole matrix but only a part of it is transformed (rows from M1
221    to M2, columns from N1 to N2). Only the elements of this submatrix are changed.
222
223    Input parameters:
224        IsForward   -   the sequence of the rotation application.
225        M1,M2       -   the range of rows to be transformed.
226        N1, N2      -   the range of columns to be transformed.
227        C,S         -   transformation coefficients.
228                        Array whose index ranges within [1..N2-N1].
229        A           -   processed matrix.
230        WORK        -   working array whose index ranges within [M1..M2].
231
232    Output parameters:
233        A           -   transformed matrix.
234
235    Utility subroutine.
236    *************************************************************************/
237    template<unsigned int Precision>
238    void applyrotationsfromtheright(bool isforward,
239        int m1,
240        int m2,
241        int n1,
242        int n2,
243        const ap::template_1d_array< amp::ampf<Precision> >& c,
244        const ap::template_1d_array< amp::ampf<Precision> >& s,
245        ap::template_2d_array< amp::ampf<Precision> >& a,
246        ap::template_1d_array< amp::ampf<Precision> >& work)
247    {
248        int j;
249        int jp1;
250        amp::ampf<Precision> ctemp;
251        amp::ampf<Precision> stemp;
252        amp::ampf<Precision> temp;
253
254
255
256        //
257        // Form A * P'
258        //
259        if( isforward )
260        {
261            if( m1!=m2 )
262            {
263
264                //
265                // Common case: M1<>M2
266                //
267                for(j=n1; j<=n2-1; j++)
268                {
269                    ctemp = c(j-n1+1);
270                    stemp = s(j-n1+1);
271                    if( ctemp!=1 || stemp!=0 )
272                    {
273                        jp1 = j+1;
274                        ap::vmove(work.getvector(m1, m2), a.getcolumn(jp1, m1, m2), ctemp);
275                        ap::vsub(work.getvector(m1, m2), a.getcolumn(j, m1, m2), stemp);
276                        ap::vmul(a.getcolumn(j, m1, m2), ctemp);
277                        ap::vadd(a.getcolumn(j, m1, m2), a.getcolumn(jp1, m1, m2), stemp);
278                        ap::vmove(a.getcolumn(jp1, m1, m2), work.getvector(m1, m2));
279                    }
280                }
281            }
282            else
283            {
284
285                //
286                // Special case: M1=M2
287                //
288                for(j=n1; j<=n2-1; j++)
289                {
290                    ctemp = c(j-n1+1);
291                    stemp = s(j-n1+1);
292                    if( ctemp!=1 || stemp!=0 )
293                    {
294                        temp = a(m1,j+1);
295                        a(m1,j+1) = ctemp*temp-stemp*a(m1,j);
296                        a(m1,j) = stemp*temp+ctemp*a(m1,j);
297                    }
298                }
299            }
300        }
301        else
302        {
303            if( m1!=m2 )
304            {
305
306                //
307                // Common case: M1<>M2
308                //
309                for(j=n2-1; j>=n1; j--)
310                {
311                    ctemp = c(j-n1+1);
312                    stemp = s(j-n1+1);
313                    if( ctemp!=1 || stemp!=0 )
314                    {
315                        jp1 = j+1;
316                        ap::vmove(work.getvector(m1, m2), a.getcolumn(jp1, m1, m2), ctemp);
317                        ap::vsub(work.getvector(m1, m2), a.getcolumn(j, m1, m2), stemp);
318                        ap::vmul(a.getcolumn(j, m1, m2), ctemp);
319                        ap::vadd(a.getcolumn(j, m1, m2), a.getcolumn(jp1, m1, m2), stemp);
320                        ap::vmove(a.getcolumn(jp1, m1, m2), work.getvector(m1, m2));
321                    }
322                }
323            }
324            else
325            {
326
327                //
328                // Special case: M1=M2
329                //
330                for(j=n2-1; j>=n1; j--)
331                {
332                    ctemp = c(j-n1+1);
333                    stemp = s(j-n1+1);
334                    if( ctemp!=1 || stemp!=0 )
335                    {
336                        temp = a(m1,j+1);
337                        a(m1,j+1) = ctemp*temp-stemp*a(m1,j);
338                        a(m1,j) = stemp*temp+ctemp*a(m1,j);
339                    }
340                }
341            }
342        }
343    }
344
345
346    /*************************************************************************
347    The subroutine generates the elementary rotation, so that:
348
349    [  CS  SN  ]  .  [ F ]  =  [ R ]
350    [ -SN  CS  ]     [ G ]     [ 0 ]
351
352    CS**2 + SN**2 = 1
353    *************************************************************************/
354    template<unsigned int Precision>
355    void generaterotation(amp::ampf<Precision> f,
356        amp::ampf<Precision> g,
357        amp::ampf<Precision>& cs,
358        amp::ampf<Precision>& sn,
359        amp::ampf<Precision>& r)
360    {
361        amp::ampf<Precision> f1;
362        amp::ampf<Precision> g1;
363
364
365        if( g==0 )
366        {
367            cs = 1;
368            sn = 0;
369            r = f;
370        }
371        else
372        {
373            if( f==0 )
374            {
375                cs = 0;
376                sn = 1;
377                r = g;
378            }
379            else
380            {
381                f1 = f;
382                g1 = g;
383                r = amp::sqrt<Precision>(amp::sqr<Precision>(f1)+amp::sqr<Precision>(g1));
384                cs = f1/r;
385                sn = g1/r;
386                if( amp::abs<Precision>(f)>amp::abs<Precision>(g) && cs<0 )
387                {
388                    cs = -cs;
389                    sn = -sn;
390                    r = -r;
391                }
392            }
393        }
394    }
395} // namespace
396
397#endif
Note: See TracBrowser for help on using the repository browser.