Free cookie consent management tool by TermsFeed Policy Generator

source: trunk/sources/HeuristicLab.LinearRegression/3.2/alglib/bidiagonal.cs @ 2154

Last change on this file since 2154 was 2154, checked in by gkronber, 15 years ago

Added linear regression plugin. #697

File size: 44.0 KB
Line 
1/*************************************************************************
2Copyright (c) 1992-2007 The University of Tennessee.  All rights reserved.
3
4Contributors:
5    * Sergey Bochkanov (ALGLIB project). Translation from FORTRAN to
6      pseudocode.
7
8See subroutines comments for additional copyrights.
9
10Redistribution and use in source and binary forms, with or without
11modification, are permitted provided that the following conditions are
12met:
13
14- Redistributions of source code must retain the above copyright
15  notice, this list of conditions and the following disclaimer.
16
17- Redistributions in binary form must reproduce the above copyright
18  notice, this list of conditions and the following disclaimer listed
19  in this license in the documentation and/or other materials
20  provided with the distribution.
21
22- Neither the name of the copyright holders nor the names of its
23  contributors may be used to endorse or promote products derived from
24  this software without specific prior written permission.
25
26THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
27"AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
28LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
29A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT
30OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
31SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
32LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
33DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
34THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
35(INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
36OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
37*************************************************************************/
38
39using System;
40
41class bidiagonal
42{
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)
98    {
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 )
114        {
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        {
133           
134            //
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 )
162                {
163                   
164                    //
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_];
179                    }
180                    t[1] = 1;
181                   
182                    //
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 )
224                {
225                   
226                    //
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_];
241                    }
242                    t[1] = 1;
243                   
244                    //
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);
248                }
249                else
250                {
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            {
466                i = i1;
467                do
468                {
469                    i1_ = (i+1) - (1);
470                    for(i_=1; i_<=m-i-1;i_++)
471                    {
472                        v[i_] = qp[i_+i1_,i];
473                    }
474                    v[1] = 1;
475                    if( fromtheright )
476                    {
477                        reflections.applyreflectionfromtheright(ref z, tauq[i], ref v, 0, zrows-1, i+1, m-1, ref work);
478                    }
479                    else
480                    {
481                        reflections.applyreflectionfromtheleft(ref z, tauq[i], ref v, i+1, m-1, 0, zcolumns-1, ref work);
482                    }
483                    i = i+istep;
484                }
485                while( i!=i2+istep );
486            }
487        }
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 )
525        {
526            return;
527        }
528       
529        //
530        // prepare PT
531        //
532        pt = new double[ptrows-1+1, n-1+1];
533        for(i=0; i<=ptrows-1; i++)
534        {
535            for(j=0; j<=n-1; j++)
536            {
537                if( i==j )
538                {
539                    pt[i,j] = 1;
540                }
541                else
542                {
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;
631            }
632            else
633            {
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            {
651                i = i1;
652                do
653                {
654                    i1_ = (i+1) - (1);
655                    for(i_=1; i_<=n-1-i;i_++)
656                    {
657                        v[i_] = qp[i,i_+i1_];
658                    }
659                    v[1] = 1;
660                    if( fromtheright )
661                    {
662                        reflections.applyreflectionfromtheright(ref z, taup[i], ref v, 0, zrows-1, i+1, n-1, ref work);
663                    }
664                    else
665                    {
666                        reflections.applyreflectionfromtheleft(ref z, taup[i], ref v, i+1, n-1, 0, zcolumns-1, ref work);
667                    }
668                    i = i+istep;
669                }
670                while( i!=i2+istep );
671            }
672        }
673        else
674        {
675           
676            //
677            // setup
678            //
679            if( fromtheright )
680            {
681                i1 = m-1;
682                i2 = 0;
683                istep = -1;
684            }
685            else
686            {
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 );
722        }
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 )
759        {
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 )
849                {
850                   
851                    //
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_];
868                    }
869                    t[1] = 1;
870                   
871                    //
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 )
914                {
915                   
916                    //
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_];
933                    }
934                    t[1] = 1;
935                   
936                    //
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;
944                }
945            }
946        }
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 )
972        {
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        {
1074           
1075            //
1076            // setup
1077            //
1078            if( fromtheright )
1079            {
1080                i1 = 1;
1081                i2 = n;
1082                istep = +1;
1083            }
1084            else
1085            {
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
1156                {
1157                    vm = m-i;
1158                    ip1 = i+1;
1159                    i1_ = (ip1) - (1);
1160                    for(i_=1; i_<=vm;i_++)
1161                    {
1162                        v[i_] = qp[i_+i1_,i];
1163                    }
1164                    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 );
1176            }
1177        }
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 )
1203        {
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;
1224                }
1225                else
1226                {
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            {
1336                i = i1;
1337                do
1338                {
1339                    vm = n-i;
1340                    ip1 = i+1;
1341                    i1_ = (ip1) - (1);
1342                    for(i_=1; i_<=vm;i_++)
1343                    {
1344                        v[i_] = qp[i,i_+i1_];
1345                    }
1346                    v[1] = 1;
1347                    if( fromtheright )
1348                    {
1349                        reflections.applyreflectionfromtheright(ref z, taup[i], ref v, 1, zrows, i+1, n, ref work);
1350                    }
1351                    else
1352                    {
1353                        reflections.applyreflectionfromtheleft(ref z, taup[i], ref v, i+1, n, 1, zcolumns, ref work);
1354                    }
1355                    i = i+istep;
1356                }
1357                while( i!=i2+istep );
1358            }
1359        }
1360        else
1361        {
1362           
1363            //
1364            // setup
1365            //
1366            if( fromtheright )
1367            {
1368                i1 = m;
1369                i2 = 1;
1370                istep = -1;
1371            }
1372            else
1373            {
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];
1453        }
1454    }
1455}
Note: See TracBrowser for help on using the repository browser.