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/bidiagonal.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 bidiagonal
     29namespace alglib
    4230{
    43     /*************************************************************************
    44     Reduction of a rectangular matrix to  bidiagonal form
    45 
    46     The algorithm reduces the rectangular matrix A to  bidiagonal form by
    47     orthogonal transformations P and Q: A = Q*B*P.
    48 
    49     Input parameters:
    50         A       -   source matrix. array[0..M-1, 0..N-1]
    51         M       -   number of rows in matrix A.
    52         N       -   number of columns in matrix A.
    53 
    54     Output parameters:
    55         A       -   matrices Q, B, P in compact form (see below).
    56         TauQ    -   scalar factors which are used to form matrix Q.
    57         TauP    -   scalar factors which are used to form matrix P.
    58 
    59     The main diagonal and one of the  secondary  diagonals  of  matrix  A  are
    60     replaced with bidiagonal  matrix  B.  Other  elements  contain  elementary
    61     reflections which form MxM matrix Q and NxN matrix P, respectively.
    62 
    63     If M>=N, B is the upper  bidiagonal  MxN  matrix  and  is  stored  in  the
    64     corresponding  elements  of  matrix  A.  Matrix  Q  is  represented  as  a
    65     product   of   elementary   reflections   Q = H(0)*H(1)*...*H(n-1),  where
    66     H(i) = 1-tau*v*v'. Here tau is a scalar which is stored  in  TauQ[i],  and
    67     vector v has the following  structure:  v(0:i-1)=0, v(i)=1, v(i+1:m-1)  is
    68     stored   in   elements   A(i+1:m-1,i).   Matrix   P  is  as  follows:  P =
    69     G(0)*G(1)*...*G(n-2), where G(i) = 1 - tau*u*u'. Tau is stored in TauP[i],
    70     u(0:i)=0, u(i+1)=1, u(i+2:n-1) is stored in elements A(i,i+2:n-1).
    71 
    72     If M<N, B is the  lower  bidiagonal  MxN  matrix  and  is  stored  in  the
    73     corresponding   elements  of  matrix  A.  Q = H(0)*H(1)*...*H(m-2),  where
    74     H(i) = 1 - tau*v*v', tau is stored in TauQ, v(0:i)=0, v(i+1)=1, v(i+2:m-1)
    75     is    stored    in   elements   A(i+2:m-1,i).    P = G(0)*G(1)*...*G(m-1),
    76     G(i) = 1-tau*u*u', tau is stored in  TauP,  u(0:i-1)=0, u(i)=1, u(i+1:n-1)
    77     is stored in A(i,i+1:n-1).
    78 
    79     EXAMPLE:
    80 
    81     m=6, n=5 (m > n):               m=5, n=6 (m < n):
    82 
    83     (  d   e   u1  u1  u1 )         (  d   u1  u1  u1  u1  u1 )
    84     (  v1  d   e   u2  u2 )         (  e   d   u2  u2  u2  u2 )
    85     (  v1  v2  d   e   u3 )         (  v1  e   d   u3  u3  u3 )
    86     (  v1  v2  v3  d   e  )         (  v1  v2  e   d   u4  u4 )
    87     (  v1  v2  v3  v4  d  )         (  v1  v2  v3  e   d   u5 )
    88     (  v1  v2  v3  v4  v5 )
    89 
    90     Here vi and ui are vectors which form H(i) and G(i), and d and e -
    91     are the diagonal and off-diagonal elements of matrix B.
    92     *************************************************************************/
    93     public static void rmatrixbd(ref double[,] a,
    94         int m,
    95         int n,
    96         ref double[] tauq,
    97         ref double[] taup)
     31    public class bidiagonal
    9832    {
    99         double[] work = new double[0];
    100         double[] t = new double[0];
    101         int minmn = 0;
    102         int maxmn = 0;
    103         int i = 0;
    104         int j = 0;
    105         double ltau = 0;
    106         int i_ = 0;
    107         int i1_ = 0;
    108 
    109        
    110         //
    111         // Prepare
    112         //
    113         if( n<=0 | m<=0 )
     33        /*************************************************************************
     34        Reduction of a rectangular matrix to  bidiagonal form
     35
     36        The algorithm reduces the rectangular matrix A to  bidiagonal form by
     37        orthogonal transformations P and Q: A = Q*B*P.
     38
     39        Input parameters:
     40            A       -   source matrix. array[0..M-1, 0..N-1]
     41            M       -   number of rows in matrix A.
     42            N       -   number of columns in matrix A.
     43
     44        Output parameters:
     45            A       -   matrices Q, B, P in compact form (see below).
     46            TauQ    -   scalar factors which are used to form matrix Q.
     47            TauP    -   scalar factors which are used to form matrix P.
     48
     49        The main diagonal and one of the  secondary  diagonals  of  matrix  A  are
     50        replaced with bidiagonal  matrix  B.  Other  elements  contain  elementary
     51        reflections which form MxM matrix Q and NxN matrix P, respectively.
     52
     53        If M>=N, B is the upper  bidiagonal  MxN  matrix  and  is  stored  in  the
     54        corresponding  elements  of  matrix  A.  Matrix  Q  is  represented  as  a
     55        product   of   elementary   reflections   Q = H(0)*H(1)*...*H(n-1),  where
     56        H(i) = 1-tau*v*v'. Here tau is a scalar which is stored  in  TauQ[i],  and
     57        vector v has the following  structure:  v(0:i-1)=0, v(i)=1, v(i+1:m-1)  is
     58        stored   in   elements   A(i+1:m-1,i).   Matrix   P  is  as  follows:  P =
     59        G(0)*G(1)*...*G(n-2), where G(i) = 1 - tau*u*u'. Tau is stored in TauP[i],
     60        u(0:i)=0, u(i+1)=1, u(i+2:n-1) is stored in elements A(i,i+2:n-1).
     61
     62        If M<N, B is the  lower  bidiagonal  MxN  matrix  and  is  stored  in  the
     63        corresponding   elements  of  matrix  A.  Q = H(0)*H(1)*...*H(m-2),  where
     64        H(i) = 1 - tau*v*v', tau is stored in TauQ, v(0:i)=0, v(i+1)=1, v(i+2:m-1)
     65        is    stored    in   elements   A(i+2:m-1,i).    P = G(0)*G(1)*...*G(m-1),
     66        G(i) = 1-tau*u*u', tau is stored in  TauP,  u(0:i-1)=0, u(i)=1, u(i+1:n-1)
     67        is stored in A(i,i+1:n-1).
     68
     69        EXAMPLE:
     70
     71        m=6, n=5 (m > n):               m=5, n=6 (m < n):
     72
     73        (  d   e   u1  u1  u1 )         (  d   u1  u1  u1  u1  u1 )
     74        (  v1  d   e   u2  u2 )         (  e   d   u2  u2  u2  u2 )
     75        (  v1  v2  d   e   u3 )         (  v1  e   d   u3  u3  u3 )
     76        (  v1  v2  v3  d   e  )         (  v1  v2  e   d   u4  u4 )
     77        (  v1  v2  v3  v4  d  )         (  v1  v2  v3  e   d   u5 )
     78        (  v1  v2  v3  v4  v5 )
     79
     80        Here vi and ui are vectors which form H(i) and G(i), and d and e -
     81        are the diagonal and off-diagonal elements of matrix B.
     82        *************************************************************************/
     83        public static void rmatrixbd(ref double[,] a,
     84            int m,
     85            int n,
     86            ref double[] tauq,
     87            ref double[] taup)
    11488        {
    115             return;
    116         }
    117         minmn = Math.Min(m, n);
    118         maxmn = Math.Max(m, n);
    119         work = new double[maxmn+1];
    120         t = new double[maxmn+1];
    121         if( m>=n )
    122         {
    123             tauq = new double[n-1+1];
    124             taup = new double[n-1+1];
    125         }
    126         else
    127         {
    128             tauq = new double[m-1+1];
    129             taup = new double[m-1+1];
    130         }
    131         if( m>=n )
    132         {
     89            double[] work = new double[0];
     90            double[] t = new double[0];
     91            int minmn = 0;
     92            int maxmn = 0;
     93            int i = 0;
     94            int j = 0;
     95            double ltau = 0;
     96            int i_ = 0;
     97            int i1_ = 0;
     98
    13399           
    134100            //
    135             // Reduce to upper bidiagonal form
    136             //
    137             for(i=0; i<=n-1; i++)
    138             {
    139                
    140                 //
    141                 // Generate elementary reflector H(i) to annihilate A(i+1:m-1,i)
    142                 //
    143                 i1_ = (i) - (1);
    144                 for(i_=1; i_<=m-i;i_++)
    145                 {
    146                     t[i_] = a[i_+i1_,i];
    147                 }
    148                 reflections.generatereflection(ref t, m-i, ref ltau);
    149                 tauq[i] = ltau;
    150                 i1_ = (1) - (i);
    151                 for(i_=i; i_<=m-1;i_++)
    152                 {
    153                     a[i_,i] = t[i_+i1_];
    154                 }
    155                 t[1] = 1;
    156                
    157                 //
    158                 // Apply H(i) to A(i:m-1,i+1:n-1) from the left
    159                 //
    160                 reflections.applyreflectionfromtheleft(ref a, ltau, ref t, i, m-1, i+1, n-1, ref work);
    161                 if( i<n-1 )
     101            // Prepare
     102            //
     103            if( n<=0 | m<=0 )
     104            {
     105                return;
     106            }
     107            minmn = Math.Min(m, n);
     108            maxmn = Math.Max(m, n);
     109            work = new double[maxmn+1];
     110            t = new double[maxmn+1];
     111            if( m>=n )
     112            {
     113                tauq = new double[n-1+1];
     114                taup = new double[n-1+1];
     115            }
     116            else
     117            {
     118                tauq = new double[m-1+1];
     119                taup = new double[m-1+1];
     120            }
     121            if( m>=n )
     122            {
     123               
     124                //
     125                // Reduce to upper bidiagonal form
     126                //
     127                for(i=0; i<=n-1; i++)
    162128                {
    163129                   
    164130                    //
    165                     // Generate elementary reflector G(i) to annihilate
    166                     // A(i,i+2:n-1)
    167                     //
    168                     i1_ = (i+1) - (1);
    169                     for(i_=1; i_<=n-i-1;i_++)
    170                     {
    171                         t[i_] = a[i,i_+i1_];
    172                     }
    173                     reflections.generatereflection(ref t, n-1-i, ref ltau);
    174                     taup[i] = ltau;
    175                     i1_ = (1) - (i+1);
    176                     for(i_=i+1; i_<=n-1;i_++)
    177                     {
    178                         a[i,i_] = t[i_+i1_];
     131                    // Generate elementary reflector H(i) to annihilate A(i+1:m-1,i)
     132                    //
     133                    i1_ = (i) - (1);
     134                    for(i_=1; i_<=m-i;i_++)
     135                    {
     136                        t[i_] = a[i_+i1_,i];
     137                    }
     138                    reflections.generatereflection(ref t, m-i, ref ltau);
     139                    tauq[i] = ltau;
     140                    i1_ = (1) - (i);
     141                    for(i_=i; i_<=m-1;i_++)
     142                    {
     143                        a[i_,i] = t[i_+i1_];
    179144                    }
    180145                    t[1] = 1;
    181146                   
    182147                    //
    183                     // Apply G(i) to A(i+1:m-1,i+1:n-1) from the right
    184                     //
    185                     reflections.applyreflectionfromtheright(ref a, ltau, ref t, i+1, m-1, i+1, n-1, ref work);
    186                 }
    187                 else
    188                 {
    189                     taup[i] = 0;
    190                 }
    191             }
    192         }
    193         else
    194         {
    195            
    196             //
    197             // Reduce to lower bidiagonal form
    198             //
    199             for(i=0; i<=m-1; i++)
    200             {
    201                
    202                 //
    203                 // Generate elementary reflector G(i) to annihilate A(i,i+1:n-1)
    204                 //
    205                 i1_ = (i) - (1);
    206                 for(i_=1; i_<=n-i;i_++)
    207                 {
    208                     t[i_] = a[i,i_+i1_];
    209                 }
    210                 reflections.generatereflection(ref t, n-i, ref ltau);
    211                 taup[i] = ltau;
    212                 i1_ = (1) - (i);
    213                 for(i_=i; i_<=n-1;i_++)
    214                 {
    215                     a[i,i_] = t[i_+i1_];
    216                 }
    217                 t[1] = 1;
    218                
    219                 //
    220                 // Apply G(i) to A(i+1:m-1,i:n-1) from the right
    221                 //
    222                 reflections.applyreflectionfromtheright(ref a, ltau, ref t, i+1, m-1, i, n-1, ref work);
    223                 if( i<m-1 )
     148                    // Apply H(i) to A(i:m-1,i+1:n-1) from the left
     149                    //
     150                    reflections.applyreflectionfromtheleft(ref a, ltau, ref t, i, m-1, i+1, n-1, ref work);
     151                    if( i<n-1 )
     152                    {
     153                       
     154                        //
     155                        // Generate elementary reflector G(i) to annihilate
     156                        // A(i,i+2:n-1)
     157                        //
     158                        i1_ = (i+1) - (1);
     159                        for(i_=1; i_<=n-i-1;i_++)
     160                        {
     161                            t[i_] = a[i,i_+i1_];
     162                        }
     163                        reflections.generatereflection(ref t, n-1-i, ref ltau);
     164                        taup[i] = ltau;
     165                        i1_ = (1) - (i+1);
     166                        for(i_=i+1; i_<=n-1;i_++)
     167                        {
     168                            a[i,i_] = t[i_+i1_];
     169                        }
     170                        t[1] = 1;
     171                       
     172                        //
     173                        // Apply G(i) to A(i+1:m-1,i+1:n-1) from the right
     174                        //
     175                        reflections.applyreflectionfromtheright(ref a, ltau, ref t, i+1, m-1, i+1, n-1, ref work);
     176                    }
     177                    else
     178                    {
     179                        taup[i] = 0;
     180                    }
     181                }
     182            }
     183            else
     184            {
     185               
     186                //
     187                // Reduce to lower bidiagonal form
     188                //
     189                for(i=0; i<=m-1; i++)
    224190                {
    225191                   
    226192                    //
    227                     // Generate elementary reflector H(i) to annihilate
    228                     // A(i+2:m-1,i)
    229                     //
    230                     i1_ = (i+1) - (1);
    231                     for(i_=1; i_<=m-1-i;i_++)
    232                     {
    233                         t[i_] = a[i_+i1_,i];
    234                     }
    235                     reflections.generatereflection(ref t, m-1-i, ref ltau);
    236                     tauq[i] = ltau;
    237                     i1_ = (1) - (i+1);
    238                     for(i_=i+1; i_<=m-1;i_++)
    239                     {
    240                         a[i_,i] = t[i_+i1_];
     193                    // Generate elementary reflector G(i) to annihilate A(i,i+1:n-1)
     194                    //
     195                    i1_ = (i) - (1);
     196                    for(i_=1; i_<=n-i;i_++)
     197                    {
     198                        t[i_] = a[i,i_+i1_];
     199                    }
     200                    reflections.generatereflection(ref t, n-i, ref ltau);
     201                    taup[i] = ltau;
     202                    i1_ = (1) - (i);
     203                    for(i_=i; i_<=n-1;i_++)
     204                    {
     205                        a[i,i_] = t[i_+i1_];
    241206                    }
    242207                    t[1] = 1;
    243208                   
    244209                    //
    245                     // Apply H(i) to A(i+1:m-1,i+1:n-1) from the left
    246                     //
    247                     reflections.applyreflectionfromtheleft(ref a, ltau, ref t, i+1, m-1, i+1, n-1, ref work);
     210                    // Apply G(i) to A(i+1:m-1,i:n-1) from the right
     211                    //
     212                    reflections.applyreflectionfromtheright(ref a, ltau, ref t, i+1, m-1, i, n-1, ref work);
     213                    if( i<m-1 )
     214                    {
     215                       
     216                        //
     217                        // Generate elementary reflector H(i) to annihilate
     218                        // A(i+2:m-1,i)
     219                        //
     220                        i1_ = (i+1) - (1);
     221                        for(i_=1; i_<=m-1-i;i_++)
     222                        {
     223                            t[i_] = a[i_+i1_,i];
     224                        }
     225                        reflections.generatereflection(ref t, m-1-i, ref ltau);
     226                        tauq[i] = ltau;
     227                        i1_ = (1) - (i+1);
     228                        for(i_=i+1; i_<=m-1;i_++)
     229                        {
     230                            a[i_,i] = t[i_+i1_];
     231                        }
     232                        t[1] = 1;
     233                       
     234                        //
     235                        // Apply H(i) to A(i+1:m-1,i+1:n-1) from the left
     236                        //
     237                        reflections.applyreflectionfromtheleft(ref a, ltau, ref t, i+1, m-1, i+1, n-1, ref work);
     238                    }
     239                    else
     240                    {
     241                        tauq[i] = 0;
     242                    }
     243                }
     244            }
     245        }
     246
     247
     248        /*************************************************************************
     249        Unpacking matrix Q which reduces a matrix to bidiagonal form.
     250
     251        Input parameters:
     252            QP          -   matrices Q and P in compact form.
     253                            Output of ToBidiagonal subroutine.
     254            M           -   number of rows in matrix A.
     255            N           -   number of columns in matrix A.
     256            TAUQ        -   scalar factors which are used to form Q.
     257                            Output of ToBidiagonal subroutine.
     258            QColumns    -   required number of columns in matrix Q.
     259                            M>=QColumns>=0.
     260
     261        Output parameters:
     262            Q           -   first QColumns columns of matrix Q.
     263                            Array[0..M-1, 0..QColumns-1]
     264                            If QColumns=0, the array is not modified.
     265
     266          -- ALGLIB --
     267             Copyright 2005 by Bochkanov Sergey
     268        *************************************************************************/
     269        public static void rmatrixbdunpackq(ref double[,] qp,
     270            int m,
     271            int n,
     272            ref double[] tauq,
     273            int qcolumns,
     274            ref double[,] q)
     275        {
     276            int i = 0;
     277            int j = 0;
     278
     279            System.Diagnostics.Debug.Assert(qcolumns<=m, "RMatrixBDUnpackQ: QColumns>M!");
     280            System.Diagnostics.Debug.Assert(qcolumns>=0, "RMatrixBDUnpackQ: QColumns<0!");
     281            if( m==0 | n==0 | qcolumns==0 )
     282            {
     283                return;
     284            }
     285           
     286            //
     287            // prepare Q
     288            //
     289            q = new double[m-1+1, qcolumns-1+1];
     290            for(i=0; i<=m-1; i++)
     291            {
     292                for(j=0; j<=qcolumns-1; j++)
     293                {
     294                    if( i==j )
     295                    {
     296                        q[i,j] = 1;
     297                    }
     298                    else
     299                    {
     300                        q[i,j] = 0;
     301                    }
     302                }
     303            }
     304           
     305            //
     306            // Calculate
     307            //
     308            rmatrixbdmultiplybyq(ref qp, m, n, ref tauq, ref q, m, qcolumns, false, false);
     309        }
     310
     311
     312        /*************************************************************************
     313        Multiplication by matrix Q which reduces matrix A to  bidiagonal form.
     314
     315        The algorithm allows pre- or post-multiply by Q or Q'.
     316
     317        Input parameters:
     318            QP          -   matrices Q and P in compact form.
     319                            Output of ToBidiagonal subroutine.
     320            M           -   number of rows in matrix A.
     321            N           -   number of columns in matrix A.
     322            TAUQ        -   scalar factors which are used to form Q.
     323                            Output of ToBidiagonal subroutine.
     324            Z           -   multiplied matrix.
     325                            array[0..ZRows-1,0..ZColumns-1]
     326            ZRows       -   number of rows in matrix Z. If FromTheRight=False,
     327                            ZRows=M, otherwise ZRows can be arbitrary.
     328            ZColumns    -   number of columns in matrix Z. If FromTheRight=True,
     329                            ZColumns=M, otherwise ZColumns can be arbitrary.
     330            FromTheRight -  pre- or post-multiply.
     331            DoTranspose -   multiply by Q or Q'.
     332
     333        Output parameters:
     334            Z           -   product of Z and Q.
     335                            Array[0..ZRows-1,0..ZColumns-1]
     336                            If ZRows=0 or ZColumns=0, the array is not modified.
     337
     338          -- ALGLIB --
     339             Copyright 2005 by Bochkanov Sergey
     340        *************************************************************************/
     341        public static void rmatrixbdmultiplybyq(ref double[,] qp,
     342            int m,
     343            int n,
     344            ref double[] tauq,
     345            ref double[,] z,
     346            int zrows,
     347            int zcolumns,
     348            bool fromtheright,
     349            bool dotranspose)
     350        {
     351            int i = 0;
     352            int i1 = 0;
     353            int i2 = 0;
     354            int istep = 0;
     355            double[] v = new double[0];
     356            double[] work = new double[0];
     357            int mx = 0;
     358            int i_ = 0;
     359            int i1_ = 0;
     360
     361            if( m<=0 | n<=0 | zrows<=0 | zcolumns<=0 )
     362            {
     363                return;
     364            }
     365            System.Diagnostics.Debug.Assert(fromtheright & zcolumns==m | !fromtheright & zrows==m, "RMatrixBDMultiplyByQ: incorrect Z size!");
     366           
     367            //
     368            // init
     369            //
     370            mx = Math.Max(m, n);
     371            mx = Math.Max(mx, zrows);
     372            mx = Math.Max(mx, zcolumns);
     373            v = new double[mx+1];
     374            work = new double[mx+1];
     375            if( m>=n )
     376            {
     377               
     378                //
     379                // setup
     380                //
     381                if( fromtheright )
     382                {
     383                    i1 = 0;
     384                    i2 = n-1;
     385                    istep = +1;
    248386                }
    249387                else
    250388                {
    251                     tauq[i] = 0;
    252                 }
    253             }
    254         }
    255     }
    256 
    257 
    258     /*************************************************************************
    259     Unpacking matrix Q which reduces a matrix to bidiagonal form.
    260 
    261     Input parameters:
    262         QP          -   matrices Q and P in compact form.
    263                         Output of ToBidiagonal subroutine.
    264         M           -   number of rows in matrix A.
    265         N           -   number of columns in matrix A.
    266         TAUQ        -   scalar factors which are used to form Q.
    267                         Output of ToBidiagonal subroutine.
    268         QColumns    -   required number of columns in matrix Q.
    269                         M>=QColumns>=0.
    270 
    271     Output parameters:
    272         Q           -   first QColumns columns of matrix Q.
    273                         Array[0..M-1, 0..QColumns-1]
    274                         If QColumns=0, the array is not modified.
    275 
    276       -- ALGLIB --
    277          Copyright 2005 by Bochkanov Sergey
    278     *************************************************************************/
    279     public static void rmatrixbdunpackq(ref double[,] qp,
    280         int m,
    281         int n,
    282         ref double[] tauq,
    283         int qcolumns,
    284         ref double[,] q)
    285     {
    286         int i = 0;
    287         int j = 0;
    288 
    289         System.Diagnostics.Debug.Assert(qcolumns<=m, "RMatrixBDUnpackQ: QColumns>M!");
    290         System.Diagnostics.Debug.Assert(qcolumns>=0, "RMatrixBDUnpackQ: QColumns<0!");
    291         if( m==0 | n==0 | qcolumns==0 )
    292         {
    293             return;
    294         }
    295        
    296         //
    297         // prepare Q
    298         //
    299         q = new double[m-1+1, qcolumns-1+1];
    300         for(i=0; i<=m-1; i++)
    301         {
    302             for(j=0; j<=qcolumns-1; j++)
    303             {
    304                 if( i==j )
    305                 {
    306                     q[i,j] = 1;
    307                 }
    308                 else
    309                 {
    310                     q[i,j] = 0;
    311                 }
    312             }
    313         }
    314        
    315         //
    316         // Calculate
    317         //
    318         rmatrixbdmultiplybyq(ref qp, m, n, ref tauq, ref q, m, qcolumns, false, false);
    319     }
    320 
    321 
    322     /*************************************************************************
    323     Multiplication by matrix Q which reduces matrix A to  bidiagonal form.
    324 
    325     The algorithm allows pre- or post-multiply by Q or Q'.
    326 
    327     Input parameters:
    328         QP          -   matrices Q and P in compact form.
    329                         Output of ToBidiagonal subroutine.
    330         M           -   number of rows in matrix A.
    331         N           -   number of columns in matrix A.
    332         TAUQ        -   scalar factors which are used to form Q.
    333                         Output of ToBidiagonal subroutine.
    334         Z           -   multiplied matrix.
    335                         array[0..ZRows-1,0..ZColumns-1]
    336         ZRows       -   number of rows in matrix Z. If FromTheRight=False,
    337                         ZRows=M, otherwise ZRows can be arbitrary.
    338         ZColumns    -   number of columns in matrix Z. If FromTheRight=True,
    339                         ZColumns=M, otherwise ZColumns can be arbitrary.
    340         FromTheRight -  pre- or post-multiply.
    341         DoTranspose -   multiply by Q or Q'.
    342 
    343     Output parameters:
    344         Z           -   product of Z and Q.
    345                         Array[0..ZRows-1,0..ZColumns-1]
    346                         If ZRows=0 or ZColumns=0, the array is not modified.
    347 
    348       -- ALGLIB --
    349          Copyright 2005 by Bochkanov Sergey
    350     *************************************************************************/
    351     public static void rmatrixbdmultiplybyq(ref double[,] qp,
    352         int m,
    353         int n,
    354         ref double[] tauq,
    355         ref double[,] z,
    356         int zrows,
    357         int zcolumns,
    358         bool fromtheright,
    359         bool dotranspose)
    360     {
    361         int i = 0;
    362         int i1 = 0;
    363         int i2 = 0;
    364         int istep = 0;
    365         double[] v = new double[0];
    366         double[] work = new double[0];
    367         int mx = 0;
    368         int i_ = 0;
    369         int i1_ = 0;
    370 
    371         if( m<=0 | n<=0 | zrows<=0 | zcolumns<=0 )
    372         {
    373             return;
    374         }
    375         System.Diagnostics.Debug.Assert(fromtheright & zcolumns==m | !fromtheright & zrows==m, "RMatrixBDMultiplyByQ: incorrect Z size!");
    376        
    377         //
    378         // init
    379         //
    380         mx = Math.Max(m, n);
    381         mx = Math.Max(mx, zrows);
    382         mx = Math.Max(mx, zcolumns);
    383         v = new double[mx+1];
    384         work = new double[mx+1];
    385         if( m>=n )
    386         {
    387            
    388             //
    389             // setup
    390             //
    391             if( fromtheright )
    392             {
    393                 i1 = 0;
    394                 i2 = n-1;
    395                 istep = +1;
    396             }
    397             else
    398             {
    399                 i1 = n-1;
    400                 i2 = 0;
    401                 istep = -1;
    402             }
    403             if( dotranspose )
    404             {
    405                 i = i1;
    406                 i1 = i2;
    407                 i2 = i;
    408                 istep = -istep;
    409             }
    410            
    411             //
    412             // Process
    413             //
    414             i = i1;
    415             do
    416             {
    417                 i1_ = (i) - (1);
    418                 for(i_=1; i_<=m-i;i_++)
    419                 {
    420                     v[i_] = qp[i_+i1_,i];
    421                 }
    422                 v[1] = 1;
    423                 if( fromtheright )
    424                 {
    425                     reflections.applyreflectionfromtheright(ref z, tauq[i], ref v, 0, zrows-1, i, m-1, ref work);
    426                 }
    427                 else
    428                 {
    429                     reflections.applyreflectionfromtheleft(ref z, tauq[i], ref v, i, m-1, 0, zcolumns-1, ref work);
    430                 }
    431                 i = i+istep;
    432             }
    433             while( i!=i2+istep );
    434         }
    435         else
    436         {
    437            
    438             //
    439             // setup
    440             //
    441             if( fromtheright )
    442             {
    443                 i1 = 0;
    444                 i2 = m-2;
    445                 istep = +1;
    446             }
    447             else
    448             {
    449                 i1 = m-2;
    450                 i2 = 0;
    451                 istep = -1;
    452             }
    453             if( dotranspose )
    454             {
    455                 i = i1;
    456                 i1 = i2;
    457                 i2 = i;
    458                 istep = -istep;
    459             }
    460            
    461             //
    462             // Process
    463             //
    464             if( m-1>0 )
    465             {
     389                    i1 = n-1;
     390                    i2 = 0;
     391                    istep = -1;
     392                }
     393                if( dotranspose )
     394                {
     395                    i = i1;
     396                    i1 = i2;
     397                    i2 = i;
     398                    istep = -istep;
     399                }
     400               
     401                //
     402                // Process
     403                //
    466404                i = i1;
    467405                do
    468406                {
    469                     i1_ = (i+1) - (1);
    470                     for(i_=1; i_<=m-i-1;i_++)
     407                    i1_ = (i) - (1);
     408                    for(i_=1; i_<=m-i;i_++)
    471409                    {
    472410                        v[i_] = qp[i_+i1_,i];
    … …  
    475413                    if( fromtheright )
    476414                    {
    477                         reflections.applyreflectionfromtheright(ref z, tauq[i], ref v, 0, zrows-1, i+1, m-1, ref work);
     415                        reflections.applyreflectionfromtheright(ref z, tauq[i], ref v, 0, zrows-1, i, m-1, ref work);
    478416                    }
    479417                    else
    480418                    {
    481                         reflections.applyreflectionfromtheleft(ref z, tauq[i], ref v, i+1, m-1, 0, zcolumns-1, ref work);
     419                        reflections.applyreflectionfromtheleft(ref z, tauq[i], ref v, i, m-1, 0, zcolumns-1, ref work);
    482420                    }
    483421                    i = i+istep;
    … …  
    485423                while( i!=i2+istep );
    486424            }
     425            else
     426            {
     427               
     428                //
     429                // setup
     430                //
     431                if( fromtheright )
     432                {
     433                    i1 = 0;
     434                    i2 = m-2;
     435                    istep = +1;
     436                }
     437                else
     438                {
     439                    i1 = m-2;
     440                    i2 = 0;
     441                    istep = -1;
     442                }
     443                if( dotranspose )
     444                {
     445                    i = i1;
     446                    i1 = i2;
     447                    i2 = i;
     448                    istep = -istep;
     449                }
     450               
     451                //
     452                // Process
     453                //
     454                if( m-1>0 )
     455                {
     456                    i = i1;
     457                    do
     458                    {
     459                        i1_ = (i+1) - (1);
     460                        for(i_=1; i_<=m-i-1;i_++)
     461                        {
     462                            v[i_] = qp[i_+i1_,i];
     463                        }
     464                        v[1] = 1;
     465                        if( fromtheright )
     466                        {
     467                            reflections.applyreflectionfromtheright(ref z, tauq[i], ref v, 0, zrows-1, i+1, m-1, ref work);
     468                        }
     469                        else
     470                        {
     471                            reflections.applyreflectionfromtheleft(ref z, tauq[i], ref v, i+1, m-1, 0, zcolumns-1, ref work);
     472                        }
     473                        i = i+istep;
     474                    }
     475                    while( i!=i2+istep );
     476                }
     477            }
    487478        }
    488     }
    489 
    490 
    491     /*************************************************************************
    492     Unpacking matrix P which reduces matrix A to bidiagonal form.
    493     The subroutine returns transposed matrix P.
    494 
    495     Input parameters:
    496         QP      -   matrices Q and P in compact form.
    497                     Output of ToBidiagonal subroutine.
    498         M       -   number of rows in matrix A.
    499         N       -   number of columns in matrix A.
    500         TAUP    -   scalar factors which are used to form P.
    501                     Output of ToBidiagonal subroutine.
    502         PTRows  -   required number of rows of matrix P^T. N >= PTRows >= 0.
    503 
    504     Output parameters:
    505         PT      -   first PTRows columns of matrix P^T
    506                     Array[0..PTRows-1, 0..N-1]
    507                     If PTRows=0, the array is not modified.
    508 
    509       -- ALGLIB --
    510          Copyright 2005-2007 by Bochkanov Sergey
    511     *************************************************************************/
    512     public static void rmatrixbdunpackpt(ref double[,] qp,
    513         int m,
    514         int n,
    515         ref double[] taup,
    516         int ptrows,
    517         ref double[,] pt)
    518     {
    519         int i = 0;
    520         int j = 0;
    521 
    522         System.Diagnostics.Debug.Assert(ptrows<=n, "RMatrixBDUnpackPT: PTRows>N!");
    523         System.Diagnostics.Debug.Assert(ptrows>=0, "RMatrixBDUnpackPT: PTRows<0!");
    524         if( m==0 | n==0 | ptrows==0 )
     479
     480
     481        /*************************************************************************
     482        Unpacking matrix P which reduces matrix A to bidiagonal form.
     483        The subroutine returns transposed matrix P.
     484
     485        Input parameters:
     486            QP      -   matrices Q and P in compact form.
     487                        Output of ToBidiagonal subroutine.
     488            M       -   number of rows in matrix A.
     489            N       -   number of columns in matrix A.
     490            TAUP    -   scalar factors which are used to form P.
     491                        Output of ToBidiagonal subroutine.
     492            PTRows  -   required number of rows of matrix P^T. N >= PTRows >= 0.
     493
     494        Output parameters:
     495            PT      -   first PTRows columns of matrix P^T
     496                        Array[0..PTRows-1, 0..N-1]
     497                        If PTRows=0, the array is not modified.
     498
     499          -- ALGLIB --
     500             Copyright 2005-2007 by Bochkanov Sergey
     501        *************************************************************************/
     502        public static void rmatrixbdunpackpt(ref double[,] qp,
     503            int m,
     504            int n,
     505            ref double[] taup,
     506            int ptrows,
     507            ref double[,] pt)
    525508        {
    526             return;
     509            int i = 0;
     510            int j = 0;
     511
     512            System.Diagnostics.Debug.Assert(ptrows<=n, "RMatrixBDUnpackPT: PTRows>N!");
     513            System.Diagnostics.Debug.Assert(ptrows>=0, "RMatrixBDUnpackPT: PTRows<0!");
     514            if( m==0 | n==0 | ptrows==0 )
     515            {
     516                return;
     517            }
     518           
     519            //
     520            // prepare PT
     521            //
     522            pt = new double[ptrows-1+1, n-1+1];
     523            for(i=0; i<=ptrows-1; i++)
     524            {
     525                for(j=0; j<=n-1; j++)
     526                {
     527                    if( i==j )
     528                    {
     529                        pt[i,j] = 1;
     530                    }
     531                    else
     532                    {
     533                        pt[i,j] = 0;
     534                    }
     535                }
     536            }
     537           
     538            //
     539            // Calculate
     540            //
     541            rmatrixbdmultiplybyp(ref qp, m, n, ref taup, ref pt, ptrows, n, true, true);
    527542        }
    528        
    529         //
    530         // prepare PT
    531         //
    532         pt = new double[ptrows-1+1, n-1+1];
    533         for(i=0; i<=ptrows-1; i++)
     543
     544
     545        /*************************************************************************
     546        Multiplication by matrix P which reduces matrix A to  bidiagonal form.
     547
     548        The algorithm allows pre- or post-multiply by P or P'.
     549
     550        Input parameters:
     551            QP          -   matrices Q and P in compact form.
     552                            Output of RMatrixBD subroutine.
     553            M           -   number of rows in matrix A.
     554            N           -   number of columns in matrix A.
     555            TAUP        -   scalar factors which are used to form P.
     556                            Output of RMatrixBD subroutine.
     557            Z           -   multiplied matrix.
     558                            Array whose indexes range within [0..ZRows-1,0..ZColumns-1].
     559            ZRows       -   number of rows in matrix Z. If FromTheRight=False,
     560                            ZRows=N, otherwise ZRows can be arbitrary.
     561            ZColumns    -   number of columns in matrix Z. If FromTheRight=True,
     562                            ZColumns=N, otherwise ZColumns can be arbitrary.
     563            FromTheRight -  pre- or post-multiply.
     564            DoTranspose -   multiply by P or P'.
     565
     566        Output parameters:
     567            Z - product of Z and P.
     568                        Array whose indexes range within [0..ZRows-1,0..ZColumns-1].
     569                        If ZRows=0 or ZColumns=0, the array is not modified.
     570
     571          -- ALGLIB --
     572             Copyright 2005-2007 by Bochkanov Sergey
     573        *************************************************************************/
     574        public static void rmatrixbdmultiplybyp(ref double[,] qp,
     575            int m,
     576            int n,
     577            ref double[] taup,
     578            ref double[,] z,
     579            int zrows,
     580            int zcolumns,
     581            bool fromtheright,
     582            bool dotranspose)
    534583        {
    535             for(j=0; j<=n-1; j++)
    536             {
    537                 if( i==j )
    538                 {
    539                     pt[i,j] = 1;
     584            int i = 0;
     585            double[] v = new double[0];
     586            double[] work = new double[0];
     587            int mx = 0;
     588            int i1 = 0;
     589            int i2 = 0;
     590            int istep = 0;
     591            int i_ = 0;
     592            int i1_ = 0;
     593
     594            if( m<=0 | n<=0 | zrows<=0 | zcolumns<=0 )
     595            {
     596                return;
     597            }
     598            System.Diagnostics.Debug.Assert(fromtheright & zcolumns==n | !fromtheright & zrows==n, "RMatrixBDMultiplyByP: incorrect Z size!");
     599           
     600            //
     601            // init
     602            //
     603            mx = Math.Max(m, n);
     604            mx = Math.Max(mx, zrows);
     605            mx = Math.Max(mx, zcolumns);
     606            v = new double[mx+1];
     607            work = new double[mx+1];
     608            v = new double[mx+1];
     609            work = new double[mx+1];
     610            if( m>=n )
     611            {
     612               
     613                //
     614                // setup
     615                //
     616                if( fromtheright )
     617                {
     618                    i1 = n-2;
     619                    i2 = 0;
     620                    istep = -1;
    540621                }
    541622                else
    542623                {
    543                     pt[i,j] = 0;
    544                 }
    545             }
    546         }
    547        
    548         //
    549         // Calculate
    550         //
    551         rmatrixbdmultiplybyp(ref qp, m, n, ref taup, ref pt, ptrows, n, true, true);
    552     }
    553 
    554 
    555     /*************************************************************************
    556     Multiplication by matrix P which reduces matrix A to  bidiagonal form.
    557 
    558     The algorithm allows pre- or post-multiply by P or P'.
    559 
    560     Input parameters:
    561         QP          -   matrices Q and P in compact form.
    562                         Output of RMatrixBD subroutine.
    563         M           -   number of rows in matrix A.
    564         N           -   number of columns in matrix A.
    565         TAUP        -   scalar factors which are used to form P.
    566                         Output of RMatrixBD subroutine.
    567         Z           -   multiplied matrix.
    568                         Array whose indexes range within [0..ZRows-1,0..ZColumns-1].
    569         ZRows       -   number of rows in matrix Z. If FromTheRight=False,
    570                         ZRows=N, otherwise ZRows can be arbitrary.
    571         ZColumns    -   number of columns in matrix Z. If FromTheRight=True,
    572                         ZColumns=N, otherwise ZColumns can be arbitrary.
    573         FromTheRight -  pre- or post-multiply.
    574         DoTranspose -   multiply by P or P'.
    575 
    576     Output parameters:
    577         Z - product of Z and P.
    578                     Array whose indexes range within [0..ZRows-1,0..ZColumns-1].
    579                     If ZRows=0 or ZColumns=0, the array is not modified.
    580 
    581       -- ALGLIB --
    582          Copyright 2005-2007 by Bochkanov Sergey
    583     *************************************************************************/
    584     public static void rmatrixbdmultiplybyp(ref double[,] qp,
    585         int m,
    586         int n,
    587         ref double[] taup,
    588         ref double[,] z,
    589         int zrows,
    590         int zcolumns,
    591         bool fromtheright,
    592         bool dotranspose)
    593     {
    594         int i = 0;
    595         double[] v = new double[0];
    596         double[] work = new double[0];
    597         int mx = 0;
    598         int i1 = 0;
    599         int i2 = 0;
    600         int istep = 0;
    601         int i_ = 0;
    602         int i1_ = 0;
    603 
    604         if( m<=0 | n<=0 | zrows<=0 | zcolumns<=0 )
    605         {
    606             return;
    607         }
    608         System.Diagnostics.Debug.Assert(fromtheright & zcolumns==n | !fromtheright & zrows==n, "RMatrixBDMultiplyByP: incorrect Z size!");
    609        
    610         //
    611         // init
    612         //
    613         mx = Math.Max(m, n);
    614         mx = Math.Max(mx, zrows);
    615         mx = Math.Max(mx, zcolumns);
    616         v = new double[mx+1];
    617         work = new double[mx+1];
    618         v = new double[mx+1];
    619         work = new double[mx+1];
    620         if( m>=n )
    621         {
    622            
    623             //
    624             // setup
    625             //
    626             if( fromtheright )
    627             {
    628                 i1 = n-2;
    629                 i2 = 0;
    630                 istep = -1;
     624                    i1 = 0;
     625                    i2 = n-2;
     626                    istep = +1;
     627                }
     628                if( !dotranspose )
     629                {
     630                    i = i1;
     631                    i1 = i2;
     632                    i2 = i;
     633                    istep = -istep;
     634                }
     635               
     636                //
     637                // Process
     638                //
     639                if( n-1>0 )
     640                {
     641                    i = i1;
     642                    do
     643                    {
     644                        i1_ = (i+1) - (1);
     645                        for(i_=1; i_<=n-1-i;i_++)
     646                        {
     647                            v[i_] = qp[i,i_+i1_];
     648                        }
     649                        v[1] = 1;
     650                        if( fromtheright )
     651                        {
     652                            reflections.applyreflectionfromtheright(ref z, taup[i], ref v, 0, zrows-1, i+1, n-1, ref work);
     653                        }
     654                        else
     655                        {
     656                            reflections.applyreflectionfromtheleft(ref z, taup[i], ref v, i+1, n-1, 0, zcolumns-1, ref work);
     657                        }
     658                        i = i+istep;
     659                    }
     660                    while( i!=i2+istep );
     661                }
    631662            }
    632663            else
    633664            {
    634                 i1 = 0;
    635                 i2 = n-2;
    636                 istep = +1;
    637             }
    638             if( !dotranspose )
    639             {
    640                 i = i1;
    641                 i1 = i2;
    642                 i2 = i;
    643                 istep = -istep;
    644             }
    645            
    646             //
    647             // Process
    648             //
    649             if( n-1>0 )
    650             {
     665               
     666                //
     667                // setup
     668                //
     669                if( fromtheright )
     670                {
     671                    i1 = m-1;
     672                    i2 = 0;
     673                    istep = -1;
     674                }
     675                else
     676                {
     677                    i1 = 0;
     678                    i2 = m-1;
     679                    istep = +1;
     680                }
     681                if( !dotranspose )
     682                {
     683                    i = i1;
     684                    i1 = i2;
     685                    i2 = i;
     686                    istep = -istep;
     687                }
     688               
     689                //
     690                // Process
     691                //
    651692                i = i1;
    652693                do
    653694                {
    654                     i1_ = (i+1) - (1);
    655                     for(i_=1; i_<=n-1-i;i_++)
     695                    i1_ = (i) - (1);
     696                    for(i_=1; i_<=n-i;i_++)
    656697                    {
    657698                        v[i_] = qp[i,i_+i1_];
    … …  
    660701                    if( fromtheright )
    661702                    {
    662                         reflections.applyreflectionfromtheright(ref z, taup[i], ref v, 0, zrows-1, i+1, n-1, ref work);
     703                        reflections.applyreflectionfromtheright(ref z, taup[i], ref v, 0, zrows-1, i, n-1, ref work);
    663704                    }
    664705                    else
    665706                    {
    666                         reflections.applyreflectionfromtheleft(ref z, taup[i], ref v, i+1, n-1, 0, zcolumns-1, ref work);
     707                        reflections.applyreflectionfromtheleft(ref z, taup[i], ref v, i, n-1, 0, zcolumns-1, ref work);
    667708                    }
    668709                    i = i+istep;
    … …  
    671712            }
    672713        }
    673         else
     714
     715
     716        /*************************************************************************
     717        Unpacking of the main and secondary diagonals of bidiagonal decomposition
     718        of matrix A.
     719
     720        Input parameters:
     721            B   -   output of RMatrixBD subroutine.
     722            M   -   number of rows in matrix B.
     723            N   -   number of columns in matrix B.
     724
     725        Output parameters:
     726            IsUpper -   True, if the matrix is upper bidiagonal.
     727                        otherwise IsUpper is False.
     728            D       -   the main diagonal.
     729                        Array whose index ranges within [0..Min(M,N)-1].
     730            E       -   the secondary diagonal (upper or lower, depending on
     731                        the value of IsUpper).
     732                        Array index ranges within [0..Min(M,N)-1], the last
     733                        element is not used.
     734
     735          -- ALGLIB --
     736             Copyright 2005-2007 by Bochkanov Sergey
     737        *************************************************************************/
     738        public static void rmatrixbdunpackdiagonals(ref double[,] b,
     739            int m,
     740            int n,
     741            ref bool isupper,
     742            ref double[] d,
     743            ref double[] e)
    674744        {
    675            
    676             //
    677             // setup
    678             //
    679             if( fromtheright )
    680             {
    681                 i1 = m-1;
    682                 i2 = 0;
    683                 istep = -1;
     745            int i = 0;
     746
     747            isupper = m>=n;
     748            if( m<=0 | n<=0 )
     749            {
     750                return;
     751            }
     752            if( isupper )
     753            {
     754                d = new double[n-1+1];
     755                e = new double[n-1+1];
     756                for(i=0; i<=n-2; i++)
     757                {
     758                    d[i] = b[i,i];
     759                    e[i] = b[i,i+1];
     760                }
     761                d[n-1] = b[n-1,n-1];
    684762            }
    685763            else
    686764            {
    687                 i1 = 0;
    688                 i2 = m-1;
    689                 istep = +1;
    690             }
    691             if( !dotranspose )
    692             {
    693                 i = i1;
    694                 i1 = i2;
    695                 i2 = i;
    696                 istep = -istep;
    697             }
    698            
    699             //
    700             // Process
    701             //
    702             i = i1;
    703             do
    704             {
    705                 i1_ = (i) - (1);
    706                 for(i_=1; i_<=n-i;i_++)
    707                 {
    708                     v[i_] = qp[i,i_+i1_];
    709                 }
    710                 v[1] = 1;
    711                 if( fromtheright )
    712                 {
    713                     reflections.applyreflectionfromtheright(ref z, taup[i], ref v, 0, zrows-1, i, n-1, ref work);
    714                 }
    715                 else
    716                 {
    717                     reflections.applyreflectionfromtheleft(ref z, taup[i], ref v, i, n-1, 0, zcolumns-1, ref work);
    718                 }
    719                 i = i+istep;
    720             }
    721             while( i!=i2+istep );
     765                d = new double[m-1+1];
     766                e = new double[m-1+1];
     767                for(i=0; i<=m-2; i++)
     768                {
     769                    d[i] = b[i,i];
     770                    e[i] = b[i+1,i];
     771                }
     772                d[m-1] = b[m-1,m-1];
     773            }
    722774        }
    723     }
    724 
    725 
    726     /*************************************************************************
    727     Unpacking of the main and secondary diagonals of bidiagonal decomposition
    728     of matrix A.
    729 
    730     Input parameters:
    731         B   -   output of RMatrixBD subroutine.
    732         M   -   number of rows in matrix B.
    733         N   -   number of columns in matrix B.
    734 
    735     Output parameters:
    736         IsUpper -   True, if the matrix is upper bidiagonal.
    737                     otherwise IsUpper is False.
    738         D       -   the main diagonal.
    739                     Array whose index ranges within [0..Min(M,N)-1].
    740         E       -   the secondary diagonal (upper or lower, depending on
    741                     the value of IsUpper).
    742                     Array index ranges within [0..Min(M,N)-1], the last
    743                     element is not used.
    744 
    745       -- ALGLIB --
    746          Copyright 2005-2007 by Bochkanov Sergey
    747     *************************************************************************/
    748     public static void rmatrixbdunpackdiagonals(ref double[,] b,
    749         int m,
    750         int n,
    751         ref bool isupper,
    752         ref double[] d,
    753         ref double[] e)
    754     {
    755         int i = 0;
    756 
    757         isupper = m>=n;
    758         if( m<=0 | n<=0 )
     775
     776
     777        /*************************************************************************
     778        Obsolete 1-based subroutine.
     779        See RMatrixBD for 0-based replacement.
     780        *************************************************************************/
     781        public static void tobidiagonal(ref double[,] a,
     782            int m,
     783            int n,
     784            ref double[] tauq,
     785            ref double[] taup)
    759786        {
    760             return;
    761         }
    762         if( isupper )
    763         {
    764             d = new double[n-1+1];
    765             e = new double[n-1+1];
    766             for(i=0; i<=n-2; i++)
    767             {
    768                 d[i] = b[i,i];
    769                 e[i] = b[i,i+1];
    770             }
    771             d[n-1] = b[n-1,n-1];
    772         }
    773         else
    774         {
    775             d = new double[m-1+1];
    776             e = new double[m-1+1];
    777             for(i=0; i<=m-2; i++)
    778             {
    779                 d[i] = b[i,i];
    780                 e[i] = b[i+1,i];
    781             }
    782             d[m-1] = b[m-1,m-1];
    783         }
    784     }
    785 
    786 
    787     /*************************************************************************
    788     Obsolete 1-based subroutine.
    789     See RMatrixBD for 0-based replacement.
    790     *************************************************************************/
    791     public static void tobidiagonal(ref double[,] a,
    792         int m,
    793         int n,
    794         ref double[] tauq,
    795         ref double[] taup)
    796     {
    797         double[] work = new double[0];
    798         double[] t = new double[0];
    799         int minmn = 0;
    800         int maxmn = 0;
    801         int i = 0;
    802         double ltau = 0;
    803         int mmip1 = 0;
    804         int nmi = 0;
    805         int ip1 = 0;
    806         int nmip1 = 0;
    807         int mmi = 0;
    808         int i_ = 0;
    809         int i1_ = 0;
    810 
    811         minmn = Math.Min(m, n);
    812         maxmn = Math.Max(m, n);
    813         work = new double[maxmn+1];
    814         t = new double[maxmn+1];
    815         taup = new double[minmn+1];
    816         tauq = new double[minmn+1];
    817         if( m>=n )
    818         {
    819            
    820             //
    821             // Reduce to upper bidiagonal form
    822             //
    823             for(i=1; i<=n; i++)
    824             {
    825                
    826                 //
    827                 // Generate elementary reflector H(i) to annihilate A(i+1:m,i)
    828                 //
    829                 mmip1 = m-i+1;
    830                 i1_ = (i) - (1);
    831                 for(i_=1; i_<=mmip1;i_++)
    832                 {
    833                     t[i_] = a[i_+i1_,i];
    834                 }
    835                 reflections.generatereflection(ref t, mmip1, ref ltau);
    836                 tauq[i] = ltau;
    837                 i1_ = (1) - (i);
    838                 for(i_=i; i_<=m;i_++)
    839                 {
    840                     a[i_,i] = t[i_+i1_];
    841                 }
    842                 t[1] = 1;
    843                
    844                 //
    845                 // Apply H(i) to A(i:m,i+1:n) from the left
    846                 //
    847                 reflections.applyreflectionfromtheleft(ref a, ltau, ref t, i, m, i+1, n, ref work);
    848                 if( i<n )
     787            double[] work = new double[0];
     788            double[] t = new double[0];
     789            int minmn = 0;
     790            int maxmn = 0;
     791            int i = 0;
     792            double ltau = 0;
     793            int mmip1 = 0;
     794            int nmi = 0;
     795            int ip1 = 0;
     796            int nmip1 = 0;
     797            int mmi = 0;
     798            int i_ = 0;
     799            int i1_ = 0;
     800
     801            minmn = Math.Min(m, n);
     802            maxmn = Math.Max(m, n);
     803            work = new double[maxmn+1];
     804            t = new double[maxmn+1];
     805            taup = new double[minmn+1];
     806            tauq = new double[minmn+1];
     807            if( m>=n )
     808            {
     809               
     810                //
     811                // Reduce to upper bidiagonal form
     812                //
     813                for(i=1; i<=n; i++)
    849814                {
    850815                   
    851816                    //
    852                     // Generate elementary reflector G(i) to annihilate
    853                     // A(i,i+2:n)
    854                     //
    855                     nmi = n-i;
    856                     ip1 = i+1;
    857                     i1_ = (ip1) - (1);
    858                     for(i_=1; i_<=nmi;i_++)
    859                     {
    860                         t[i_] = a[i,i_+i1_];
    861                     }
    862                     reflections.generatereflection(ref t, nmi, ref ltau);
    863                     taup[i] = ltau;
    864                     i1_ = (1) - (ip1);
    865                     for(i_=ip1; i_<=n;i_++)
    866                     {
    867                         a[i,i_] = t[i_+i1_];
     817                    // Generate elementary reflector H(i) to annihilate A(i+1:m,i)
     818                    //
     819                    mmip1 = m-i+1;
     820                    i1_ = (i) - (1);
     821                    for(i_=1; i_<=mmip1;i_++)
     822                    {
     823                        t[i_] = a[i_+i1_,i];
     824                    }
     825                    reflections.generatereflection(ref t, mmip1, ref ltau);
     826                    tauq[i] = ltau;
     827                    i1_ = (1) - (i);
     828                    for(i_=i; i_<=m;i_++)
     829                    {
     830                        a[i_,i] = t[i_+i1_];
    868831                    }
    869832                    t[1] = 1;
    870833                   
    871834                    //
    872                     // Apply G(i) to A(i+1:m,i+1:n) from the right
    873                     //
    874                     reflections.applyreflectionfromtheright(ref a, ltau, ref t, i+1, m, i+1, n, ref work);
    875                 }
    876                 else
    877                 {
    878                     taup[i] = 0;
    879                 }
    880             }
    881         }
    882         else
    883         {
    884            
    885             //
    886             // Reduce to lower bidiagonal form
    887             //
    888             for(i=1; i<=m; i++)
    889             {
    890                
    891                 //
    892                 // Generate elementary reflector G(i) to annihilate A(i,i+1:n)
    893                 //
    894                 nmip1 = n-i+1;
    895                 i1_ = (i) - (1);
    896                 for(i_=1; i_<=nmip1;i_++)
    897                 {
    898                     t[i_] = a[i,i_+i1_];
    899                 }
    900                 reflections.generatereflection(ref t, nmip1, ref ltau);
    901                 taup[i] = ltau;
    902                 i1_ = (1) - (i);
    903                 for(i_=i; i_<=n;i_++)
    904                 {
    905                     a[i,i_] = t[i_+i1_];
    906                 }
    907                 t[1] = 1;
    908                
    909                 //
    910                 // Apply G(i) to A(i+1:m,i:n) from the right
    911                 //
    912                 reflections.applyreflectionfromtheright(ref a, ltau, ref t, i+1, m, i, n, ref work);
    913                 if( i<m )
     835                    // Apply H(i) to A(i:m,i+1:n) from the left
     836                    //
     837                    reflections.applyreflectionfromtheleft(ref a, ltau, ref t, i, m, i+1, n, ref work);
     838                    if( i<n )
     839                    {
     840                       
     841                        //
     842                        // Generate elementary reflector G(i) to annihilate
     843                        // A(i,i+2:n)
     844                        //
     845                        nmi = n-i;
     846                        ip1 = i+1;
     847                        i1_ = (ip1) - (1);
     848                        for(i_=1; i_<=nmi;i_++)
     849                        {
     850                            t[i_] = a[i,i_+i1_];
     851                        }
     852                        reflections.generatereflection(ref t, nmi, ref ltau);
     853                        taup[i] = ltau;
     854                        i1_ = (1) - (ip1);
     855                        for(i_=ip1; i_<=n;i_++)
     856                        {
     857                            a[i,i_] = t[i_+i1_];
     858                        }
     859                        t[1] = 1;
     860                       
     861                        //
     862                        // Apply G(i) to A(i+1:m,i+1:n) from the right
     863                        //
     864                        reflections.applyreflectionfromtheright(ref a, ltau, ref t, i+1, m, i+1, n, ref work);
     865                    }
     866                    else
     867                    {
     868                        taup[i] = 0;
     869                    }
     870                }
     871            }
     872            else
     873            {
     874               
     875                //
     876                // Reduce to lower bidiagonal form
     877                //
     878                for(i=1; i<=m; i++)
    914879                {
    915880                   
    916881                    //
    917                     // Generate elementary reflector H(i) to annihilate
    918                     // A(i+2:m,i)
    919                     //
    920                     mmi = m-i;
    921                     ip1 = i+1;
    922                     i1_ = (ip1) - (1);
    923                     for(i_=1; i_<=mmi;i_++)
    924                     {
    925                         t[i_] = a[i_+i1_,i];
    926                     }
    927                     reflections.generatereflection(ref t, mmi, ref ltau);
    928                     tauq[i] = ltau;
    929                     i1_ = (1) - (ip1);
    930                     for(i_=ip1; i_<=m;i_++)
    931                     {
    932                         a[i_,i] = t[i_+i1_];
     882                    // Generate elementary reflector G(i) to annihilate A(i,i+1:n)
     883                    //
     884                    nmip1 = n-i+1;
     885                    i1_ = (i) - (1);
     886                    for(i_=1; i_<=nmip1;i_++)
     887                    {
     888                        t[i_] = a[i,i_+i1_];
     889                    }
     890                    reflections.generatereflection(ref t, nmip1, ref ltau);
     891                    taup[i] = ltau;
     892                    i1_ = (1) - (i);
     893                    for(i_=i; i_<=n;i_++)
     894                    {
     895                        a[i,i_] = t[i_+i1_];
    933896                    }
    934897                    t[1] = 1;
    935898                   
    936899                    //
    937                     // Apply H(i) to A(i+1:m,i+1:n) from the left
    938                     //
    939                     reflections.applyreflectionfromtheleft(ref a, ltau, ref t, i+1, m, i+1, n, ref work);
    940                 }
    941                 else
    942                 {
    943                     tauq[i] = 0;
     900                    // Apply G(i) to A(i+1:m,i:n) from the right
     901                    //
     902                    reflections.applyreflectionfromtheright(ref a, ltau, ref t, i+1, m, i, n, ref work);
     903                    if( i<m )
     904                    {
     905                       
     906                        //
     907                        // Generate elementary reflector H(i) to annihilate
     908                        // A(i+2:m,i)
     909                        //
     910                        mmi = m-i;
     911                        ip1 = i+1;
     912                        i1_ = (ip1) - (1);
     913                        for(i_=1; i_<=mmi;i_++)
     914                        {
     915                            t[i_] = a[i_+i1_,i];
     916                        }
     917                        reflections.generatereflection(ref t, mmi, ref ltau);
     918                        tauq[i] = ltau;
     919                        i1_ = (1) - (ip1);
     920                        for(i_=ip1; i_<=m;i_++)
     921                        {
     922                            a[i_,i] = t[i_+i1_];
     923                        }
     924                        t[1] = 1;
     925                       
     926                        //
     927                        // Apply H(i) to A(i+1:m,i+1:n) from the left
     928                        //
     929                        reflections.applyreflectionfromtheleft(ref a, ltau, ref t, i+1, m, i+1, n, ref work);
     930                    }
     931                    else
     932                    {
     933                        tauq[i] = 0;
     934                    }
    944935                }
    945936            }
    946937        }
    947     }
    948 
    949 
    950     /*************************************************************************
    951     Obsolete 1-based subroutine.
    952     See RMatrixBDUnpackQ for 0-based replacement.
    953     *************************************************************************/
    954     public static void unpackqfrombidiagonal(ref double[,] qp,
    955         int m,
    956         int n,
    957         ref double[] tauq,
    958         int qcolumns,
    959         ref double[,] q)
    960     {
    961         int i = 0;
    962         int j = 0;
    963         int ip1 = 0;
    964         double[] v = new double[0];
    965         double[] work = new double[0];
    966         int vm = 0;
    967         int i_ = 0;
    968         int i1_ = 0;
    969 
    970         System.Diagnostics.Debug.Assert(qcolumns<=m, "UnpackQFromBidiagonal: QColumns>M!");
    971         if( m==0 | n==0 | qcolumns==0 )
     938
     939
     940        /*************************************************************************
     941        Obsolete 1-based subroutine.
     942        See RMatrixBDUnpackQ for 0-based replacement.
     943        *************************************************************************/
     944        public static void unpackqfrombidiagonal(ref double[,] qp,
     945            int m,
     946            int n,
     947            ref double[] tauq,
     948            int qcolumns,
     949            ref double[,] q)
    972950        {
    973             return;
    974         }
    975        
    976         //
    977         // init
    978         //
    979         q = new double[m+1, qcolumns+1];
    980         v = new double[m+1];
    981         work = new double[qcolumns+1];
    982        
    983         //
    984         // prepare Q
    985         //
    986         for(i=1; i<=m; i++)
    987         {
    988             for(j=1; j<=qcolumns; j++)
    989             {
    990                 if( i==j )
    991                 {
    992                     q[i,j] = 1;
    993                 }
    994                 else
    995                 {
    996                     q[i,j] = 0;
    997                 }
    998             }
    999         }
    1000         if( m>=n )
    1001         {
    1002             for(i=Math.Min(n, qcolumns); i>=1; i--)
    1003             {
    1004                 vm = m-i+1;
    1005                 i1_ = (i) - (1);
    1006                 for(i_=1; i_<=vm;i_++)
    1007                 {
    1008                     v[i_] = qp[i_+i1_,i];
    1009                 }
    1010                 v[1] = 1;
    1011                 reflections.applyreflectionfromtheleft(ref q, tauq[i], ref v, i, m, 1, qcolumns, ref work);
    1012             }
    1013         }
    1014         else
    1015         {
    1016             for(i=Math.Min(m-1, qcolumns-1); i>=1; i--)
    1017             {
    1018                 vm = m-i;
    1019                 ip1 = i+1;
    1020                 i1_ = (ip1) - (1);
    1021                 for(i_=1; i_<=vm;i_++)
    1022                 {
    1023                     v[i_] = qp[i_+i1_,i];
    1024                 }
    1025                 v[1] = 1;
    1026                 reflections.applyreflectionfromtheleft(ref q, tauq[i], ref v, i+1, m, 1, qcolumns, ref work);
    1027             }
    1028         }
    1029     }
    1030 
    1031 
    1032     /*************************************************************************
    1033     Obsolete 1-based subroutine.
    1034     See RMatrixBDMultiplyByQ for 0-based replacement.
    1035     *************************************************************************/
    1036     public static void multiplybyqfrombidiagonal(ref double[,] qp,
    1037         int m,
    1038         int n,
    1039         ref double[] tauq,
    1040         ref double[,] z,
    1041         int zrows,
    1042         int zcolumns,
    1043         bool fromtheright,
    1044         bool dotranspose)
    1045     {
    1046         int i = 0;
    1047         int ip1 = 0;
    1048         int i1 = 0;
    1049         int i2 = 0;
    1050         int istep = 0;
    1051         double[] v = new double[0];
    1052         double[] work = new double[0];
    1053         int vm = 0;
    1054         int mx = 0;
    1055         int i_ = 0;
    1056         int i1_ = 0;
    1057 
    1058         if( m<=0 | n<=0 | zrows<=0 | zcolumns<=0 )
    1059         {
    1060             return;
    1061         }
    1062         System.Diagnostics.Debug.Assert(fromtheright & zcolumns==m | !fromtheright & zrows==m, "MultiplyByQFromBidiagonal: incorrect Z size!");
    1063        
    1064         //
    1065         // init
    1066         //
    1067         mx = Math.Max(m, n);
    1068         mx = Math.Max(mx, zrows);
    1069         mx = Math.Max(mx, zcolumns);
    1070         v = new double[mx+1];
    1071         work = new double[mx+1];
    1072         if( m>=n )
    1073         {
     951            int i = 0;
     952            int j = 0;
     953            int ip1 = 0;
     954            double[] v = new double[0];
     955            double[] work = new double[0];
     956            int vm = 0;
     957            int i_ = 0;
     958            int i1_ = 0;
     959
     960            System.Diagnostics.Debug.Assert(qcolumns<=m, "UnpackQFromBidiagonal: QColumns>M!");
     961            if( m==0 | n==0 | qcolumns==0 )
     962            {
     963                return;
     964            }
    1074965           
    1075966            //
    1076             // setup
    1077             //
    1078             if( fromtheright )
    1079             {
    1080                 i1 = 1;
    1081                 i2 = n;
    1082                 istep = +1;
     967            // init
     968            //
     969            q = new double[m+1, qcolumns+1];
     970            v = new double[m+1];
     971            work = new double[qcolumns+1];
     972           
     973            //
     974            // prepare Q
     975            //
     976            for(i=1; i<=m; i++)
     977            {
     978                for(j=1; j<=qcolumns; j++)
     979                {
     980                    if( i==j )
     981                    {
     982                        q[i,j] = 1;
     983                    }
     984                    else
     985                    {
     986                        q[i,j] = 0;
     987                    }
     988                }
     989            }
     990            if( m>=n )
     991            {
     992                for(i=Math.Min(n, qcolumns); i>=1; i--)
     993                {
     994                    vm = m-i+1;
     995                    i1_ = (i) - (1);
     996                    for(i_=1; i_<=vm;i_++)
     997                    {
     998                        v[i_] = qp[i_+i1_,i];
     999                    }
     1000                    v[1] = 1;
     1001                    reflections.applyreflectionfromtheleft(ref q, tauq[i], ref v, i, m, 1, qcolumns, ref work);
     1002                }
    10831003            }
    10841004            else
    10851005            {
    1086                 i1 = n;
    1087                 i2 = 1;
    1088                 istep = -1;
    1089             }
    1090             if( dotranspose )
    1091             {
    1092                 i = i1;
    1093                 i1 = i2;
    1094                 i2 = i;
    1095                 istep = -istep;
    1096             }
    1097            
    1098             //
    1099             // Process
    1100             //
    1101             i = i1;
    1102             do
    1103             {
    1104                 vm = m-i+1;
    1105                 i1_ = (i) - (1);
    1106                 for(i_=1; i_<=vm;i_++)
    1107                 {
    1108                     v[i_] = qp[i_+i1_,i];
    1109                 }
    1110                 v[1] = 1;
    1111                 if( fromtheright )
    1112                 {
    1113                     reflections.applyreflectionfromtheright(ref z, tauq[i], ref v, 1, zrows, i, m, ref work);
    1114                 }
    1115                 else
    1116                 {
    1117                     reflections.applyreflectionfromtheleft(ref z, tauq[i], ref v, i, m, 1, zcolumns, ref work);
    1118                 }
    1119                 i = i+istep;
    1120             }
    1121             while( i!=i2+istep );
    1122         }
    1123         else
    1124         {
    1125            
    1126             //
    1127             // setup
    1128             //
    1129             if( fromtheright )
    1130             {
    1131                 i1 = 1;
    1132                 i2 = m-1;
    1133                 istep = +1;
    1134             }
    1135             else
    1136             {
    1137                 i1 = m-1;
    1138                 i2 = 1;
    1139                 istep = -1;
    1140             }
    1141             if( dotranspose )
    1142             {
    1143                 i = i1;
    1144                 i1 = i2;
    1145                 i2 = i;
    1146                 istep = -istep;
    1147             }
    1148            
    1149             //
    1150             // Process
    1151             //
    1152             if( m-1>0 )
    1153             {
    1154                 i = i1;
    1155                 do
     1006                for(i=Math.Min(m-1, qcolumns-1); i>=1; i--)
    11561007                {
    11571008                    vm = m-i;
    … …  
    11631014                    }
    11641015                    v[1] = 1;
    1165                     if( fromtheright )
    1166                     {
    1167                         reflections.applyreflectionfromtheright(ref z, tauq[i], ref v, 1, zrows, i+1, m, ref work);
    1168                     }
    1169                     else
    1170                     {
    1171                         reflections.applyreflectionfromtheleft(ref z, tauq[i], ref v, i+1, m, 1, zcolumns, ref work);
    1172                     }
    1173                     i = i+istep;
    1174                 }
    1175                 while( i!=i2+istep );
     1016                    reflections.applyreflectionfromtheleft(ref q, tauq[i], ref v, i+1, m, 1, qcolumns, ref work);
     1017                }
    11761018            }
    11771019        }
    1178     }
    1179 
    1180 
    1181     /*************************************************************************
    1182     Obsolete 1-based subroutine.
    1183     See RMatrixBDUnpackPT for 0-based replacement.
    1184     *************************************************************************/
    1185     public static void unpackptfrombidiagonal(ref double[,] qp,
    1186         int m,
    1187         int n,
    1188         ref double[] taup,
    1189         int ptrows,
    1190         ref double[,] pt)
    1191     {
    1192         int i = 0;
    1193         int j = 0;
    1194         int ip1 = 0;
    1195         double[] v = new double[0];
    1196         double[] work = new double[0];
    1197         int vm = 0;
    1198         int i_ = 0;
    1199         int i1_ = 0;
    1200 
    1201         System.Diagnostics.Debug.Assert(ptrows<=n, "UnpackPTFromBidiagonal: PTRows>N!");
    1202         if( m==0 | n==0 | ptrows==0 )
     1020
     1021
     1022        /*************************************************************************
     1023        Obsolete 1-based subroutine.
     1024        See RMatrixBDMultiplyByQ for 0-based replacement.
     1025        *************************************************************************/
     1026        public static void multiplybyqfrombidiagonal(ref double[,] qp,
     1027            int m,
     1028            int n,
     1029            ref double[] tauq,
     1030            ref double[,] z,
     1031            int zrows,
     1032            int zcolumns,
     1033            bool fromtheright,
     1034            bool dotranspose)
    12031035        {
    1204             return;
    1205         }
    1206        
    1207         //
    1208         // init
    1209         //
    1210         pt = new double[ptrows+1, n+1];
    1211         v = new double[n+1];
    1212         work = new double[ptrows+1];
    1213        
    1214         //
    1215         // prepare PT
    1216         //
    1217         for(i=1; i<=ptrows; i++)
    1218         {
    1219             for(j=1; j<=n; j++)
    1220             {
    1221                 if( i==j )
    1222                 {
    1223                     pt[i,j] = 1;
     1036            int i = 0;
     1037            int ip1 = 0;
     1038            int i1 = 0;
     1039            int i2 = 0;
     1040            int istep = 0;
     1041            double[] v = new double[0];
     1042            double[] work = new double[0];
     1043            int vm = 0;
     1044            int mx = 0;
     1045            int i_ = 0;
     1046            int i1_ = 0;
     1047
     1048            if( m<=0 | n<=0 | zrows<=0 | zcolumns<=0 )
     1049            {
     1050                return;
     1051            }
     1052            System.Diagnostics.Debug.Assert(fromtheright & zcolumns==m | !fromtheright & zrows==m, "MultiplyByQFromBidiagonal: incorrect Z size!");
     1053           
     1054            //
     1055            // init
     1056            //
     1057            mx = Math.Max(m, n);
     1058            mx = Math.Max(mx, zrows);
     1059            mx = Math.Max(mx, zcolumns);
     1060            v = new double[mx+1];
     1061            work = new double[mx+1];
     1062            if( m>=n )
     1063            {
     1064               
     1065                //
     1066                // setup
     1067                //
     1068                if( fromtheright )
     1069                {
     1070                    i1 = 1;
     1071                    i2 = n;
     1072                    istep = +1;
    12241073                }
    12251074                else
    12261075                {
    1227                     pt[i,j] = 0;
    1228                 }
    1229             }
    1230         }
    1231         if( m>=n )
    1232         {
    1233             for(i=Math.Min(n-1, ptrows-1); i>=1; i--)
    1234             {
    1235                 vm = n-i;
    1236                 ip1 = i+1;
    1237                 i1_ = (ip1) - (1);
    1238                 for(i_=1; i_<=vm;i_++)
    1239                 {
    1240                     v[i_] = qp[i,i_+i1_];
    1241                 }
    1242                 v[1] = 1;
    1243                 reflections.applyreflectionfromtheright(ref pt, taup[i], ref v, 1, ptrows, i+1, n, ref work);
    1244             }
    1245         }
    1246         else
    1247         {
    1248             for(i=Math.Min(m, ptrows); i>=1; i--)
    1249             {
    1250                 vm = n-i+1;
    1251                 i1_ = (i) - (1);
    1252                 for(i_=1; i_<=vm;i_++)
    1253                 {
    1254                     v[i_] = qp[i,i_+i1_];
    1255                 }
    1256                 v[1] = 1;
    1257                 reflections.applyreflectionfromtheright(ref pt, taup[i], ref v, 1, ptrows, i, n, ref work);
    1258             }
    1259         }
    1260     }
    1261 
    1262 
    1263     /*************************************************************************
    1264     Obsolete 1-based subroutine.
    1265     See RMatrixBDMultiplyByP for 0-based replacement.
    1266     *************************************************************************/
    1267     public static void multiplybypfrombidiagonal(ref double[,] qp,
    1268         int m,
    1269         int n,
    1270         ref double[] taup,
    1271         ref double[,] z,
    1272         int zrows,
    1273         int zcolumns,
    1274         bool fromtheright,
    1275         bool dotranspose)
    1276     {
    1277         int i = 0;
    1278         int ip1 = 0;
    1279         double[] v = new double[0];
    1280         double[] work = new double[0];
    1281         int vm = 0;
    1282         int mx = 0;
    1283         int i1 = 0;
    1284         int i2 = 0;
    1285         int istep = 0;
    1286         int i_ = 0;
    1287         int i1_ = 0;
    1288 
    1289         if( m<=0 | n<=0 | zrows<=0 | zcolumns<=0 )
    1290         {
    1291             return;
    1292         }
    1293         System.Diagnostics.Debug.Assert(fromtheright & zcolumns==n | !fromtheright & zrows==n, "MultiplyByQFromBidiagonal: incorrect Z size!");
    1294        
    1295         //
    1296         // init
    1297         //
    1298         mx = Math.Max(m, n);
    1299         mx = Math.Max(mx, zrows);
    1300         mx = Math.Max(mx, zcolumns);
    1301         v = new double[mx+1];
    1302         work = new double[mx+1];
    1303         v = new double[mx+1];
    1304         work = new double[mx+1];
    1305         if( m>=n )
    1306         {
    1307            
    1308             //
    1309             // setup
    1310             //
    1311             if( fromtheright )
    1312             {
    1313                 i1 = n-1;
    1314                 i2 = 1;
    1315                 istep = -1;
    1316             }
    1317             else
    1318             {
    1319                 i1 = 1;
    1320                 i2 = n-1;
    1321                 istep = +1;
    1322             }
    1323             if( !dotranspose )
    1324             {
    1325                 i = i1;
    1326                 i1 = i2;
    1327                 i2 = i;
    1328                 istep = -istep;
    1329             }
    1330            
    1331             //
    1332             // Process
    1333             //
    1334             if( n-1>0 )
    1335             {
     1076                    i1 = n;
     1077                    i2 = 1;
     1078                    istep = -1;
     1079                }
     1080                if( dotranspose )
     1081                {
     1082                    i = i1;
     1083                    i1 = i2;
     1084                    i2 = i;
     1085                    istep = -istep;
     1086                }
     1087               
     1088                //
     1089                // Process
     1090                //
    13361091                i = i1;
    13371092                do
     1093                {
     1094                    vm = m-i+1;
     1095                    i1_ = (i) - (1);
     1096                    for(i_=1; i_<=vm;i_++)
     1097                    {
     1098                        v[i_] = qp[i_+i1_,i];
     1099                    }
     1100                    v[1] = 1;
     1101                    if( fromtheright )
     1102                    {
     1103                        reflections.applyreflectionfromtheright(ref z, tauq[i], ref v, 1, zrows, i, m, ref work);
     1104                    }
     1105                    else
     1106                    {
     1107                        reflections.applyreflectionfromtheleft(ref z, tauq[i], ref v, i, m, 1, zcolumns, ref work);
     1108                    }
     1109                    i = i+istep;
     1110                }
     1111                while( i!=i2+istep );
     1112            }
     1113            else
     1114            {
     1115               
     1116                //
     1117                // setup
     1118                //
     1119                if( fromtheright )
     1120                {
     1121                    i1 = 1;
     1122                    i2 = m-1;
     1123                    istep = +1;
     1124                }
     1125                else
     1126                {
     1127                    i1 = m-1;
     1128                    i2 = 1;
     1129                    istep = -1;
     1130                }
     1131                if( dotranspose )
     1132                {
     1133                    i = i1;
     1134                    i1 = i2;
     1135                    i2 = i;
     1136                    istep = -istep;
     1137                }
     1138               
     1139                //
     1140                // Process
     1141                //
     1142                if( m-1>0 )
     1143                {
     1144                    i = i1;
     1145                    do
     1146                    {
     1147                        vm = m-i;
     1148                        ip1 = i+1;
     1149                        i1_ = (ip1) - (1);
     1150                        for(i_=1; i_<=vm;i_++)
     1151                        {
     1152                            v[i_] = qp[i_+i1_,i];
     1153                        }
     1154                        v[1] = 1;
     1155                        if( fromtheright )
     1156                        {
     1157                            reflections.applyreflectionfromtheright(ref z, tauq[i], ref v, 1, zrows, i+1, m, ref work);
     1158                        }
     1159                        else
     1160                        {
     1161                            reflections.applyreflectionfromtheleft(ref z, tauq[i], ref v, i+1, m, 1, zcolumns, ref work);
     1162                        }
     1163                        i = i+istep;
     1164                    }
     1165                    while( i!=i2+istep );
     1166                }
     1167            }
     1168        }
     1169
     1170
     1171        /*************************************************************************
     1172        Obsolete 1-based subroutine.
     1173        See RMatrixBDUnpackPT for 0-based replacement.
     1174        *************************************************************************/
     1175        public static void unpackptfrombidiagonal(ref double[,] qp,
     1176            int m,
     1177            int n,
     1178            ref double[] taup,
     1179            int ptrows,
     1180            ref double[,] pt)
     1181        {
     1182            int i = 0;
     1183            int j = 0;
     1184            int ip1 = 0;
     1185            double[] v = new double[0];
     1186            double[] work = new double[0];
     1187            int vm = 0;
     1188            int i_ = 0;
     1189            int i1_ = 0;
     1190
     1191            System.Diagnostics.Debug.Assert(ptrows<=n, "UnpackPTFromBidiagonal: PTRows>N!");
     1192            if( m==0 | n==0 | ptrows==0 )
     1193            {
     1194                return;
     1195            }
     1196           
     1197            //
     1198            // init
     1199            //
     1200            pt = new double[ptrows+1, n+1];
     1201            v = new double[n+1];
     1202            work = new double[ptrows+1];
     1203           
     1204            //
     1205            // prepare PT
     1206            //
     1207            for(i=1; i<=ptrows; i++)
     1208            {
     1209                for(j=1; j<=n; j++)
     1210                {
     1211                    if( i==j )
     1212                    {
     1213                        pt[i,j] = 1;
     1214                    }
     1215                    else
     1216                    {
     1217                        pt[i,j] = 0;
     1218                    }
     1219                }
     1220            }
     1221            if( m>=n )
     1222            {
     1223                for(i=Math.Min(n-1, ptrows-1); i>=1; i--)
    13381224                {
    13391225                    vm = n-i;
    … …  
    13451231                    }
    13461232                    v[1] = 1;
     1233                    reflections.applyreflectionfromtheright(ref pt, taup[i], ref v, 1, ptrows, i+1, n, ref work);
     1234                }
     1235            }
     1236            else
     1237            {
     1238                for(i=Math.Min(m, ptrows); i>=1; i--)
     1239                {
     1240                    vm = n-i+1;
     1241                    i1_ = (i) - (1);
     1242                    for(i_=1; i_<=vm;i_++)
     1243                    {
     1244                        v[i_] = qp[i,i_+i1_];
     1245                    }
     1246                    v[1] = 1;
     1247                    reflections.applyreflectionfromtheright(ref pt, taup[i], ref v, 1, ptrows, i, n, ref work);
     1248                }
     1249            }
     1250        }
     1251
     1252
     1253        /*************************************************************************
     1254        Obsolete 1-based subroutine.
     1255        See RMatrixBDMultiplyByP for 0-based replacement.
     1256        *************************************************************************/
     1257        public static void multiplybypfrombidiagonal(ref double[,] qp,
     1258            int m,
     1259            int n,
     1260            ref double[] taup,
     1261            ref double[,] z,
     1262            int zrows,
     1263            int zcolumns,
     1264            bool fromtheright,
     1265            bool dotranspose)
     1266        {
     1267            int i = 0;
     1268            int ip1 = 0;
     1269            double[] v = new double[0];
     1270            double[] work = new double[0];
     1271            int vm = 0;
     1272            int mx = 0;
     1273            int i1 = 0;
     1274            int i2 = 0;
     1275            int istep = 0;
     1276            int i_ = 0;
     1277            int i1_ = 0;
     1278
     1279            if( m<=0 | n<=0 | zrows<=0 | zcolumns<=0 )
     1280            {
     1281                return;
     1282            }
     1283            System.Diagnostics.Debug.Assert(fromtheright & zcolumns==n | !fromtheright & zrows==n, "MultiplyByQFromBidiagonal: incorrect Z size!");
     1284           
     1285            //
     1286            // init
     1287            //
     1288            mx = Math.Max(m, n);
     1289            mx = Math.Max(mx, zrows);
     1290            mx = Math.Max(mx, zcolumns);
     1291            v = new double[mx+1];
     1292            work = new double[mx+1];
     1293            v = new double[mx+1];
     1294            work = new double[mx+1];
     1295            if( m>=n )
     1296            {
     1297               
     1298                //
     1299                // setup
     1300                //
     1301                if( fromtheright )
     1302                {
     1303                    i1 = n-1;
     1304                    i2 = 1;
     1305                    istep = -1;
     1306                }
     1307                else
     1308                {
     1309                    i1 = 1;
     1310                    i2 = n-1;
     1311                    istep = +1;
     1312                }
     1313                if( !dotranspose )
     1314                {
     1315                    i = i1;
     1316                    i1 = i2;
     1317                    i2 = i;
     1318                    istep = -istep;
     1319                }
     1320               
     1321                //
     1322                // Process
     1323                //
     1324                if( n-1>0 )
     1325                {
     1326                    i = i1;
     1327                    do
     1328                    {
     1329                        vm = n-i;
     1330                        ip1 = i+1;
     1331                        i1_ = (ip1) - (1);
     1332                        for(i_=1; i_<=vm;i_++)
     1333                        {
     1334                            v[i_] = qp[i,i_+i1_];
     1335                        }
     1336                        v[1] = 1;
     1337                        if( fromtheright )
     1338                        {
     1339                            reflections.applyreflectionfromtheright(ref z, taup[i], ref v, 1, zrows, i+1, n, ref work);
     1340                        }
     1341                        else
     1342                        {
     1343                            reflections.applyreflectionfromtheleft(ref z, taup[i], ref v, i+1, n, 1, zcolumns, ref work);
     1344                        }
     1345                        i = i+istep;
     1346                    }
     1347                    while( i!=i2+istep );
     1348                }
     1349            }
     1350            else
     1351            {
     1352               
     1353                //
     1354                // setup
     1355                //
     1356                if( fromtheright )
     1357                {
     1358                    i1 = m;
     1359                    i2 = 1;
     1360                    istep = -1;
     1361                }
     1362                else
     1363                {
     1364                    i1 = 1;
     1365                    i2 = m;
     1366                    istep = +1;
     1367                }
     1368                if( !dotranspose )
     1369                {
     1370                    i = i1;
     1371                    i1 = i2;
     1372                    i2 = i;
     1373                    istep = -istep;
     1374                }
     1375               
     1376                //
     1377                // Process
     1378                //
     1379                i = i1;
     1380                do
     1381                {
     1382                    vm = n-i+1;
     1383                    i1_ = (i) - (1);
     1384                    for(i_=1; i_<=vm;i_++)
     1385                    {
     1386                        v[i_] = qp[i,i_+i1_];
     1387                    }
     1388                    v[1] = 1;
    13471389                    if( fromtheright )
    13481390                    {
    1349                         reflections.applyreflectionfromtheright(ref z, taup[i], ref v, 1, zrows, i+1, n, ref work);
     1391                        reflections.applyreflectionfromtheright(ref z, taup[i], ref v, 1, zrows, i, n, ref work);
    13501392                    }
    13511393                    else
    13521394                    {
    1353                         reflections.applyreflectionfromtheleft(ref z, taup[i], ref v, i+1, n, 1, zcolumns, ref work);
     1395                        reflections.applyreflectionfromtheleft(ref z, taup[i], ref v, i, n, 1, zcolumns, ref work);
    13541396                    }
    13551397                    i = i+istep;
    … …  
    13581400            }
    13591401        }
    1360         else
     1402
     1403
     1404        /*************************************************************************
     1405        Obsolete 1-based subroutine.
     1406        See RMatrixBDUnpackDiagonals for 0-based replacement.
     1407        *************************************************************************/
     1408        public static void unpackdiagonalsfrombidiagonal(ref double[,] b,
     1409            int m,
     1410            int n,
     1411            ref bool isupper,
     1412            ref double[] d,
     1413            ref double[] e)
    13611414        {
    1362            
    1363             //
    1364             // setup
    1365             //
    1366             if( fromtheright )
    1367             {
    1368                 i1 = m;
    1369                 i2 = 1;
    1370                 istep = -1;
     1415            int i = 0;
     1416
     1417            isupper = m>=n;
     1418            if( m==0 | n==0 )
     1419            {
     1420                return;
     1421            }
     1422            if( isupper )
     1423            {
     1424                d = new double[n+1];
     1425                e = new double[n+1];
     1426                for(i=1; i<=n-1; i++)
     1427                {
     1428                    d[i] = b[i,i];
     1429                    e[i] = b[i,i+1];
     1430                }
     1431                d[n] = b[n,n];
    13711432            }
    13721433            else
    13731434            {
    1374                 i1 = 1;
    1375                 i2 = m;
    1376                 istep = +1;
    1377             }
    1378             if( !dotranspose )
    1379             {
    1380                 i = i1;
    1381                 i1 = i2;
    1382                 i2 = i;
    1383                 istep = -istep;
    1384             }
    1385            
    1386             //
    1387             // Process
    1388             //
    1389             i = i1;
    1390             do
    1391             {
    1392                 vm = n-i+1;
    1393                 i1_ = (i) - (1);
    1394                 for(i_=1; i_<=vm;i_++)
    1395                 {
    1396                     v[i_] = qp[i,i_+i1_];
    1397                 }
    1398                 v[1] = 1;
    1399                 if( fromtheright )
    1400                 {
    1401                     reflections.applyreflectionfromtheright(ref z, taup[i], ref v, 1, zrows, i, n, ref work);
    1402                 }
    1403                 else
    1404                 {
    1405                     reflections.applyreflectionfromtheleft(ref z, taup[i], ref v, i, n, 1, zcolumns, ref work);
    1406                 }
    1407                 i = i+istep;
    1408             }
    1409             while( i!=i2+istep );
    1410         }
    1411     }
    1412 
    1413 
    1414     /*************************************************************************
    1415     Obsolete 1-based subroutine.
    1416     See RMatrixBDUnpackDiagonals for 0-based replacement.
    1417     *************************************************************************/
    1418     public static void unpackdiagonalsfrombidiagonal(ref double[,] b,
    1419         int m,
    1420         int n,
    1421         ref bool isupper,
    1422         ref double[] d,
    1423         ref double[] e)
    1424     {
    1425         int i = 0;
    1426 
    1427         isupper = m>=n;
    1428         if( m==0 | n==0 )
    1429         {
    1430             return;
    1431         }
    1432         if( isupper )
    1433         {
    1434             d = new double[n+1];
    1435             e = new double[n+1];
    1436             for(i=1; i<=n-1; i++)
    1437             {
    1438                 d[i] = b[i,i];
    1439                 e[i] = b[i,i+1];
    1440             }
    1441             d[n] = b[n,n];
    1442         }
    1443         else
    1444         {
    1445             d = new double[m+1];
    1446             e = new double[m+1];
    1447             for(i=1; i<=m-1; i++)
    1448             {
    1449                 d[i] = b[i,i];
    1450                 e[i] = b[i+1,i];
    1451             }
    1452             d[m] = b[m,m];
     1435                d = new double[m+1];
     1436                e = new double[m+1];
     1437                for(i=1; i<=m-1; i++)
     1438                {
     1439                    d[i] = b[i,i];
     1440                    e[i] = b[i+1,i];
     1441                }
     1442                d[m] = b[m,m];
     1443            }
    14531444        }
    14541445    }
Note: See TracChangeset for help on using the changeset viewer.