Free cookie consent management tool by TermsFeed Policy Generator

source: trunk/sources/HeuristicLab.ExtLibs/HeuristicLab.ALGLIB/2.3.0/ALGLIB-2.3.0/bidiagonal.cs @ 2806

Last change on this file since 2806 was 2806, checked in by gkronber, 14 years ago

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