My Project
svd.h
Go to the documentation of this file.
1/*************************************************************************
2Copyright (c) 2005-2007, Sergey Bochkanov (ALGLIB project).
3
4Redistribution and use in source and binary forms, with or without
5modification, are permitted provided that the following conditions are
6met:
7
8- Redistributions of source code must retain the above copyright
9 notice, this list of conditions and the following disclaimer.
10
11- Redistributions in binary form must reproduce the above copyright
12 notice, this list of conditions and the following disclaimer listed
13 in this license in the documentation and/or other materials
14 provided with the distribution.
15
16- Neither the name of the copyright holders nor the names of its
17 contributors may be used to endorse or promote products derived from
18 this software without specific prior written permission.
19
20THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
21"AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
22LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
23A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT
24OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
25SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
26LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
27DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
28THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
29(INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
30OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
31*************************************************************************/
32
33#ifndef _svd_h
34#define _svd_h
35
36/*MAKEHEADER*/
37#include "ap.h"
38
39/*MAKEHEADER*/
40#include "amp.h"
41
42/*MAKEHEADER*/
43#include "reflections.h"
44
45/*MAKEHEADER*/
46#include "bidiagonal.h"
47
48/*MAKEHEADER*/
49#include "qr.h"
50
51/*MAKEHEADER*/
52#include "lq.h"
53
54/*MAKEHEADER*/
55#include "blas.h"
56
57/*MAKEHEADER*/
58#include "rotations.h"
59
60/*MAKEHEADER*/
61#include "bdsvd.h"
62
63namespace svd
64{
65 template<unsigned int Precision>
66 bool rmatrixsvd(ap::template_2d_array< amp::ampf<Precision> > a,
67 int m,
68 int n,
69 int uneeded,
70 int vtneeded,
71 int additionalmemory,
72 ap::template_1d_array< amp::ampf<Precision> >& w,
73 ap::template_2d_array< amp::ampf<Precision> >& u,
74 ap::template_2d_array< amp::ampf<Precision> >& vt);
75 template<unsigned int Precision>
76 bool svddecomposition(ap::template_2d_array< amp::ampf<Precision> > a,
77 int m,
78 int n,
79 int uneeded,
80 int vtneeded,
81 int additionalmemory,
82 ap::template_1d_array< amp::ampf<Precision> >& w,
83 ap::template_2d_array< amp::ampf<Precision> >& u,
84 ap::template_2d_array< amp::ampf<Precision> >& vt);
85
86
87 /*************************************************************************
88 Singular value decomposition of a rectangular matrix.
89
90 The algorithm calculates the singular value decomposition of a matrix of
91 size MxN: A = U * S * V^T
92
93 The algorithm finds the singular values and, optionally, matrices U and V^T.
94 The algorithm can find both first min(M,N) columns of matrix U and rows of
95 matrix V^T (singular vectors), and matrices U and V^T wholly (of sizes MxM
96 and NxN respectively).
97
98 Take into account that the subroutine does not return matrix V but V^T.
99
100 Input parameters:
101 A - matrix to be decomposed.
102 Array whose indexes range within [0..M-1, 0..N-1].
103 M - number of rows in matrix A.
104 N - number of columns in matrix A.
105 UNeeded - 0, 1 or 2. See the description of the parameter U.
106 VTNeeded - 0, 1 or 2. See the description of the parameter VT.
107 AdditionalMemory -
108 If the parameter:
109 * equals 0, the algorithm doesn't use additional
110 memory (lower requirements, lower performance).
111 * equals 1, the algorithm uses additional
112 memory of size min(M,N)*min(M,N) of real numbers.
113 It often speeds up the algorithm.
114 * equals 2, the algorithm uses additional
115 memory of size M*min(M,N) of real numbers.
116 It allows to get a maximum performance.
117 The recommended value of the parameter is 2.
118
119 Output parameters:
120 W - contains singular values in descending order.
121 U - if UNeeded=0, U isn't changed, the left singular vectors
122 are not calculated.
123 if Uneeded=1, U contains left singular vectors (first
124 min(M,N) columns of matrix U). Array whose indexes range
125 within [0..M-1, 0..Min(M,N)-1].
126 if UNeeded=2, U contains matrix U wholly. Array whose
127 indexes range within [0..M-1, 0..M-1].
128 VT - if VTNeeded=0, VT isn't changed, the right singular vectors
129 are not calculated.
130 if VTNeeded=1, VT contains right singular vectors (first
131 min(M,N) rows of matrix V^T). Array whose indexes range
132 within [0..min(M,N)-1, 0..N-1].
133 if VTNeeded=2, VT contains matrix V^T wholly. Array whose
134 indexes range within [0..N-1, 0..N-1].
135
136 -- ALGLIB --
137 Copyright 2005 by Bochkanov Sergey
138 *************************************************************************/
139 template<unsigned int Precision>
140 bool rmatrixsvd(ap::template_2d_array< amp::ampf<Precision> > a,
141 int m,
142 int n,
143 int uneeded,
144 int vtneeded,
145 int additionalmemory,
146 ap::template_1d_array< amp::ampf<Precision> >& w,
147 ap::template_2d_array< amp::ampf<Precision> >& u,
148 ap::template_2d_array< amp::ampf<Precision> >& vt)
149 {
150 bool result;
151 ap::template_1d_array< amp::ampf<Precision> > tauq;
152 ap::template_1d_array< amp::ampf<Precision> > taup;
153 ap::template_1d_array< amp::ampf<Precision> > tau;
154 ap::template_1d_array< amp::ampf<Precision> > e;
155 ap::template_1d_array< amp::ampf<Precision> > work;
156 ap::template_2d_array< amp::ampf<Precision> > t2;
157 bool isupper;
158 int minmn;
159 int ncu;
160 int nrvt;
161 int nru;
162 int ncvt;
163 int i;
164 int j;
165 int im1;
166 amp::ampf<Precision> sm;
167
168
169 result = true;
170 if( m==0 || n==0 )
171 {
172 return result;
173 }
174 ap::ap_error::make_assertion(uneeded>=0 && uneeded<=2);
175 ap::ap_error::make_assertion(vtneeded>=0 && vtneeded<=2);
176 ap::ap_error::make_assertion(additionalmemory>=0 && additionalmemory<=2);
177
178 //
179 // initialize
180 //
181 minmn = ap::minint(m, n);
182 w.setbounds(1, minmn);
183 ncu = 0;
184 nru = 0;
185 if( uneeded==1 )
186 {
187 nru = m;
188 ncu = minmn;
189 u.setbounds(0, nru-1, 0, ncu-1);
190 }
191 if( uneeded==2 )
192 {
193 nru = m;
194 ncu = m;
195 u.setbounds(0, nru-1, 0, ncu-1);
196 }
197 nrvt = 0;
198 ncvt = 0;
199 if( vtneeded==1 )
200 {
201 nrvt = minmn;
202 ncvt = n;
203 vt.setbounds(0, nrvt-1, 0, ncvt-1);
204 }
205 if( vtneeded==2 )
206 {
207 nrvt = n;
208 ncvt = n;
209 vt.setbounds(0, nrvt-1, 0, ncvt-1);
210 }
211
212 //
213 // M much larger than N
214 // Use bidiagonal reduction with QR-decomposition
215 //
216 if( m>amp::ampf<Precision>("1.6")*n )
217 {
218 if( uneeded==0 )
219 {
220
221 //
222 // No left singular vectors to be computed
223 //
224 qr::rmatrixqr<Precision>(a, m, n, tau);
225 for(i=0; i<=n-1; i++)
226 {
227 for(j=0; j<=i-1; j++)
228 {
229 a(i,j) = 0;
230 }
231 }
232 bidiagonal::rmatrixbd<Precision>(a, n, n, tauq, taup);
233 bidiagonal::rmatrixbdunpackpt<Precision>(a, n, n, taup, nrvt, vt);
234 bidiagonal::rmatrixbdunpackdiagonals<Precision>(a, n, n, isupper, w, e);
235 result = bdsvd::rmatrixbdsvd<Precision>(w, e, n, isupper, false, u, 0, a, 0, vt, ncvt);
236 return result;
237 }
238 else
239 {
240
241 //
242 // Left singular vectors (may be full matrix U) to be computed
243 //
244 qr::rmatrixqr<Precision>(a, m, n, tau);
245 qr::rmatrixqrunpackq<Precision>(a, m, n, tau, ncu, u);
246 for(i=0; i<=n-1; i++)
247 {
248 for(j=0; j<=i-1; j++)
249 {
250 a(i,j) = 0;
251 }
252 }
253 bidiagonal::rmatrixbd<Precision>(a, n, n, tauq, taup);
254 bidiagonal::rmatrixbdunpackpt<Precision>(a, n, n, taup, nrvt, vt);
255 bidiagonal::rmatrixbdunpackdiagonals<Precision>(a, n, n, isupper, w, e);
256 if( additionalmemory<1 )
257 {
258
259 //
260 // No additional memory can be used
261 //
262 bidiagonal::rmatrixbdmultiplybyq<Precision>(a, n, n, tauq, u, m, n, true, false);
263 result = bdsvd::rmatrixbdsvd<Precision>(w, e, n, isupper, false, u, m, a, 0, vt, ncvt);
264 }
265 else
266 {
267
268 //
269 // Large U. Transforming intermediate matrix T2
270 //
271 work.setbounds(1, ap::maxint(m, n));
272 bidiagonal::rmatrixbdunpackq<Precision>(a, n, n, tauq, n, t2);
273 blas::copymatrix<Precision>(u, 0, m-1, 0, n-1, a, 0, m-1, 0, n-1);
274 blas::inplacetranspose<Precision>(t2, 0, n-1, 0, n-1, work);
275 result = bdsvd::rmatrixbdsvd<Precision>(w, e, n, isupper, false, u, 0, t2, n, vt, ncvt);
276 blas::matrixmatrixmultiply<Precision>(a, 0, m-1, 0, n-1, false, t2, 0, n-1, 0, n-1, true, amp::ampf<Precision>("1.0"), u, 0, m-1, 0, n-1, amp::ampf<Precision>("0.0"), work);
277 }
278 return result;
279 }
280 }
281
282 //
283 // N much larger than M
284 // Use bidiagonal reduction with LQ-decomposition
285 //
286 if( n>amp::ampf<Precision>("1.6")*m )
287 {
288 if( vtneeded==0 )
289 {
290
291 //
292 // No right singular vectors to be computed
293 //
294 lq::rmatrixlq<Precision>(a, m, n, tau);
295 for(i=0; i<=m-1; i++)
296 {
297 for(j=i+1; j<=m-1; j++)
298 {
299 a(i,j) = 0;
300 }
301 }
302 bidiagonal::rmatrixbd<Precision>(a, m, m, tauq, taup);
303 bidiagonal::rmatrixbdunpackq<Precision>(a, m, m, tauq, ncu, u);
304 bidiagonal::rmatrixbdunpackdiagonals<Precision>(a, m, m, isupper, w, e);
305 work.setbounds(1, m);
306 blas::inplacetranspose<Precision>(u, 0, nru-1, 0, ncu-1, work);
307 result = bdsvd::rmatrixbdsvd<Precision>(w, e, m, isupper, false, a, 0, u, nru, vt, 0);
308 blas::inplacetranspose<Precision>(u, 0, nru-1, 0, ncu-1, work);
309 return result;
310 }
311 else
312 {
313
314 //
315 // Right singular vectors (may be full matrix VT) to be computed
316 //
317 lq::rmatrixlq<Precision>(a, m, n, tau);
318 lq::rmatrixlqunpackq<Precision>(a, m, n, tau, nrvt, vt);
319 for(i=0; i<=m-1; i++)
320 {
321 for(j=i+1; j<=m-1; j++)
322 {
323 a(i,j) = 0;
324 }
325 }
326 bidiagonal::rmatrixbd<Precision>(a, m, m, tauq, taup);
327 bidiagonal::rmatrixbdunpackq<Precision>(a, m, m, tauq, ncu, u);
328 bidiagonal::rmatrixbdunpackdiagonals<Precision>(a, m, m, isupper, w, e);
329 work.setbounds(1, ap::maxint(m, n));
330 blas::inplacetranspose<Precision>(u, 0, nru-1, 0, ncu-1, work);
331 if( additionalmemory<1 )
332 {
333
334 //
335 // No additional memory available
336 //
337 bidiagonal::rmatrixbdmultiplybyp<Precision>(a, m, m, taup, vt, m, n, false, true);
338 result = bdsvd::rmatrixbdsvd<Precision>(w, e, m, isupper, false, a, 0, u, nru, vt, n);
339 }
340 else
341 {
342
343 //
344 // Large VT. Transforming intermediate matrix T2
345 //
346 bidiagonal::rmatrixbdunpackpt<Precision>(a, m, m, taup, m, t2);
347 result = bdsvd::rmatrixbdsvd<Precision>(w, e, m, isupper, false, a, 0, u, nru, t2, m);
348 blas::copymatrix<Precision>(vt, 0, m-1, 0, n-1, a, 0, m-1, 0, n-1);
349 blas::matrixmatrixmultiply<Precision>(t2, 0, m-1, 0, m-1, false, a, 0, m-1, 0, n-1, false, amp::ampf<Precision>("1.0"), vt, 0, m-1, 0, n-1, amp::ampf<Precision>("0.0"), work);
350 }
351 blas::inplacetranspose<Precision>(u, 0, nru-1, 0, ncu-1, work);
352 return result;
353 }
354 }
355
356 //
357 // M<=N
358 // We can use inplace transposition of U to get rid of columnwise operations
359 //
360 if( m<=n )
361 {
362 bidiagonal::rmatrixbd<Precision>(a, m, n, tauq, taup);
363 bidiagonal::rmatrixbdunpackq<Precision>(a, m, n, tauq, ncu, u);
364 bidiagonal::rmatrixbdunpackpt<Precision>(a, m, n, taup, nrvt, vt);
365 bidiagonal::rmatrixbdunpackdiagonals<Precision>(a, m, n, isupper, w, e);
366 work.setbounds(1, m);
367 blas::inplacetranspose<Precision>(u, 0, nru-1, 0, ncu-1, work);
368 result = bdsvd::rmatrixbdsvd<Precision>(w, e, minmn, isupper, false, a, 0, u, nru, vt, ncvt);
369 blas::inplacetranspose<Precision>(u, 0, nru-1, 0, ncu-1, work);
370 return result;
371 }
372
373 //
374 // Simple bidiagonal reduction
375 //
376 bidiagonal::rmatrixbd<Precision>(a, m, n, tauq, taup);
377 bidiagonal::rmatrixbdunpackq<Precision>(a, m, n, tauq, ncu, u);
378 bidiagonal::rmatrixbdunpackpt<Precision>(a, m, n, taup, nrvt, vt);
379 bidiagonal::rmatrixbdunpackdiagonals<Precision>(a, m, n, isupper, w, e);
380 if( additionalmemory<2 || uneeded==0 )
381 {
382
383 //
384 // We cant use additional memory or there is no need in such operations
385 //
386 result = bdsvd::rmatrixbdsvd<Precision>(w, e, minmn, isupper, false, u, nru, a, 0, vt, ncvt);
387 }
388 else
389 {
390
391 //
392 // We can use additional memory
393 //
394 t2.setbounds(0, minmn-1, 0, m-1);
395 blas::copyandtranspose<Precision>(u, 0, m-1, 0, minmn-1, t2, 0, minmn-1, 0, m-1);
396 result = bdsvd::rmatrixbdsvd<Precision>(w, e, minmn, isupper, false, u, 0, t2, m, vt, ncvt);
397 blas::copyandtranspose<Precision>(t2, 0, minmn-1, 0, m-1, u, 0, m-1, 0, minmn-1);
398 }
399 return result;
400 }
401
402
403 /*************************************************************************
404 Obsolete 1-based subroutine.
405 See RMatrixSVD for 0-based replacement.
406 *************************************************************************/
407 template<unsigned int Precision>
408 bool svddecomposition(ap::template_2d_array< amp::ampf<Precision> > a,
409 int m,
410 int n,
411 int uneeded,
412 int vtneeded,
413 int additionalmemory,
414 ap::template_1d_array< amp::ampf<Precision> >& w,
415 ap::template_2d_array< amp::ampf<Precision> >& u,
416 ap::template_2d_array< amp::ampf<Precision> >& vt)
417 {
418 bool result;
419 ap::template_1d_array< amp::ampf<Precision> > tauq;
420 ap::template_1d_array< amp::ampf<Precision> > taup;
421 ap::template_1d_array< amp::ampf<Precision> > tau;
422 ap::template_1d_array< amp::ampf<Precision> > e;
423 ap::template_1d_array< amp::ampf<Precision> > work;
424 ap::template_2d_array< amp::ampf<Precision> > t2;
425 bool isupper;
426 int minmn;
427 int ncu;
428 int nrvt;
429 int nru;
430 int ncvt;
431 int i;
432 int j;
433 int im1;
434 amp::ampf<Precision> sm;
435
436
437 result = true;
438 if( m==0 || n==0 )
439 {
440 return result;
441 }
442 ap::ap_error::make_assertion(uneeded>=0 && uneeded<=2);
443 ap::ap_error::make_assertion(vtneeded>=0 && vtneeded<=2);
444 ap::ap_error::make_assertion(additionalmemory>=0 && additionalmemory<=2);
445
446 //
447 // initialize
448 //
449 minmn = ap::minint(m, n);
450 w.setbounds(1, minmn);
451 ncu = 0;
452 nru = 0;
453 if( uneeded==1 )
454 {
455 nru = m;
456 ncu = minmn;
457 u.setbounds(1, nru, 1, ncu);
458 }
459 if( uneeded==2 )
460 {
461 nru = m;
462 ncu = m;
463 u.setbounds(1, nru, 1, ncu);
464 }
465 nrvt = 0;
466 ncvt = 0;
467 if( vtneeded==1 )
468 {
469 nrvt = minmn;
470 ncvt = n;
471 vt.setbounds(1, nrvt, 1, ncvt);
472 }
473 if( vtneeded==2 )
474 {
475 nrvt = n;
476 ncvt = n;
477 vt.setbounds(1, nrvt, 1, ncvt);
478 }
479
480 //
481 // M much larger than N
482 // Use bidiagonal reduction with QR-decomposition
483 //
484 if( m>amp::ampf<Precision>("1.6")*n )
485 {
486 if( uneeded==0 )
487 {
488
489 //
490 // No left singular vectors to be computed
491 //
492 qr::qrdecomposition<Precision>(a, m, n, tau);
493 for(i=2; i<=n; i++)
494 {
495 for(j=1; j<=i-1; j++)
496 {
497 a(i,j) = 0;
498 }
499 }
500 bidiagonal::tobidiagonal<Precision>(a, n, n, tauq, taup);
501 bidiagonal::unpackptfrombidiagonal<Precision>(a, n, n, taup, nrvt, vt);
502 bidiagonal::unpackdiagonalsfrombidiagonal<Precision>(a, n, n, isupper, w, e);
503 result = bdsvd::bidiagonalsvddecomposition<Precision>(w, e, n, isupper, false, u, 0, a, 0, vt, ncvt);
504 return result;
505 }
506 else
507 {
508
509 //
510 // Left singular vectors (may be full matrix U) to be computed
511 //
512 qr::qrdecomposition<Precision>(a, m, n, tau);
513 qr::unpackqfromqr<Precision>(a, m, n, tau, ncu, u);
514 for(i=2; i<=n; i++)
515 {
516 for(j=1; j<=i-1; j++)
517 {
518 a(i,j) = 0;
519 }
520 }
521 bidiagonal::tobidiagonal<Precision>(a, n, n, tauq, taup);
522 bidiagonal::unpackptfrombidiagonal<Precision>(a, n, n, taup, nrvt, vt);
523 bidiagonal::unpackdiagonalsfrombidiagonal<Precision>(a, n, n, isupper, w, e);
524 if( additionalmemory<1 )
525 {
526
527 //
528 // No additional memory can be used
529 //
530 bidiagonal::multiplybyqfrombidiagonal<Precision>(a, n, n, tauq, u, m, n, true, false);
531 result = bdsvd::bidiagonalsvddecomposition<Precision>(w, e, n, isupper, false, u, m, a, 0, vt, ncvt);
532 }
533 else
534 {
535
536 //
537 // Large U. Transforming intermediate matrix T2
538 //
539 work.setbounds(1, ap::maxint(m, n));
540 bidiagonal::unpackqfrombidiagonal<Precision>(a, n, n, tauq, n, t2);
541 blas::copymatrix<Precision>(u, 1, m, 1, n, a, 1, m, 1, n);
542 blas::inplacetranspose<Precision>(t2, 1, n, 1, n, work);
543 result = bdsvd::bidiagonalsvddecomposition<Precision>(w, e, n, isupper, false, u, 0, t2, n, vt, ncvt);
544 blas::matrixmatrixmultiply<Precision>(a, 1, m, 1, n, false, t2, 1, n, 1, n, true, amp::ampf<Precision>("1.0"), u, 1, m, 1, n, amp::ampf<Precision>("0.0"), work);
545 }
546 return result;
547 }
548 }
549
550 //
551 // N much larger than M
552 // Use bidiagonal reduction with LQ-decomposition
553 //
554 if( n>amp::ampf<Precision>("1.6")*m )
555 {
556 if( vtneeded==0 )
557 {
558
559 //
560 // No right singular vectors to be computed
561 //
562 lq::lqdecomposition<Precision>(a, m, n, tau);
563 for(i=1; i<=m-1; i++)
564 {
565 for(j=i+1; j<=m; j++)
566 {
567 a(i,j) = 0;
568 }
569 }
570 bidiagonal::tobidiagonal<Precision>(a, m, m, tauq, taup);
571 bidiagonal::unpackqfrombidiagonal<Precision>(a, m, m, tauq, ncu, u);
572 bidiagonal::unpackdiagonalsfrombidiagonal<Precision>(a, m, m, isupper, w, e);
573 work.setbounds(1, m);
574 blas::inplacetranspose<Precision>(u, 1, nru, 1, ncu, work);
575 result = bdsvd::bidiagonalsvddecomposition<Precision>(w, e, m, isupper, false, a, 0, u, nru, vt, 0);
576 blas::inplacetranspose<Precision>(u, 1, nru, 1, ncu, work);
577 return result;
578 }
579 else
580 {
581
582 //
583 // Right singular vectors (may be full matrix VT) to be computed
584 //
585 lq::lqdecomposition<Precision>(a, m, n, tau);
586 lq::unpackqfromlq<Precision>(a, m, n, tau, nrvt, vt);
587 for(i=1; i<=m-1; i++)
588 {
589 for(j=i+1; j<=m; j++)
590 {
591 a(i,j) = 0;
592 }
593 }
594 bidiagonal::tobidiagonal<Precision>(a, m, m, tauq, taup);
595 bidiagonal::unpackqfrombidiagonal<Precision>(a, m, m, tauq, ncu, u);
596 bidiagonal::unpackdiagonalsfrombidiagonal<Precision>(a, m, m, isupper, w, e);
597 work.setbounds(1, ap::maxint(m, n));
598 blas::inplacetranspose<Precision>(u, 1, nru, 1, ncu, work);
599 if( additionalmemory<1 )
600 {
601
602 //
603 // No additional memory available
604 //
605 bidiagonal::multiplybypfrombidiagonal<Precision>(a, m, m, taup, vt, m, n, false, true);
606 result = bdsvd::bidiagonalsvddecomposition<Precision>(w, e, m, isupper, false, a, 0, u, nru, vt, n);
607 }
608 else
609 {
610
611 //
612 // Large VT. Transforming intermediate matrix T2
613 //
614 bidiagonal::unpackptfrombidiagonal<Precision>(a, m, m, taup, m, t2);
615 result = bdsvd::bidiagonalsvddecomposition<Precision>(w, e, m, isupper, false, a, 0, u, nru, t2, m);
616 blas::copymatrix<Precision>(vt, 1, m, 1, n, a, 1, m, 1, n);
617 blas::matrixmatrixmultiply<Precision>(t2, 1, m, 1, m, false, a, 1, m, 1, n, false, amp::ampf<Precision>("1.0"), vt, 1, m, 1, n, amp::ampf<Precision>("0.0"), work);
618 }
619 blas::inplacetranspose<Precision>(u, 1, nru, 1, ncu, work);
620 return result;
621 }
622 }
623
624 //
625 // M<=N
626 // We can use inplace transposition of U to get rid of columnwise operations
627 //
628 if( m<=n )
629 {
630 bidiagonal::tobidiagonal<Precision>(a, m, n, tauq, taup);
631 bidiagonal::unpackqfrombidiagonal<Precision>(a, m, n, tauq, ncu, u);
632 bidiagonal::unpackptfrombidiagonal<Precision>(a, m, n, taup, nrvt, vt);
633 bidiagonal::unpackdiagonalsfrombidiagonal<Precision>(a, m, n, isupper, w, e);
634 work.setbounds(1, m);
635 blas::inplacetranspose<Precision>(u, 1, nru, 1, ncu, work);
636 result = bdsvd::bidiagonalsvddecomposition<Precision>(w, e, minmn, isupper, false, a, 0, u, nru, vt, ncvt);
637 blas::inplacetranspose<Precision>(u, 1, nru, 1, ncu, work);
638 return result;
639 }
640
641 //
642 // Simple bidiagonal reduction
643 //
644 bidiagonal::tobidiagonal<Precision>(a, m, n, tauq, taup);
645 bidiagonal::unpackqfrombidiagonal<Precision>(a, m, n, tauq, ncu, u);
646 bidiagonal::unpackptfrombidiagonal<Precision>(a, m, n, taup, nrvt, vt);
647 bidiagonal::unpackdiagonalsfrombidiagonal<Precision>(a, m, n, isupper, w, e);
648 if( additionalmemory<2 || uneeded==0 )
649 {
650
651 //
652 // We cant use additional memory or there is no need in such operations
653 //
654 result = bdsvd::bidiagonalsvddecomposition<Precision>(w, e, minmn, isupper, false, u, nru, a, 0, vt, ncvt);
655 }
656 else
657 {
658
659 //
660 // We can use additional memory
661 //
662 t2.setbounds(1, minmn, 1, m);
663 blas::copyandtranspose<Precision>(u, 1, m, 1, minmn, t2, 1, minmn, 1, m);
664 result = bdsvd::bidiagonalsvddecomposition<Precision>(w, e, minmn, isupper, false, u, 0, t2, m, vt, ncvt);
665 blas::copyandtranspose<Precision>(t2, 1, minmn, 1, m, u, 1, m, 1, minmn);
666 }
667 return result;
668 }
669} // namespace
670
671#endif
m
int m
Definition: cfEzgcd.cc:128
i
int i
Definition: cfEzgcd.cc:132
tau
void tau(int **points, int sizePoints, int k)
Definition: cfNewtonPolygon.cc:461
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::setbounds
void setbounds(int iLow, int iHigh)
Definition: ap.h:735
ap::template_2d_array
Definition: ap.h:807
ap::template_2d_array::setbounds
void setbounds(int iLow1, int iHigh1, int iLow2, int iHigh2)
Definition: ap.h:890
result
return result
Definition: facAbsBiFact.cc:75
w
const CanonicalForm & w
Definition: facAbsFact.cc:51
j
int j
Definition: facHensel.cc:110
ap::maxint
int maxint(int m1, int m2)
Definition: ap.cpp:162
ap::minint
int minint(int m1, int m2)
Definition: ap.cpp:167
svd
Definition: svd.h:64
svd::rmatrixsvd
bool rmatrixsvd(ap::template_2d_array< amp::ampf< Precision > > a, int m, int n, int uneeded, int vtneeded, int additionalmemory, ap::template_1d_array< amp::ampf< Precision > > &w, ap::template_2d_array< amp::ampf< Precision > > &u, ap::template_2d_array< amp::ampf< Precision > > &vt)
Definition: svd.h:140
svd::svddecomposition
bool svddecomposition(ap::template_2d_array< amp::ampf< Precision > > a, int m, int n, int uneeded, int vtneeded, int additionalmemory, ap::template_1d_array< amp::ampf< Precision > > &w, ap::template_2d_array< amp::ampf< Precision > > &u, ap::template_2d_array< amp::ampf< Precision > > &vt)
Definition: svd.h:408
Generated on Thu Apr 27 2023 08:16:02 for My Project by   doxygen 1.9.4 for Singular