Free cookie consent management tool by TermsFeed Policy Generator

Ignore:
Timestamp:
10/15/09 13:22:07 (15 years ago)
Author:
gkronber
Message:

Moved ALGLIB code into a separate plugin. #783

Location:
trunk/sources/ALGLIB
Files:
1 added
1 copied

Legend:

Unmodified
Added
Removed

  • trunk/sources/ALGLIB/reflections.cs

    r2426 r2430  
    88See subroutines comments for additional copyrights.
    99
    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.
     10>>> SOURCE LICENSE >>>
     11This program is free software; you can redistribute it and/or modify
     12it under the terms of the GNU General Public License as published by
     13the Free Software Foundation (www.fsf.org); either version 2 of the
     14License, or (at your option) any later version.
     15
     16This program is distributed in the hope that it will be useful,
     17but WITHOUT ANY WARRANTY; without even the implied warranty of
     18MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
     19GNU General Public License for more details.
     20
     21A copy of the GNU General Public License is available at
     22http://www.fsf.org/licensing/licenses
     23
     24>>> END OF LICENSE >>>
    3725*************************************************************************/
    3826
    3927using System;
    4028
    41 class reflections
     29namespace alglib
    4230{
    43     /*************************************************************************
    44     Generation of an elementary reflection transformation
    45 
    46     The subroutine generates elementary reflection H of order N, so that, for
    47     a given X, the following equality holds true:
    48 
    49         ( X(1) )   ( Beta )
    50     H * (  ..  ) = (  0   )
    51         ( X(n) )   (  0   )
    52 
    53     where
    54                   ( V(1) )
    55     H = 1 - Tau * (  ..  ) * ( V(1), ..., V(n) )
    56                   ( V(n) )
    57 
    58     where the first component of vector V equals 1.
    59 
    60     Input parameters:
    61         X   -   vector. Array whose index ranges within [1..N].
    62         N   -   reflection order.
    63 
    64     Output parameters:
    65         X   -   components from 2 to N are replaced with vector V.
    66                 The first component is replaced with parameter Beta.
    67         Tau -   scalar value Tau. If X is a null vector, Tau equals 0,
    68                 otherwise 1 <= Tau <= 2.
    69 
    70     This subroutine is the modification of the DLARFG subroutines from
    71     the LAPACK library. It has a similar functionality except for the
    72     fact that it doesn’t handle errors when the intermediate results
    73     cause an overflow.
    74 
    75 
    76     MODIFICATIONS:
    77         24.12.2005 sign(Alpha) was replaced with an analogous to the Fortran SIGN code.
    78 
    79       -- LAPACK auxiliary routine (version 3.0) --
    80          Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
    81          Courant Institute, Argonne National Lab, and Rice University
    82          September 30, 1994
    83     *************************************************************************/
    84     public static void generatereflection(ref double[] x,
    85         int n,
    86         ref double tau)
     31    public class reflections
    8732    {
    88         int j = 0;
    89         double alpha = 0;
    90         double xnorm = 0;
    91         double v = 0;
    92         double beta = 0;
    93         double mx = 0;
    94         int i_ = 0;
    95 
    96        
    97         //
    98         // Executable Statements ..
    99         //
    100         if( n<=1 )
     33        /*************************************************************************
     34        Generation of an elementary reflection transformation
     35
     36        The subroutine generates elementary reflection H of order N, so that, for
     37        a given X, the following equality holds true:
     38
     39            ( X(1) )   ( Beta )
     40        H * (  ..  ) = (  0   )
     41            ( X(n) )   (  0   )
     42
     43        where
     44                      ( V(1) )
     45        H = 1 - Tau * (  ..  ) * ( V(1), ..., V(n) )
     46                      ( V(n) )
     47
     48        where the first component of vector V equals 1.
     49
     50        Input parameters:
     51            X   -   vector. Array whose index ranges within [1..N].
     52            N   -   reflection order.
     53
     54        Output parameters:
     55            X   -   components from 2 to N are replaced with vector V.
     56                    The first component is replaced with parameter Beta.
     57            Tau -   scalar value Tau. If X is a null vector, Tau equals 0,
     58                    otherwise 1 <= Tau <= 2.
     59
     60        This subroutine is the modification of the DLARFG subroutines from
     61        the LAPACK library.
     62
     63        MODIFICATIONS:
     64            24.12.2005 sign(Alpha) was replaced with an analogous to the Fortran SIGN code.
     65
     66          -- LAPACK auxiliary routine (version 3.0) --
     67             Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
     68             Courant Institute, Argonne National Lab, and Rice University
     69             September 30, 1994
     70        *************************************************************************/
     71        public static void generatereflection(ref double[] x,
     72            int n,
     73            ref double tau)
    10174        {
    102             tau = 0;
    103             return;
     75            int j = 0;
     76            double alpha = 0;
     77            double xnorm = 0;
     78            double v = 0;
     79            double beta = 0;
     80            double mx = 0;
     81            double s = 0;
     82            int i_ = 0;
     83
     84            if( n<=1 )
     85            {
     86                tau = 0;
     87                return;
     88            }
     89           
     90            //
     91            // Scale if needed (to avoid overflow/underflow during intermediate
     92            // calculations).
     93            //
     94            mx = 0;
     95            for(j=1; j<=n; j++)
     96            {
     97                mx = Math.Max(Math.Abs(x[j]), mx);
     98            }
     99            s = 1;
     100            if( mx!=0 )
     101            {
     102                if( mx<=AP.Math.MinRealNumber/AP.Math.MachineEpsilon )
     103                {
     104                    s = AP.Math.MinRealNumber/AP.Math.MachineEpsilon;
     105                    v = 1/s;
     106                    for(i_=1; i_<=n;i_++)
     107                    {
     108                        x[i_] = v*x[i_];
     109                    }
     110                    mx = mx*v;
     111                }
     112                else
     113                {
     114                    if( mx>=AP.Math.MaxRealNumber*AP.Math.MachineEpsilon )
     115                    {
     116                        s = AP.Math.MaxRealNumber*AP.Math.MachineEpsilon;
     117                        v = 1/s;
     118                        for(i_=1; i_<=n;i_++)
     119                        {
     120                            x[i_] = v*x[i_];
     121                        }
     122                        mx = mx*v;
     123                    }
     124                }
     125            }
     126           
     127            //
     128            // XNORM = DNRM2( N-1, X, INCX )
     129            //
     130            alpha = x[1];
     131            xnorm = 0;
     132            if( mx!=0 )
     133            {
     134                for(j=2; j<=n; j++)
     135                {
     136                    xnorm = xnorm+AP.Math.Sqr(x[j]/mx);
     137                }
     138                xnorm = Math.Sqrt(xnorm)*mx;
     139            }
     140            if( xnorm==0 )
     141            {
     142               
     143                //
     144                // H  =  I
     145                //
     146                tau = 0;
     147                x[1] = x[1]*s;
     148                return;
     149            }
     150           
     151            //
     152            // general case
     153            //
     154            mx = Math.Max(Math.Abs(alpha), Math.Abs(xnorm));
     155            beta = -(mx*Math.Sqrt(AP.Math.Sqr(alpha/mx)+AP.Math.Sqr(xnorm/mx)));
     156            if( alpha<0 )
     157            {
     158                beta = -beta;
     159            }
     160            tau = (beta-alpha)/beta;
     161            v = 1/(alpha-beta);
     162            for(i_=2; i_<=n;i_++)
     163            {
     164                x[i_] = v*x[i_];
     165            }
     166            x[1] = beta;
     167           
     168            //
     169            // Scale back outputs
     170            //
     171            x[1] = x[1]*s;
    104172        }
    105        
    106         //
    107         // XNORM = DNRM2( N-1, X, INCX )
    108         //
    109         alpha = x[1];
    110         mx = 0;
    111         for(j=2; j<=n; j++)
     173
     174
     175        /*************************************************************************
     176        Application of an elementary reflection to a rectangular matrix of size MxN
     177
     178        The algorithm pre-multiplies the matrix by an elementary reflection transformation
     179        which is given by column V and scalar Tau (see the description of the
     180        GenerateReflection procedure). Not the whole matrix but only a part of it
     181        is transformed (rows from M1 to M2, columns from N1 to N2). Only the elements
     182        of this submatrix are changed.
     183
     184        Input parameters:
     185            C       -   matrix to be transformed.
     186            Tau     -   scalar defining the transformation.
     187            V       -   column defining the transformation.
     188                        Array whose index ranges within [1..M2-M1+1].
     189            M1, M2  -   range of rows to be transformed.
     190            N1, N2  -   range of columns to be transformed.
     191            WORK    -   working array whose indexes goes from N1 to N2.
     192
     193        Output parameters:
     194            C       -   the result of multiplying the input matrix C by the
     195                        transformation matrix which is given by Tau and V.
     196                        If N1>N2 or M1>M2, C is not modified.
     197
     198          -- LAPACK auxiliary routine (version 3.0) --
     199             Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
     200             Courant Institute, Argonne National Lab, and Rice University
     201             September 30, 1994
     202        *************************************************************************/
     203        public static void applyreflectionfromtheleft(ref double[,] c,
     204            double tau,
     205            ref double[] v,
     206            int m1,
     207            int m2,
     208            int n1,
     209            int n2,
     210            ref double[] work)
    112211        {
    113             mx = Math.Max(Math.Abs(x[j]), mx);
     212            double t = 0;
     213            int i = 0;
     214            int vm = 0;
     215            int i_ = 0;
     216
     217            if( tau==0 | n1>n2 | m1>m2 )
     218            {
     219                return;
     220            }
     221           
     222            //
     223            // w := C' * v
     224            //
     225            vm = m2-m1+1;
     226            for(i=n1; i<=n2; i++)
     227            {
     228                work[i] = 0;
     229            }
     230            for(i=m1; i<=m2; i++)
     231            {
     232                t = v[i+1-m1];
     233                for(i_=n1; i_<=n2;i_++)
     234                {
     235                    work[i_] = work[i_] + t*c[i,i_];
     236                }
     237            }
     238           
     239            //
     240            // C := C - tau * v * w'
     241            //
     242            for(i=m1; i<=m2; i++)
     243            {
     244                t = v[i-m1+1]*tau;
     245                for(i_=n1; i_<=n2;i_++)
     246                {
     247                    c[i,i_] = c[i,i_] - t*work[i_];
     248                }
     249            }
    114250        }
    115         xnorm = 0;
    116         if( mx!=0 )
     251
     252
     253        /*************************************************************************
     254        Application of an elementary reflection to a rectangular matrix of size MxN
     255
     256        The algorithm post-multiplies the matrix by an elementary reflection transformation
     257        which is given by column V and scalar Tau (see the description of the
     258        GenerateReflection procedure). Not the whole matrix but only a part of it
     259        is transformed (rows from M1 to M2, columns from N1 to N2). Only the
     260        elements of this submatrix are changed.
     261
     262        Input parameters:
     263            C       -   matrix to be transformed.
     264            Tau     -   scalar defining the transformation.
     265            V       -   column defining the transformation.
     266                        Array whose index ranges within [1..N2-N1+1].
     267            M1, M2  -   range of rows to be transformed.
     268            N1, N2  -   range of columns to be transformed.
     269            WORK    -   working array whose indexes goes from M1 to M2.
     270
     271        Output parameters:
     272            C       -   the result of multiplying the input matrix C by the
     273                        transformation matrix which is given by Tau and V.
     274                        If N1>N2 or M1>M2, C is not modified.
     275
     276          -- LAPACK auxiliary routine (version 3.0) --
     277             Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
     278             Courant Institute, Argonne National Lab, and Rice University
     279             September 30, 1994
     280        *************************************************************************/
     281        public static void applyreflectionfromtheright(ref double[,] c,
     282            double tau,
     283            ref double[] v,
     284            int m1,
     285            int m2,
     286            int n1,
     287            int n2,
     288            ref double[] work)
    117289        {
    118             for(j=2; j<=n; j++)
    119             {
    120                 xnorm = xnorm+AP.Math.Sqr(x[j]/mx);
    121             }
    122             xnorm = Math.Sqrt(xnorm)*mx;
    123         }
    124         if( xnorm==0 )
    125         {
    126            
    127             //
    128             // H  =  I
    129             //
    130             tau = 0;
    131             return;
    132         }
    133        
    134         //
    135         // general case
    136         //
    137         mx = Math.Max(Math.Abs(alpha), Math.Abs(xnorm));
    138         beta = -(mx*Math.Sqrt(AP.Math.Sqr(alpha/mx)+AP.Math.Sqr(xnorm/mx)));
    139         if( alpha<0 )
    140         {
    141             beta = -beta;
    142         }
    143         tau = (beta-alpha)/beta;
    144         v = 1/(alpha-beta);
    145         for(i_=2; i_<=n;i_++)
    146         {
    147             x[i_] = v*x[i_];
    148         }
    149         x[1] = beta;
    150     }
    151 
    152 
    153     /*************************************************************************
    154     Application of an elementary reflection to a rectangular matrix of size MxN
    155 
    156     The algorithm pre-multiplies the matrix by an elementary reflection transformation
    157     which is given by column V and scalar Tau (see the description of the
    158     GenerateReflection procedure). Not the whole matrix but only a part of it
    159     is transformed (rows from M1 to M2, columns from N1 to N2). Only the elements
    160     of this submatrix are changed.
    161 
    162     Input parameters:
    163         C       -   matrix to be transformed.
    164         Tau     -   scalar defining the transformation.
    165         V       -   column defining the transformation.
    166                     Array whose index ranges within [1..M2-M1+1].
    167         M1, M2  -   range of rows to be transformed.
    168         N1, N2  -   range of columns to be transformed.
    169         WORK    -   working array whose indexes goes from N1 to N2.
    170 
    171     Output parameters:
    172         C       -   the result of multiplying the input matrix C by the
    173                     transformation matrix which is given by Tau and V.
    174                     If N1>N2 or M1>M2, C is not modified.
    175 
    176       -- LAPACK auxiliary routine (version 3.0) --
    177          Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
    178          Courant Institute, Argonne National Lab, and Rice University
    179          September 30, 1994
    180     *************************************************************************/
    181     public static void applyreflectionfromtheleft(ref double[,] c,
    182         double tau,
    183         ref double[] v,
    184         int m1,
    185         int m2,
    186         int n1,
    187         int n2,
    188         ref double[] work)
    189     {
    190         double t = 0;
    191         int i = 0;
    192         int vm = 0;
    193         int i_ = 0;
    194 
    195         if( tau==0 | n1>n2 | m1>m2 )
    196         {
    197             return;
    198         }
    199        
    200         //
    201         // w := C' * v
    202         //
    203         vm = m2-m1+1;
    204         for(i=n1; i<=n2; i++)
    205         {
    206             work[i] = 0;
    207         }
    208         for(i=m1; i<=m2; i++)
    209         {
    210             t = v[i+1-m1];
    211             for(i_=n1; i_<=n2;i_++)
    212             {
    213                 work[i_] = work[i_] + t*c[i,i_];
    214             }
    215         }
    216        
    217         //
    218         // C := C - tau * v * w'
    219         //
    220         for(i=m1; i<=m2; i++)
    221         {
    222             t = v[i-m1+1]*tau;
    223             for(i_=n1; i_<=n2;i_++)
    224             {
    225                 c[i,i_] = c[i,i_] - t*work[i_];
     290            double t = 0;
     291            int i = 0;
     292            int vm = 0;
     293            int i_ = 0;
     294            int i1_ = 0;
     295
     296            if( tau==0 | n1>n2 | m1>m2 )
     297            {
     298                return;
     299            }
     300           
     301            //
     302            // w := C * v
     303            //
     304            vm = n2-n1+1;
     305            for(i=m1; i<=m2; i++)
     306            {
     307                i1_ = (1)-(n1);
     308                t = 0.0;
     309                for(i_=n1; i_<=n2;i_++)
     310                {
     311                    t += c[i,i_]*v[i_+i1_];
     312                }
     313                work[i] = t;
     314            }
     315           
     316            //
     317            // C := C - w * v'
     318            //
     319            for(i=m1; i<=m2; i++)
     320            {
     321                t = work[i]*tau;
     322                i1_ = (1) - (n1);
     323                for(i_=n1; i_<=n2;i_++)
     324                {
     325                    c[i,i_] = c[i,i_] - t*v[i_+i1_];
     326                }
    226327            }
    227328        }
    228329    }
    229 
    230 
    231     /*************************************************************************
    232     Application of an elementary reflection to a rectangular matrix of size MxN
    233 
    234     The algorithm post-multiplies the matrix by an elementary reflection transformation
    235     which is given by column V and scalar Tau (see the description of the
    236     GenerateReflection procedure). Not the whole matrix but only a part of it
    237     is transformed (rows from M1 to M2, columns from N1 to N2). Only the
    238     elements of this submatrix are changed.
    239 
    240     Input parameters:
    241         C       -   matrix to be transformed.
    242         Tau     -   scalar defining the transformation.
    243         V       -   column defining the transformation.
    244                     Array whose index ranges within [1..N2-N1+1].
    245         M1, M2  -   range of rows to be transformed.
    246         N1, N2  -   range of columns to be transformed.
    247         WORK    -   working array whose indexes goes from M1 to M2.
    248 
    249     Output parameters:
    250         C       -   the result of multiplying the input matrix C by the
    251                     transformation matrix which is given by Tau and V.
    252                     If N1>N2 or M1>M2, C is not modified.
    253 
    254       -- LAPACK auxiliary routine (version 3.0) --
    255          Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
    256          Courant Institute, Argonne National Lab, and Rice University
    257          September 30, 1994
    258     *************************************************************************/
    259     public static void applyreflectionfromtheright(ref double[,] c,
    260         double tau,
    261         ref double[] v,
    262         int m1,
    263         int m2,
    264         int n1,
    265         int n2,
    266         ref double[] work)
    267     {
    268         double t = 0;
    269         int i = 0;
    270         int vm = 0;
    271         int i_ = 0;
    272         int i1_ = 0;
    273 
    274         if( tau==0 | n1>n2 | m1>m2 )
    275         {
    276             return;
    277         }
    278        
    279         //
    280         // w := C * v
    281         //
    282         vm = n2-n1+1;
    283         for(i=m1; i<=m2; i++)
    284         {
    285             i1_ = (1)-(n1);
    286             t = 0.0;
    287             for(i_=n1; i_<=n2;i_++)
    288             {
    289                 t += c[i,i_]*v[i_+i1_];
    290             }
    291             work[i] = t;
    292         }
    293        
    294         //
    295         // C := C - w * v'
    296         //
    297         for(i=m1; i<=m2; i++)
    298         {
    299             t = work[i]*tau;
    300             i1_ = (1) - (n1);
    301             for(i_=n1; i_<=n2;i_++)
    302             {
    303                 c[i,i_] = c[i,i_] - t*v[i_+i1_];
    304             }
    305         }
    306     }
    307 
    308 
    309     private static void testreflections()
    310     {
    311         int i = 0;
    312         int j = 0;
    313         int n = 0;
    314         int m = 0;
    315         int maxmn = 0;
    316         double[] x = new double[0];
    317         double[] v = new double[0];
    318         double[] work = new double[0];
    319         double[,] h = new double[0,0];
    320         double[,] a = new double[0,0];
    321         double[,] b = new double[0,0];
    322         double[,] c = new double[0,0];
    323         double tmp = 0;
    324         double beta = 0;
    325         double tau = 0;
    326         double err = 0;
    327         double mer = 0;
    328         double mel = 0;
    329         double meg = 0;
    330         int pass = 0;
    331         int passcount = 0;
    332         int i_ = 0;
    333 
    334         passcount = 1000;
    335         mer = 0;
    336         mel = 0;
    337         meg = 0;
    338         for(pass=1; pass<=passcount; pass++)
    339         {
    340            
    341             //
    342             // Task
    343             //
    344             n = 1+AP.Math.RandomInteger(10);
    345             m = 1+AP.Math.RandomInteger(10);
    346             maxmn = Math.Max(m, n);
    347            
    348             //
    349             // Initialize
    350             //
    351             x = new double[maxmn+1];
    352             v = new double[maxmn+1];
    353             work = new double[maxmn+1];
    354             h = new double[maxmn+1, maxmn+1];
    355             a = new double[maxmn+1, maxmn+1];
    356             b = new double[maxmn+1, maxmn+1];
    357             c = new double[maxmn+1, maxmn+1];
    358            
    359             //
    360             // GenerateReflection
    361             //
    362             for(i=1; i<=n; i++)
    363             {
    364                 x[i] = 2*AP.Math.RandomReal()-1;
    365                 v[i] = x[i];
    366             }
    367             generatereflection(ref v, n, ref tau);
    368             beta = v[1];
    369             v[1] = 1;
    370             for(i=1; i<=n; i++)
    371             {
    372                 for(j=1; j<=n; j++)
    373                 {
    374                     if( i==j )
    375                     {
    376                         h[i,j] = 1-tau*v[i]*v[j];
    377                     }
    378                     else
    379                     {
    380                         h[i,j] = -(tau*v[i]*v[j]);
    381                     }
    382                 }
    383             }
    384             err = 0;
    385             for(i=1; i<=n; i++)
    386             {
    387                 tmp = 0.0;
    388                 for(i_=1; i_<=n;i_++)
    389                 {
    390                     tmp += h[i,i_]*x[i_];
    391                 }
    392                 if( i==1 )
    393                 {
    394                     err = Math.Max(err, Math.Abs(tmp-beta));
    395                 }
    396                 else
    397                 {
    398                     err = Math.Max(err, Math.Abs(tmp));
    399                 }
    400             }
    401             meg = Math.Max(meg, err);
    402            
    403             //
    404             // ApplyReflectionFromTheLeft
    405             //
    406             for(i=1; i<=m; i++)
    407             {
    408                 x[i] = 2*AP.Math.RandomReal()-1;
    409                 v[i] = x[i];
    410             }
    411             for(i=1; i<=m; i++)
    412             {
    413                 for(j=1; j<=n; j++)
    414                 {
    415                     a[i,j] = 2*AP.Math.RandomReal()-1;
    416                     b[i,j] = a[i,j];
    417                 }
    418             }
    419             generatereflection(ref v, m, ref tau);
    420             beta = v[1];
    421             v[1] = 1;
    422             applyreflectionfromtheleft(ref b, tau, ref v, 1, m, 1, n, ref work);
    423             for(i=1; i<=m; i++)
    424             {
    425                 for(j=1; j<=m; j++)
    426                 {
    427                     if( i==j )
    428                     {
    429                         h[i,j] = 1-tau*v[i]*v[j];
    430                     }
    431                     else
    432                     {
    433                         h[i,j] = -(tau*v[i]*v[j]);
    434                     }
    435                 }
    436             }
    437             for(i=1; i<=m; i++)
    438             {
    439                 for(j=1; j<=n; j++)
    440                 {
    441                     tmp = 0.0;
    442                     for(i_=1; i_<=m;i_++)
    443                     {
    444                         tmp += h[i,i_]*a[i_,j];
    445                     }
    446                     c[i,j] = tmp;
    447                 }
    448             }
    449             err = 0;
    450             for(i=1; i<=m; i++)
    451             {
    452                 for(j=1; j<=n; j++)
    453                 {
    454                     err = Math.Max(err, Math.Abs(b[i,j]-c[i,j]));
    455                 }
    456             }
    457             mel = Math.Max(mel, err);
    458            
    459             //
    460             // ApplyReflectionFromTheRight
    461             //
    462             for(i=1; i<=n; i++)
    463             {
    464                 x[i] = 2*AP.Math.RandomReal()-1;
    465                 v[i] = x[i];
    466             }
    467             for(i=1; i<=m; i++)
    468             {
    469                 for(j=1; j<=n; j++)
    470                 {
    471                     a[i,j] = 2*AP.Math.RandomReal()-1;
    472                     b[i,j] = a[i,j];
    473                 }
    474             }
    475             generatereflection(ref v, n, ref tau);
    476             beta = v[1];
    477             v[1] = 1;
    478             applyreflectionfromtheright(ref b, tau, ref v, 1, m, 1, n, ref work);
    479             for(i=1; i<=n; i++)
    480             {
    481                 for(j=1; j<=n; j++)
    482                 {
    483                     if( i==j )
    484                     {
    485                         h[i,j] = 1-tau*v[i]*v[j];
    486                     }
    487                     else
    488                     {
    489                         h[i,j] = -(tau*v[i]*v[j]);
    490                     }
    491                 }
    492             }
    493             for(i=1; i<=m; i++)
    494             {
    495                 for(j=1; j<=n; j++)
    496                 {
    497                     tmp = 0.0;
    498                     for(i_=1; i_<=n;i_++)
    499                     {
    500                         tmp += a[i,i_]*h[i_,j];
    501                     }
    502                     c[i,j] = tmp;
    503                 }
    504             }
    505             err = 0;
    506             for(i=1; i<=m; i++)
    507             {
    508                 for(j=1; j<=n; j++)
    509                 {
    510                     err = Math.Max(err, Math.Abs(b[i,j]-c[i,j]));
    511                 }
    512             }
    513             mer = Math.Max(mer, err);
    514         }
    515        
    516         //
    517         // Overflow crash test
    518         //
    519         x = new double[10+1];
    520         v = new double[10+1];
    521         for(i=1; i<=10; i++)
    522         {
    523             v[i] = AP.Math.MaxRealNumber*0.01*(2*AP.Math.RandomReal()-1);
    524         }
    525         generatereflection(ref v, 10, ref tau);
    526         System.Console.Write("TESTING REFLECTIONS");
    527         System.Console.WriteLine();
    528         System.Console.Write("Pass count is ");
    529         System.Console.Write("{0,0:d}",passcount);
    530         System.Console.WriteLine();
    531         System.Console.Write("Generate     absolute error is       ");
    532         System.Console.Write("{0,5:E3}",meg);
    533         System.Console.WriteLine();
    534         System.Console.Write("Apply(Left)  absolute error is       ");
    535         System.Console.Write("{0,5:E3}",mel);
    536         System.Console.WriteLine();
    537         System.Console.Write("Apply(Right) absolute error is       ");
    538         System.Console.Write("{0,5:E3}",mer);
    539         System.Console.WriteLine();
    540         System.Console.Write("Overflow crash test passed");
    541         System.Console.WriteLine();
    542     }
    543330}
Note: See TracChangeset for help on using the changeset viewer.