My Project
Loading...
Searching...
No Matches
bidiagonal.h
Go to the documentation of this file.
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
m
int m
Definition: cfEzgcd.cc:128
i
int i
Definition: cfEzgcd.cc:132
b
CanonicalForm b
Definition: cfModGcd.cc:4103
amp::ampf
Definition: amp.h:82
ap::ap_error::make_assertion
static void make_assertion(bool bClause)
Definition: ap.h:49
ap::template_1d_array
Definition: ap.h:657
ap::template_1d_array::getvector
raw_vector< T > getvector(int iStart, int iEnd)
Definition: ap.h:776
ap::template_1d_array::setbounds
void setbounds(int iLow, int iHigh)
Definition: ap.h:735
ap::template_2d_array
Definition: ap.h:807
v
const Variable & v
< [in] a sqrfree bivariate poly
Definition: facBivar.h:39
j
int j
Definition: facHensel.cc:110
ap::maxint
int maxint(int m1, int m2)
Definition: ap.cpp:162
ap::vmove
void vmove(raw_vector< T > vdst, const_raw_vector< T > vsrc)
Definition: ap.h:237
ap::minint
int minint(int m1, int m2)
Definition: ap.cpp:167
bidiagonal::unpackptfrombidiagonal
void unpackptfrombidiagonal(const ap::template_2d_array< amp::ampf< Precision > > &qp, int m, int n, const ap::template_1d_array< amp::ampf< Precision > > &taup, int ptrows, ap::template_2d_array< amp::ampf< Precision > > &pt)
Definition: bidiagonal.h:1196
bidiagonal::multiplybyqfrombidiagonal
void multiplybyqfrombidiagonal(const ap::template_2d_array< amp::ampf< Precision > > &qp, int m, int n, const ap::template_1d_array< amp::ampf< Precision > > &tauq, ap::template_2d_array< amp::ampf< Precision > > &z, int zrows, int zcolumns, bool fromtheright, bool dotranspose)
Definition: bidiagonal.h:1055
bidiagonal::rmatrixbdunpackpt
void rmatrixbdunpackpt(const ap::template_2d_array< amp::ampf< Precision > > &qp, int m, int n, const ap::template_1d_array< amp::ampf< Precision > > &taup, int ptrows, ap::template_2d_array< amp::ampf< Precision > > &pt)
Definition: bidiagonal.h:575
bidiagonal::unpackqfrombidiagonal
void unpackqfrombidiagonal(const ap::template_2d_array< amp::ampf< Precision > > &qp, int m, int n, const ap::template_1d_array< amp::ampf< Precision > > &tauq, int qcolumns, ap::template_2d_array< amp::ampf< Precision > > &q)
Definition: bidiagonal.h:981
bidiagonal::rmatrixbdunpackq
void rmatrixbdunpackq(const ap::template_2d_array< amp::ampf< Precision > > &qp, int m, int n, const ap::template_1d_array< amp::ampf< Precision > > &tauq, int qcolumns, ap::template_2d_array< amp::ampf< Precision > > &q)
Definition: bidiagonal.h:348
bidiagonal::rmatrixbdmultiplybyq
void rmatrixbdmultiplybyq(const ap::template_2d_array< amp::ampf< Precision > > &qp, int m, int n, const ap::template_1d_array< amp::ampf< Precision > > &tauq, ap::template_2d_array< amp::ampf< Precision > > &z, int zrows, int zcolumns, bool fromtheright, bool dotranspose)
Definition: bidiagonal.h:422
bidiagonal::unpackdiagonalsfrombidiagonal
void unpackdiagonalsfrombidiagonal(const ap::template_2d_array< amp::ampf< Precision > > &b, int m, int n, bool &isupper, ap::template_1d_array< amp::ampf< Precision > > &d, ap::template_1d_array< amp::ampf< Precision > > &e)
Definition: bidiagonal.h:1413
bidiagonal::rmatrixbdunpackdiagonals
void rmatrixbdunpackdiagonals(const ap::template_2d_array< amp::ampf< Precision > > &b, int m, int n, bool &isupper, ap::template_1d_array< amp::ampf< Precision > > &d, ap::template_1d_array< amp::ampf< Precision > > &e)
Definition: bidiagonal.h:805
bidiagonal::tobidiagonal
void tobidiagonal(ap::template_2d_array< amp::ampf< Precision > > &a, int m, int n, ap::template_1d_array< amp::ampf< Precision > > &tauq, ap::template_1d_array< amp::ampf< Precision > > &taup)
Definition: bidiagonal.h:850
bidiagonal::rmatrixbd
void rmatrixbd(ap::template_2d_array< amp::ampf< Precision > > &a, int m, int n, ap::template_1d_array< amp::ampf< Precision > > &tauq, ap::template_1d_array< amp::ampf< Precision > > &taup)
Definition: bidiagonal.h:194
bidiagonal::rmatrixbdmultiplybyp
void rmatrixbdmultiplybyp(const ap::template_2d_array< amp::ampf< Precision > > &qp, int m, int n, const ap::template_1d_array< amp::ampf< Precision > > &taup, ap::template_2d_array< amp::ampf< Precision > > &z, int zrows, int zcolumns, bool fromtheright, bool dotranspose)
Definition: bidiagonal.h:649
bidiagonal::multiplybypfrombidiagonal
void multiplybypfrombidiagonal(const ap::template_2d_array< amp::ampf< Precision > > &qp, int m, int n, const ap::template_1d_array< amp::ampf< Precision > > &taup, ap::template_2d_array< amp::ampf< Precision > > &z, int zrows, int zcolumns, bool fromtheright, bool dotranspose)
Definition: bidiagonal.h:1270
Generated on Mon Feb 27 2023 10:53:50 for My Project by   doxygen 1.9.5 for Singular