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/lq.cs

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