source: git/Singular/svd/bidiagonal.h @ aad4ca4

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