"Fossies" - the Fresh Open Source Software Archive

Member "singular-4.2.1/Singular/svd/reflections.h" (9 Jun 2021, 10758 Bytes) of package /linux/misc/singular-4.2.1.tar.gz:


As a special service "Fossies" has tried to format the requested source page into HTML format using (guessed) C and C++ source code syntax highlighting (style: standard) with prefixed line numbers and code folding option. Alternatively you can here view or download the uninterpreted source code file. For more information about "reflections.h" see the Fossies "Dox" file reference documentation and the latest Fossies "Diffs" side-by-side code changes report: 4.2.0p3_vs_4.2.1.

    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 #ifndef _reflections_h
   40 #define _reflections_h
   41 
   42 #include "ap.h"
   43 #include "amp.h"
   44 namespace 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