Free cookie consent management tool by TermsFeed Policy Generator

source: trunk/sources/HeuristicLab.ExtLibs/HeuristicLab.ALGLIB/2.1.2.2591/ALGLIB-2.1.2.2591/bidiagonal.cs @ 3031

Last change on this file since 3031 was 2645, checked in by mkommend, 15 years ago

extracted external libraries and adapted dependent plugins (ticket #837)
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            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        public static void tobidiagonal(ref double[,] a,
778            int m,
779            int n,
780            ref double[] tauq,
781            ref double[] taup)
782        {
783            double[] work = new double[0];
784            double[] t = new double[0];
785            int minmn = 0;
786            int maxmn = 0;
787            int i = 0;
788            double ltau = 0;
789            int mmip1 = 0;
790            int nmi = 0;
791            int ip1 = 0;
792            int nmip1 = 0;
793            int mmi = 0;
794            int i_ = 0;
795            int i1_ = 0;
796
797            minmn = Math.Min(m, n);
798            maxmn = Math.Max(m, n);
799            work = new double[maxmn+1];
800            t = new double[maxmn+1];
801            taup = new double[minmn+1];
802            tauq = new double[minmn+1];
803            if( m>=n )
804            {
805               
806                //
807                // Reduce to upper bidiagonal form
808                //
809                for(i=1; i<=n; i++)
810                {
811                   
812                    //
813                    // Generate elementary reflector H(i) to annihilate A(i+1:m,i)
814                    //
815                    mmip1 = m-i+1;
816                    i1_ = (i) - (1);
817                    for(i_=1; i_<=mmip1;i_++)
818                    {
819                        t[i_] = a[i_+i1_,i];
820                    }
821                    reflections.generatereflection(ref t, mmip1, ref ltau);
822                    tauq[i] = ltau;
823                    i1_ = (1) - (i);
824                    for(i_=i; i_<=m;i_++)
825                    {
826                        a[i_,i] = t[i_+i1_];
827                    }
828                    t[1] = 1;
829                   
830                    //
831                    // Apply H(i) to A(i:m,i+1:n) from the left
832                    //
833                    reflections.applyreflectionfromtheleft(ref a, ltau, ref t, i, m, i+1, n, ref work);
834                    if( i<n )
835                    {
836                       
837                        //
838                        // Generate elementary reflector G(i) to annihilate
839                        // A(i,i+2:n)
840                        //
841                        nmi = n-i;
842                        ip1 = i+1;
843                        i1_ = (ip1) - (1);
844                        for(i_=1; i_<=nmi;i_++)
845                        {
846                            t[i_] = a[i,i_+i1_];
847                        }
848                        reflections.generatereflection(ref t, nmi, ref ltau);
849                        taup[i] = ltau;
850                        i1_ = (1) - (ip1);
851                        for(i_=ip1; i_<=n;i_++)
852                        {
853                            a[i,i_] = t[i_+i1_];
854                        }
855                        t[1] = 1;
856                       
857                        //
858                        // Apply G(i) to A(i+1:m,i+1:n) from the right
859                        //
860                        reflections.applyreflectionfromtheright(ref a, ltau, ref t, i+1, m, i+1, n, ref work);
861                    }
862                    else
863                    {
864                        taup[i] = 0;
865                    }
866                }
867            }
868            else
869            {
870               
871                //
872                // Reduce to lower bidiagonal form
873                //
874                for(i=1; i<=m; i++)
875                {
876                   
877                    //
878                    // Generate elementary reflector G(i) to annihilate A(i,i+1:n)
879                    //
880                    nmip1 = n-i+1;
881                    i1_ = (i) - (1);
882                    for(i_=1; i_<=nmip1;i_++)
883                    {
884                        t[i_] = a[i,i_+i1_];
885                    }
886                    reflections.generatereflection(ref t, nmip1, ref ltau);
887                    taup[i] = ltau;
888                    i1_ = (1) - (i);
889                    for(i_=i; i_<=n;i_++)
890                    {
891                        a[i,i_] = t[i_+i1_];
892                    }
893                    t[1] = 1;
894                   
895                    //
896                    // Apply G(i) to A(i+1:m,i:n) from the right
897                    //
898                    reflections.applyreflectionfromtheright(ref a, ltau, ref t, i+1, m, i, n, ref work);
899                    if( i<m )
900                    {
901                       
902                        //
903                        // Generate elementary reflector H(i) to annihilate
904                        // A(i+2:m,i)
905                        //
906                        mmi = m-i;
907                        ip1 = i+1;
908                        i1_ = (ip1) - (1);
909                        for(i_=1; i_<=mmi;i_++)
910                        {
911                            t[i_] = a[i_+i1_,i];
912                        }
913                        reflections.generatereflection(ref t, mmi, ref ltau);
914                        tauq[i] = ltau;
915                        i1_ = (1) - (ip1);
916                        for(i_=ip1; i_<=m;i_++)
917                        {
918                            a[i_,i] = t[i_+i1_];
919                        }
920                        t[1] = 1;
921                       
922                        //
923                        // Apply H(i) to A(i+1:m,i+1:n) from the left
924                        //
925                        reflections.applyreflectionfromtheleft(ref a, ltau, ref t, i+1, m, i+1, n, ref work);
926                    }
927                    else
928                    {
929                        tauq[i] = 0;
930                    }
931                }
932            }
933        }
934
935
936        public static void unpackqfrombidiagonal(ref double[,] qp,
937            int m,
938            int n,
939            ref double[] tauq,
940            int qcolumns,
941            ref double[,] q)
942        {
943            int i = 0;
944            int j = 0;
945            int ip1 = 0;
946            double[] v = new double[0];
947            double[] work = new double[0];
948            int vm = 0;
949            int i_ = 0;
950            int i1_ = 0;
951
952            System.Diagnostics.Debug.Assert(qcolumns<=m, "UnpackQFromBidiagonal: QColumns>M!");
953            if( m==0 | n==0 | qcolumns==0 )
954            {
955                return;
956            }
957           
958            //
959            // init
960            //
961            q = new double[m+1, qcolumns+1];
962            v = new double[m+1];
963            work = new double[qcolumns+1];
964           
965            //
966            // prepare Q
967            //
968            for(i=1; i<=m; i++)
969            {
970                for(j=1; j<=qcolumns; j++)
971                {
972                    if( i==j )
973                    {
974                        q[i,j] = 1;
975                    }
976                    else
977                    {
978                        q[i,j] = 0;
979                    }
980                }
981            }
982            if( m>=n )
983            {
984                for(i=Math.Min(n, qcolumns); i>=1; i--)
985                {
986                    vm = m-i+1;
987                    i1_ = (i) - (1);
988                    for(i_=1; i_<=vm;i_++)
989                    {
990                        v[i_] = qp[i_+i1_,i];
991                    }
992                    v[1] = 1;
993                    reflections.applyreflectionfromtheleft(ref q, tauq[i], ref v, i, m, 1, qcolumns, ref work);
994                }
995            }
996            else
997            {
998                for(i=Math.Min(m-1, qcolumns-1); i>=1; i--)
999                {
1000                    vm = m-i;
1001                    ip1 = i+1;
1002                    i1_ = (ip1) - (1);
1003                    for(i_=1; i_<=vm;i_++)
1004                    {
1005                        v[i_] = qp[i_+i1_,i];
1006                    }
1007                    v[1] = 1;
1008                    reflections.applyreflectionfromtheleft(ref q, tauq[i], ref v, i+1, m, 1, qcolumns, ref work);
1009                }
1010            }
1011        }
1012
1013
1014        public static void multiplybyqfrombidiagonal(ref double[,] qp,
1015            int m,
1016            int n,
1017            ref double[] tauq,
1018            ref double[,] z,
1019            int zrows,
1020            int zcolumns,
1021            bool fromtheright,
1022            bool dotranspose)
1023        {
1024            int i = 0;
1025            int ip1 = 0;
1026            int i1 = 0;
1027            int i2 = 0;
1028            int istep = 0;
1029            double[] v = new double[0];
1030            double[] work = new double[0];
1031            int vm = 0;
1032            int mx = 0;
1033            int i_ = 0;
1034            int i1_ = 0;
1035
1036            if( m<=0 | n<=0 | zrows<=0 | zcolumns<=0 )
1037            {
1038                return;
1039            }
1040            System.Diagnostics.Debug.Assert(fromtheright & zcolumns==m | !fromtheright & zrows==m, "MultiplyByQFromBidiagonal: incorrect Z size!");
1041           
1042            //
1043            // init
1044            //
1045            mx = Math.Max(m, n);
1046            mx = Math.Max(mx, zrows);
1047            mx = Math.Max(mx, zcolumns);
1048            v = new double[mx+1];
1049            work = new double[mx+1];
1050            if( m>=n )
1051            {
1052               
1053                //
1054                // setup
1055                //
1056                if( fromtheright )
1057                {
1058                    i1 = 1;
1059                    i2 = n;
1060                    istep = +1;
1061                }
1062                else
1063                {
1064                    i1 = n;
1065                    i2 = 1;
1066                    istep = -1;
1067                }
1068                if( dotranspose )
1069                {
1070                    i = i1;
1071                    i1 = i2;
1072                    i2 = i;
1073                    istep = -istep;
1074                }
1075               
1076                //
1077                // Process
1078                //
1079                i = i1;
1080                do
1081                {
1082                    vm = m-i+1;
1083                    i1_ = (i) - (1);
1084                    for(i_=1; i_<=vm;i_++)
1085                    {
1086                        v[i_] = qp[i_+i1_,i];
1087                    }
1088                    v[1] = 1;
1089                    if( fromtheright )
1090                    {
1091                        reflections.applyreflectionfromtheright(ref z, tauq[i], ref v, 1, zrows, i, m, ref work);
1092                    }
1093                    else
1094                    {
1095                        reflections.applyreflectionfromtheleft(ref z, tauq[i], ref v, i, m, 1, zcolumns, ref work);
1096                    }
1097                    i = i+istep;
1098                }
1099                while( i!=i2+istep );
1100            }
1101            else
1102            {
1103               
1104                //
1105                // setup
1106                //
1107                if( fromtheright )
1108                {
1109                    i1 = 1;
1110                    i2 = m-1;
1111                    istep = +1;
1112                }
1113                else
1114                {
1115                    i1 = m-1;
1116                    i2 = 1;
1117                    istep = -1;
1118                }
1119                if( dotranspose )
1120                {
1121                    i = i1;
1122                    i1 = i2;
1123                    i2 = i;
1124                    istep = -istep;
1125                }
1126               
1127                //
1128                // Process
1129                //
1130                if( m-1>0 )
1131                {
1132                    i = i1;
1133                    do
1134                    {
1135                        vm = m-i;
1136                        ip1 = i+1;
1137                        i1_ = (ip1) - (1);
1138                        for(i_=1; i_<=vm;i_++)
1139                        {
1140                            v[i_] = qp[i_+i1_,i];
1141                        }
1142                        v[1] = 1;
1143                        if( fromtheright )
1144                        {
1145                            reflections.applyreflectionfromtheright(ref z, tauq[i], ref v, 1, zrows, i+1, m, ref work);
1146                        }
1147                        else
1148                        {
1149                            reflections.applyreflectionfromtheleft(ref z, tauq[i], ref v, i+1, m, 1, zcolumns, ref work);
1150                        }
1151                        i = i+istep;
1152                    }
1153                    while( i!=i2+istep );
1154                }
1155            }
1156        }
1157
1158
1159        public static void unpackptfrombidiagonal(ref double[,] qp,
1160            int m,
1161            int n,
1162            ref double[] taup,
1163            int ptrows,
1164            ref double[,] pt)
1165        {
1166            int i = 0;
1167            int j = 0;
1168            int ip1 = 0;
1169            double[] v = new double[0];
1170            double[] work = new double[0];
1171            int vm = 0;
1172            int i_ = 0;
1173            int i1_ = 0;
1174
1175            System.Diagnostics.Debug.Assert(ptrows<=n, "UnpackPTFromBidiagonal: PTRows>N!");
1176            if( m==0 | n==0 | ptrows==0 )
1177            {
1178                return;
1179            }
1180           
1181            //
1182            // init
1183            //
1184            pt = new double[ptrows+1, n+1];
1185            v = new double[n+1];
1186            work = new double[ptrows+1];
1187           
1188            //
1189            // prepare PT
1190            //
1191            for(i=1; i<=ptrows; i++)
1192            {
1193                for(j=1; j<=n; j++)
1194                {
1195                    if( i==j )
1196                    {
1197                        pt[i,j] = 1;
1198                    }
1199                    else
1200                    {
1201                        pt[i,j] = 0;
1202                    }
1203                }
1204            }
1205            if( m>=n )
1206            {
1207                for(i=Math.Min(n-1, ptrows-1); i>=1; i--)
1208                {
1209                    vm = n-i;
1210                    ip1 = i+1;
1211                    i1_ = (ip1) - (1);
1212                    for(i_=1; i_<=vm;i_++)
1213                    {
1214                        v[i_] = qp[i,i_+i1_];
1215                    }
1216                    v[1] = 1;
1217                    reflections.applyreflectionfromtheright(ref pt, taup[i], ref v, 1, ptrows, i+1, n, ref work);
1218                }
1219            }
1220            else
1221            {
1222                for(i=Math.Min(m, ptrows); i>=1; i--)
1223                {
1224                    vm = n-i+1;
1225                    i1_ = (i) - (1);
1226                    for(i_=1; i_<=vm;i_++)
1227                    {
1228                        v[i_] = qp[i,i_+i1_];
1229                    }
1230                    v[1] = 1;
1231                    reflections.applyreflectionfromtheright(ref pt, taup[i], ref v, 1, ptrows, i, n, ref work);
1232                }
1233            }
1234        }
1235
1236
1237        public static void multiplybypfrombidiagonal(ref double[,] qp,
1238            int m,
1239            int n,
1240            ref double[] taup,
1241            ref double[,] z,
1242            int zrows,
1243            int zcolumns,
1244            bool fromtheright,
1245            bool dotranspose)
1246        {
1247            int i = 0;
1248            int ip1 = 0;
1249            double[] v = new double[0];
1250            double[] work = new double[0];
1251            int vm = 0;
1252            int mx = 0;
1253            int i1 = 0;
1254            int i2 = 0;
1255            int istep = 0;
1256            int i_ = 0;
1257            int i1_ = 0;
1258
1259            if( m<=0 | n<=0 | zrows<=0 | zcolumns<=0 )
1260            {
1261                return;
1262            }
1263            System.Diagnostics.Debug.Assert(fromtheright & zcolumns==n | !fromtheright & zrows==n, "MultiplyByQFromBidiagonal: incorrect Z size!");
1264           
1265            //
1266            // init
1267            //
1268            mx = Math.Max(m, n);
1269            mx = Math.Max(mx, zrows);
1270            mx = Math.Max(mx, zcolumns);
1271            v = new double[mx+1];
1272            work = new double[mx+1];
1273            v = new double[mx+1];
1274            work = new double[mx+1];
1275            if( m>=n )
1276            {
1277               
1278                //
1279                // setup
1280                //
1281                if( fromtheright )
1282                {
1283                    i1 = n-1;
1284                    i2 = 1;
1285                    istep = -1;
1286                }
1287                else
1288                {
1289                    i1 = 1;
1290                    i2 = n-1;
1291                    istep = +1;
1292                }
1293                if( !dotranspose )
1294                {
1295                    i = i1;
1296                    i1 = i2;
1297                    i2 = i;
1298                    istep = -istep;
1299                }
1300               
1301                //
1302                // Process
1303                //
1304                if( n-1>0 )
1305                {
1306                    i = i1;
1307                    do
1308                    {
1309                        vm = n-i;
1310                        ip1 = i+1;
1311                        i1_ = (ip1) - (1);
1312                        for(i_=1; i_<=vm;i_++)
1313                        {
1314                            v[i_] = qp[i,i_+i1_];
1315                        }
1316                        v[1] = 1;
1317                        if( fromtheright )
1318                        {
1319                            reflections.applyreflectionfromtheright(ref z, taup[i], ref v, 1, zrows, i+1, n, ref work);
1320                        }
1321                        else
1322                        {
1323                            reflections.applyreflectionfromtheleft(ref z, taup[i], ref v, i+1, n, 1, zcolumns, ref work);
1324                        }
1325                        i = i+istep;
1326                    }
1327                    while( i!=i2+istep );
1328                }
1329            }
1330            else
1331            {
1332               
1333                //
1334                // setup
1335                //
1336                if( fromtheright )
1337                {
1338                    i1 = m;
1339                    i2 = 1;
1340                    istep = -1;
1341                }
1342                else
1343                {
1344                    i1 = 1;
1345                    i2 = m;
1346                    istep = +1;
1347                }
1348                if( !dotranspose )
1349                {
1350                    i = i1;
1351                    i1 = i2;
1352                    i2 = i;
1353                    istep = -istep;
1354                }
1355               
1356                //
1357                // Process
1358                //
1359                i = i1;
1360                do
1361                {
1362                    vm = n-i+1;
1363                    i1_ = (i) - (1);
1364                    for(i_=1; i_<=vm;i_++)
1365                    {
1366                        v[i_] = qp[i,i_+i1_];
1367                    }
1368                    v[1] = 1;
1369                    if( fromtheright )
1370                    {
1371                        reflections.applyreflectionfromtheright(ref z, taup[i], ref v, 1, zrows, i, n, ref work);
1372                    }
1373                    else
1374                    {
1375                        reflections.applyreflectionfromtheleft(ref z, taup[i], ref v, i, n, 1, zcolumns, ref work);
1376                    }
1377                    i = i+istep;
1378                }
1379                while( i!=i2+istep );
1380            }
1381        }
1382
1383
1384        public static void unpackdiagonalsfrombidiagonal(ref double[,] b,
1385            int m,
1386            int n,
1387            ref bool isupper,
1388            ref double[] d,
1389            ref double[] e)
1390        {
1391            int i = 0;
1392
1393            isupper = m>=n;
1394            if( m==0 | n==0 )
1395            {
1396                return;
1397            }
1398            if( isupper )
1399            {
1400                d = new double[n+1];
1401                e = new double[n+1];
1402                for(i=1; i<=n-1; i++)
1403                {
1404                    d[i] = b[i,i];
1405                    e[i] = b[i,i+1];
1406                }
1407                d[n] = b[n,n];
1408            }
1409            else
1410            {
1411                d = new double[m+1];
1412                e = new double[m+1];
1413                for(i=1; i<=m-1; i++)
1414                {
1415                    d[i] = b[i,i];
1416                    e[i] = b[i+1,i];
1417                }
1418                d[m] = b[m,m];
1419            }
1420        }
1421    }
1422}
Note: See TracBrowser for help on using the repository browser.