Free cookie consent management tool by TermsFeed Policy Generator

source: branches/Persistence Test/ALGLIB/bidiagonal.cs @ 3579

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

Moved ALGLIB code into a separate plugin. #783
File size: 48.5 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
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 >>>
25*************************************************************************/
26
27using System;
28
29namespace alglib
30{
31    public class bidiagonal
32    {
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)
88        {
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
99           
100            //
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++)
128                {
129                   
130                    //
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_];
144                    }
145                    t[1] = 1;
146                   
147                    //
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++)
190                {
191                   
192                    //
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_];
206                    }
207                    t[1] = 1;
208                   
209                    //
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;
386                }
387                else
388                {
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                //
404                i = i1;
405                do
406                {
407                    i1_ = (i) - (1);
408                    for(i_=1; i_<=m-i;i_++)
409                    {
410                        v[i_] = qp[i_+i1_,i];
411                    }
412                    v[1] = 1;
413                    if( fromtheright )
414                    {
415                        reflections.applyreflectionfromtheright(ref z, tauq[i], ref v, 0, zrows-1, i, m-1, ref work);
416                    }
417                    else
418                    {
419                        reflections.applyreflectionfromtheleft(ref z, tauq[i], ref v, i, m-1, 0, zcolumns-1, ref work);
420                    }
421                    i = i+istep;
422                }
423                while( i!=i2+istep );
424            }
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            }
478        }
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)
508        {
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);
542        }
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)
583        {
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;
621                }
622                else
623                {
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                }
662            }
663            else
664            {
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                //
692                i = i1;
693                do
694                {
695                    i1_ = (i) - (1);
696                    for(i_=1; i_<=n-i;i_++)
697                    {
698                        v[i_] = qp[i,i_+i1_];
699                    }
700                    v[1] = 1;
701                    if( fromtheright )
702                    {
703                        reflections.applyreflectionfromtheright(ref z, taup[i], ref v, 0, zrows-1, i, n-1, ref work);
704                    }
705                    else
706                    {
707                        reflections.applyreflectionfromtheleft(ref z, taup[i], ref v, i, n-1, 0, zcolumns-1, ref work);
708                    }
709                    i = i+istep;
710                }
711                while( i!=i2+istep );
712            }
713        }
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)
744        {
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];
762            }
763            else
764            {
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            }
774        }
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)
786        {
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++)
814                {
815                   
816                    //
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_];
831                    }
832                    t[1] = 1;
833                   
834                    //
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++)
879                {
880                   
881                    //
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_];
896                    }
897                    t[1] = 1;
898                   
899                    //
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                    }
935                }
936            }
937        }
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)
950        {
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            }
965           
966            //
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                }
1003            }
1004            else
1005            {
1006                for(i=Math.Min(m-1, qcolumns-1); i>=1; i--)
1007                {
1008                    vm = m-i;
1009                    ip1 = i+1;
1010                    i1_ = (ip1) - (1);
1011                    for(i_=1; i_<=vm;i_++)
1012                    {
1013                        v[i_] = qp[i_+i1_,i];
1014                    }
1015                    v[1] = 1;
1016                    reflections.applyreflectionfromtheleft(ref q, tauq[i], ref v, i+1, m, 1, qcolumns, ref work);
1017                }
1018            }
1019        }
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)
1035        {
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;
1073                }
1074                else
1075                {
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                //
1091                i = i1;
1092                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--)
1224                {
1225                    vm = n-i;
1226                    ip1 = i+1;
1227                    i1_ = (ip1) - (1);
1228                    for(i_=1; i_<=vm;i_++)
1229                    {
1230                        v[i_] = qp[i,i_+i1_];
1231                    }
1232                    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;
1389                    if( fromtheright )
1390                    {
1391                        reflections.applyreflectionfromtheright(ref z, taup[i], ref v, 1, zrows, i, n, ref work);
1392                    }
1393                    else
1394                    {
1395                        reflections.applyreflectionfromtheleft(ref z, taup[i], ref v, i, n, 1, zcolumns, ref work);
1396                    }
1397                    i = i+istep;
1398                }
1399                while( i!=i2+istep );
1400            }
1401        }
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)
1414        {
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];
1432            }
1433            else
1434            {
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            }
1444        }
1445    }
1446}
Note: See TracBrowser for help on using the repository browser.