GNU Octave 10.1.0
A high-level interpreted language, primarily intended for numerical computations, mostly compatible with Matlab
 
Loading...
Searching...
No Matches
gsvd.cc
Go to the documentation of this file.
1////////////////////////////////////////////////////////////////////////
2//
3// Copyright (C) 1997-2025 The Octave Project Developers
4//
5// See the file COPYRIGHT.md in the top-level directory of this
6// distribution or <https://octave.org/copyright/>.
7//
8// This file is part of Octave.
9//
10// Octave is free software: you can redistribute it and/or modify it
11// under the terms of the GNU General Public License as published by
12// the Free Software Foundation, either version 3 of the License, or
13// (at your option) any later version.
14//
15// Octave is distributed in the hope that it will be useful, but
16// WITHOUT ANY WARRANTY; without even the implied warranty of
17// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
18// GNU General Public License for more details.
19//
20// You should have received a copy of the GNU General Public License
21// along with Octave; see the file COPYING. If not, see
22// <https://www.gnu.org/licenses/>.
23//
24////////////////////////////////////////////////////////////////////////
25
26#if defined (HAVE_CONFIG_H)
27# include "config.h"
28#endif
29
30#include <algorithm>
31#include <unordered_map>
32
33#include "CMatrix.h"
34#include "dDiagMatrix.h"
35#include "dMatrix.h"
36#include "fCMatrix.h"
37#include "fDiagMatrix.h"
38#include "fMatrix.h"
39#include "gsvd.h"
40#include "lo-error.h"
41#include "lo-lapack-proto.h"
42#include "oct-locbuf.h"
43#include "oct-shlib.h"
44
45OCTAVE_BEGIN_NAMESPACE(octave)
46
47static std::unordered_map<std::string, void *> gsvd_fcn;
48
49static bool have_DGGSVD3 = false;
50static bool gsvd_initialized = false;
51
52/* Hack to stringize results of F77_FUNC macro. */
53#define xSTRINGIZE(x) #x
54#define STRINGIZE(x) xSTRINGIZE(x)
55
56static void
57initialize_gsvd ()
58{
59 if (gsvd_initialized)
60 return;
61
62 dynamic_library libs ("");
63 if (! libs)
64 (*current_liboctave_error_handler)
65 ("gsvd: unable to query LAPACK library");
66
67 have_DGGSVD3 = (libs.search (STRINGIZE (F77_FUNC (dggsvd3, DGGSVD3)))
68 != nullptr);
69
70 if (have_DGGSVD3)
71 {
72 gsvd_fcn["dg"] = libs.search (STRINGIZE (F77_FUNC (dggsvd3, DGGSVD3)));
73 gsvd_fcn["sg"] = libs.search (STRINGIZE (F77_FUNC (sggsvd3, SGGSVD3)));
74 gsvd_fcn["zg"] = libs.search (STRINGIZE (F77_FUNC (zggsvd3, ZGGSVD3)));
75 gsvd_fcn["cg"] = libs.search (STRINGIZE (F77_FUNC (cggsvd3, CGGSVD3)));
76 }
77 else
78 {
79 gsvd_fcn["dg"] = libs.search (STRINGIZE (F77_FUNC (dggsvd, DGGSVD)));
80 gsvd_fcn["sg"] = libs.search (STRINGIZE (F77_FUNC (sggsvd, SGGSVD)));
81 gsvd_fcn["zg"] = libs.search (STRINGIZE (F77_FUNC (zggsvd, ZGGSVD)));
82 gsvd_fcn["cg"] = libs.search (STRINGIZE (F77_FUNC (cggsvd, CGGSVD)));
83 }
84
85 gsvd_initialized = true;
86}
87
88/* Clean up macro namespace as soon as we are done using them */
89#undef xSTRINGIZE
90#undef STRINGIZE
91
92template<class T1>
93struct real_ggsvd_ptr
94{
95 typedef F77_RET_T (*type)
96 (F77_CONST_CHAR_ARG_DECL, // JOBU
97 F77_CONST_CHAR_ARG_DECL, // JOBV
98 F77_CONST_CHAR_ARG_DECL, // JOBQ
99 const F77_INT&, // M
100 const F77_INT&, // N
101 const F77_INT&, // P
102 F77_INT&, // K
103 F77_INT&, // L
104 T1 *, // A(LDA,N)
105 const F77_INT&, // LDA
106 T1 *, // B(LDB,N)
107 const F77_INT&, // LDB
108 T1 *, // ALPHA(N)
109 T1 *, // BETA(N)
110 T1 *, // U(LDU,M)
111 const F77_INT&, // LDU
112 T1 *, // V(LDV,P)
113 const F77_INT&, // LDV
114 T1 *, // Q(LDQ,N)
115 const F77_INT&, // LDQ
116 T1 *, // WORK
117 F77_INT *, // IWORK(N)
118 F77_INT& // INFO
119 F77_CHAR_ARG_LEN_DECL
120 F77_CHAR_ARG_LEN_DECL
121 F77_CHAR_ARG_LEN_DECL);
122};
123
124template<class T1>
125struct real_ggsvd3_ptr
126{
127 typedef F77_RET_T (*type)
128 (F77_CONST_CHAR_ARG_DECL, // JOBU
129 F77_CONST_CHAR_ARG_DECL, // JOBV
130 F77_CONST_CHAR_ARG_DECL, // JOBQ
131 const F77_INT&, // M
132 const F77_INT&, // N
133 const F77_INT&, // P
134 F77_INT&, // K
135 F77_INT&, // L
136 T1 *, // A(LDA,N)
137 const F77_INT&, // LDA
138 T1 *, // B(LDB,N)
139 const F77_INT&, // LDB
140 T1 *, // ALPHA(N)
141 T1 *, // BETA(N)
142 T1 *, // U(LDU,M)
143 const F77_INT&, // LDU
144 T1 *, // V(LDV,P)
145 const F77_INT&, // LDV
146 T1 *, // Q(LDQ,N)
147 const F77_INT&, // LDQ
148 T1 *, // WORK
149 const F77_INT&, // LWORK
150 F77_INT *, // IWORK(N)
151 F77_INT& // INFO
152 F77_CHAR_ARG_LEN_DECL
153 F77_CHAR_ARG_LEN_DECL
154 F77_CHAR_ARG_LEN_DECL);
155};
156
157template<class T1, class T2>
158struct comp_ggsvd_ptr
159{
160 typedef F77_RET_T (*type)
161 (F77_CONST_CHAR_ARG_DECL, // JOBU
162 F77_CONST_CHAR_ARG_DECL, // JOBV
163 F77_CONST_CHAR_ARG_DECL, // JOBQ
164 const F77_INT&, // M
165 const F77_INT&, // N
166 const F77_INT&, // P
167 F77_INT&, // K
168 F77_INT&, // L
169 T1 *, // A(LDA,N)
170 const F77_INT&, // LDA
171 T1 *, // B(LDB,N)
172 const F77_INT&, // LDB
173 T2 *, // ALPHA(N)
174 T2 *, // BETA(N)
175 T1 *, // U(LDU,M)
176 const F77_INT&, // LDU
177 T1 *, // V(LDV,P)
178 const F77_INT&, // LDV
179 T1 *, // Q(LDQ,N)
180 const F77_INT&, // LDQ
181 T1 *, // WORK
182 T2 *, // RWORK
183 F77_INT *, // IWORK(N)
184 F77_INT& // INFO
185 F77_CHAR_ARG_LEN_DECL
186 F77_CHAR_ARG_LEN_DECL
187 F77_CHAR_ARG_LEN_DECL);
188};
189
190template<class T1, class T2>
191struct comp_ggsvd3_ptr
192{
193 typedef F77_RET_T (*type)
194 (F77_CONST_CHAR_ARG_DECL, // JOBU
195 F77_CONST_CHAR_ARG_DECL, // JOBV
196 F77_CONST_CHAR_ARG_DECL, // JOBQ
197 const F77_INT&, // M
198 const F77_INT&, // N
199 const F77_INT&, // P
200 F77_INT&, // K
201 F77_INT&, // L
202 T1 *, // A(LDA,N)
203 const F77_INT&, // LDA
204 T1 *, // B(LDB,N)
205 const F77_INT&, // LDB
206 T2 *, // ALPHA(N)
207 T2 *, // BETA(N)
208 T1 *, // U(LDU,M)
209 const F77_INT&, // LDU
210 T1 *, // V(LDV,P)
211 const F77_INT&, // LDV
212 T1 *, // Q(LDQ,N)
213 const F77_INT&, // LDQ
214 T1 *, // WORK
215 const F77_INT&, // LWORK
216 T2 *, // RWORK
217 F77_INT *, // IWORK(N)
218 F77_INT& // INFO
219 F77_CHAR_ARG_LEN_DECL
220 F77_CHAR_ARG_LEN_DECL
221 F77_CHAR_ARG_LEN_DECL);
222};
223
224// template specializations
225typedef real_ggsvd_ptr<F77_DBLE>::type dggsvd_type;
226typedef real_ggsvd3_ptr<F77_DBLE>::type dggsvd3_type;
227typedef real_ggsvd_ptr<F77_REAL>::type sggsvd_type;
228typedef real_ggsvd3_ptr<F77_REAL>::type sggsvd3_type;
229typedef comp_ggsvd_ptr<F77_DBLE_CMPLX, F77_DBLE>::type zggsvd_type;
230typedef comp_ggsvd3_ptr<F77_DBLE_CMPLX, F77_DBLE>::type zggsvd3_type;
231typedef comp_ggsvd_ptr<F77_CMPLX, F77_REAL>::type cggsvd_type;
232typedef comp_ggsvd3_ptr<F77_CMPLX, F77_REAL>::type cggsvd3_type;
233
234OCTAVE_BEGIN_NAMESPACE(math)
235
236template <>
237void
238gsvd<Matrix>::ggsvd (char& jobu, char& jobv, char& jobq, F77_INT m,
239 F77_INT n, F77_INT p, F77_INT& k, F77_INT& l,
240 double *tmp_dataA, F77_INT m1, double *tmp_dataB,
241 F77_INT p1, Matrix& alpha, Matrix& beta, double *u,
242 F77_INT nrow_u, double *v, F77_INT nrow_v, double *q,
243 F77_INT nrow_q, double *work, F77_INT lwork,
244 F77_INT *iwork, F77_INT& info)
245{
246 if (! gsvd_initialized)
247 initialize_gsvd ();
248
249 if (have_DGGSVD3)
250 {
251 dggsvd3_type f_ptr = reinterpret_cast<dggsvd3_type> (gsvd_fcn["dg"]);
252 f_ptr (F77_CONST_CHAR_ARG2 (&jobu, 1),
253 F77_CONST_CHAR_ARG2 (&jobv, 1),
254 F77_CONST_CHAR_ARG2 (&jobq, 1),
255 m, n, p, k, l, tmp_dataA, m1, tmp_dataB, p1,
256 alpha.rwdata (), beta.rwdata (),
257 u, nrow_u, v, nrow_v, q, nrow_q,
258 work, lwork, iwork, info
259 F77_CHAR_ARG_LEN (1)
260 F77_CHAR_ARG_LEN (1)
261 F77_CHAR_ARG_LEN (1));
262 }
263 else
264 {
265 dggsvd_type f_ptr = reinterpret_cast<dggsvd_type> (gsvd_fcn["dg"]);
266 f_ptr (F77_CONST_CHAR_ARG2 (&jobu, 1),
267 F77_CONST_CHAR_ARG2 (&jobv, 1),
268 F77_CONST_CHAR_ARG2 (&jobq, 1),
269 m, n, p, k, l, tmp_dataA, m1, tmp_dataB, p1,
270 alpha.rwdata (), beta.rwdata (),
271 u, nrow_u, v, nrow_v, q, nrow_q,
272 work, iwork, info
273 F77_CHAR_ARG_LEN (1)
274 F77_CHAR_ARG_LEN (1)
275 F77_CHAR_ARG_LEN (1));
276 }
277}
278
279template <>
280void
281gsvd<FloatMatrix>::ggsvd (char& jobu, char& jobv, char& jobq, F77_INT m,
282 F77_INT n, F77_INT p, F77_INT& k, F77_INT& l,
283 float *tmp_dataA, F77_INT m1, float *tmp_dataB,
284 F77_INT p1, FloatMatrix& alpha,
285 FloatMatrix& beta, float *u, F77_INT nrow_u,
286 float *v, F77_INT nrow_v, float *q,
287 F77_INT nrow_q, float *work,
288 F77_INT lwork, F77_INT *iwork, F77_INT& info)
289{
290 if (! gsvd_initialized)
291 initialize_gsvd ();
292
293 if (have_DGGSVD3)
294 {
295 sggsvd3_type f_ptr = reinterpret_cast<sggsvd3_type> (gsvd_fcn["sg"]);
296 f_ptr (F77_CONST_CHAR_ARG2 (&jobu, 1),
297 F77_CONST_CHAR_ARG2 (&jobv, 1),
298 F77_CONST_CHAR_ARG2 (&jobq, 1),
299 m, n, p, k, l, tmp_dataA, m1, tmp_dataB, p1,
300 alpha.rwdata (), beta.rwdata (),
301 u, nrow_u, v, nrow_v, q, nrow_q,
302 work, lwork, iwork, info
303 F77_CHAR_ARG_LEN (1)
304 F77_CHAR_ARG_LEN (1)
305 F77_CHAR_ARG_LEN (1));
306 }
307 else
308 {
309 sggsvd_type f_ptr = reinterpret_cast<sggsvd_type> (gsvd_fcn["sg"]);
310 f_ptr (F77_CONST_CHAR_ARG2 (&jobu, 1),
311 F77_CONST_CHAR_ARG2 (&jobv, 1),
312 F77_CONST_CHAR_ARG2 (&jobq, 1),
313 m, n, p, k, l, tmp_dataA, m1, tmp_dataB, p1,
314 alpha.rwdata (), beta.rwdata (),
315 u, nrow_u, v, nrow_v, q, nrow_q,
316 work, iwork, info
317 F77_CHAR_ARG_LEN (1)
318 F77_CHAR_ARG_LEN (1)
319 F77_CHAR_ARG_LEN (1));
320 }
321}
322
323template <>
324void
325gsvd<ComplexMatrix>::ggsvd (char& jobu, char& jobv, char& jobq,
326 F77_INT m, F77_INT n, F77_INT p, F77_INT& k,
327 F77_INT& l, Complex *tmp_dataA, F77_INT m1,
328 Complex *tmp_dataB, F77_INT p1, Matrix& alpha,
329 Matrix& beta, Complex *u, F77_INT nrow_u,
330 Complex *v, F77_INT nrow_v, Complex *q,
331 F77_INT nrow_q, Complex *work,
332 F77_INT lwork, F77_INT *iwork, F77_INT& info)
333{
334 if (! gsvd_initialized)
335 initialize_gsvd ();
336
337 OCTAVE_LOCAL_BUFFER(double, rwork, 2*n);
338
339 if (have_DGGSVD3)
340 {
341 zggsvd3_type f_ptr = reinterpret_cast<zggsvd3_type> (gsvd_fcn["zg"]);
342 f_ptr (F77_CONST_CHAR_ARG2 (&jobu, 1),
343 F77_CONST_CHAR_ARG2 (&jobv, 1),
344 F77_CONST_CHAR_ARG2 (&jobq, 1),
345 m, n, p, k, l,
346 F77_DBLE_CMPLX_ARG (tmp_dataA), m1,
347 F77_DBLE_CMPLX_ARG (tmp_dataB), p1,
348 alpha.rwdata (), beta.rwdata (),
349 F77_DBLE_CMPLX_ARG (u), nrow_u,
350 F77_DBLE_CMPLX_ARG (v), nrow_v,
351 F77_DBLE_CMPLX_ARG (q), nrow_q,
352 F77_DBLE_CMPLX_ARG (work),
353 lwork, rwork, iwork, info
354 F77_CHAR_ARG_LEN (1)
355 F77_CHAR_ARG_LEN (1)
356 F77_CHAR_ARG_LEN (1));
357 }
358 else
359 {
360 zggsvd_type f_ptr = reinterpret_cast<zggsvd_type> (gsvd_fcn["zg"]);
361 f_ptr (F77_CONST_CHAR_ARG2 (&jobu, 1),
362 F77_CONST_CHAR_ARG2 (&jobv, 1),
363 F77_CONST_CHAR_ARG2 (&jobq, 1),
364 m, n, p, k, l,
365 F77_DBLE_CMPLX_ARG (tmp_dataA), m1,
366 F77_DBLE_CMPLX_ARG (tmp_dataB), p1,
367 alpha.rwdata (), beta.rwdata (),
368 F77_DBLE_CMPLX_ARG (u), nrow_u,
369 F77_DBLE_CMPLX_ARG (v), nrow_v,
370 F77_DBLE_CMPLX_ARG (q), nrow_q,
371 F77_DBLE_CMPLX_ARG (work),
372 rwork, iwork, info
373 F77_CHAR_ARG_LEN (1)
374 F77_CHAR_ARG_LEN (1)
375 F77_CHAR_ARG_LEN (1));
376 }
377}
378
379template <>
380void
381gsvd<FloatComplexMatrix>::ggsvd (char& jobu, char& jobv, char& jobq,
382 F77_INT m, F77_INT n, F77_INT p,
383 F77_INT& k, F77_INT& l,
384 FloatComplex *tmp_dataA, F77_INT m1,
385 FloatComplex *tmp_dataB, F77_INT p1,
386 FloatMatrix& alpha, FloatMatrix& beta,
387 FloatComplex *u, F77_INT nrow_u,
388 FloatComplex *v, F77_INT nrow_v,
389 FloatComplex *q, F77_INT nrow_q,
390 FloatComplex *work, F77_INT lwork,
391 F77_INT *iwork, F77_INT& info)
392{
393 if (! gsvd_initialized)
394 initialize_gsvd ();
395
396 OCTAVE_LOCAL_BUFFER(float, rwork, 2*n);
397
398 if (have_DGGSVD3)
399 {
400 cggsvd3_type f_ptr = reinterpret_cast<cggsvd3_type> (gsvd_fcn["cg"]);
401 f_ptr (F77_CONST_CHAR_ARG2 (&jobu, 1),
402 F77_CONST_CHAR_ARG2 (&jobv, 1),
403 F77_CONST_CHAR_ARG2 (&jobq, 1),
404 m, n, p, k, l,
405 F77_CMPLX_ARG (tmp_dataA), m1,
406 F77_CMPLX_ARG (tmp_dataB), p1,
407 alpha.rwdata (), beta.rwdata (),
408 F77_CMPLX_ARG (u), nrow_u,
409 F77_CMPLX_ARG (v), nrow_v,
410 F77_CMPLX_ARG (q), nrow_q,
411 F77_CMPLX_ARG (work), lwork,
412 rwork, iwork, info
413 F77_CHAR_ARG_LEN (1)
414 F77_CHAR_ARG_LEN (1)
415 F77_CHAR_ARG_LEN (1));
416 }
417 else
418 {
419 cggsvd_type f_ptr = reinterpret_cast<cggsvd_type> (gsvd_fcn["cg"]);
420 f_ptr (F77_CONST_CHAR_ARG2 (&jobu, 1),
421 F77_CONST_CHAR_ARG2 (&jobv, 1),
422 F77_CONST_CHAR_ARG2 (&jobq, 1),
423 m, n, p, k, l,
424 F77_CMPLX_ARG (tmp_dataA), m1,
425 F77_CMPLX_ARG (tmp_dataB), p1,
426 alpha.rwdata (), beta.rwdata (),
427 F77_CMPLX_ARG (u), nrow_u,
428 F77_CMPLX_ARG (v), nrow_v,
429 F77_CMPLX_ARG (q), nrow_q,
430 F77_CMPLX_ARG (work),
431 rwork, iwork, info
432 F77_CHAR_ARG_LEN (1)
433 F77_CHAR_ARG_LEN (1)
434 F77_CHAR_ARG_LEN (1));
435 }
436}
437
438template <typename T>
439T
440gsvd<T>::left_singular_matrix_A () const
441{
442 if (m_type == gsvd::Type::sigma_only)
443 (*current_liboctave_error_handler)
444 ("gsvd: U not computed because type == gsvd::sigma_only");
445
446 return m_left_smA;
447}
448
449template <typename T>
450T
451gsvd<T>::left_singular_matrix_B () const
452{
453 if (m_type == gsvd::Type::sigma_only)
454 (*current_liboctave_error_handler)
455 ("gsvd: V not computed because type == gsvd::sigma_only");
456
457 return m_left_smB;
458}
459
460template <typename T>
461T
462gsvd<T>::right_singular_matrix () const
463{
464 if (m_type == gsvd::Type::sigma_only)
465 (*current_liboctave_error_handler)
466 ("gsvd: X not computed because type == gsvd::sigma_only");
467
468 return m_right_sm;
469}
470
471template <typename T>
472gsvd<T>::gsvd (const T& a, const T& b, gsvd::Type gsvd_type)
473{
474 if (a.isempty () || b.isempty ())
475 (*current_liboctave_error_handler)
476 ("gsvd: A and B cannot be empty matrices");
477
478 F77_INT info;
479
480 F77_INT m = to_f77_int (a.rows ());
481 F77_INT n = to_f77_int (a.cols ());
482 F77_INT p = to_f77_int (b.rows ());
483
484 T atmp = a;
485 P *tmp_dataA = atmp.rwdata ();
486
487 T btmp = b;
488 P *tmp_dataB = btmp.rwdata ();
489
490 char jobu = 'U';
491 char jobv = 'V';
492 char jobq = 'Q';
493
494 F77_INT nrow_u = m;
495 F77_INT nrow_v = p;
496 F77_INT nrow_q = n;
497
498 F77_INT k, l;
499
500 switch (gsvd_type)
501 {
502 case gsvd<T>::Type::sigma_only:
503
504 // FIXME: In LAPACK 3.0, problem below seems to be fixed,
505 // so we now set jobX = 'N'.
506 //
507 // For calculating sigma_only, both jobu and jobv should be 'N', but
508 // there seems to be a bug in dgesvd from LAPACK V2.0. To
509 // demonstrate the bug, set both jobu and jobv to 'N' and find the
510 // singular values of [eye(3), eye(3)]. The result is
511 // [-sqrt(2), -sqrt(2), -sqrt(2)].
512
513 jobu = jobv = jobq = 'N';
514 nrow_u = nrow_v = nrow_q = 1;
515 break;
516
517 default:
518 break;
519 }
520
521 m_type = gsvd_type;
522
523 if (jobu != 'N')
524 m_left_smA.resize (nrow_u, m);
525
526 P *u = m_left_smA.rwdata ();
527
528 if (jobv != 'N')
529 m_left_smB.resize (nrow_v, p);
530
531 P *v = m_left_smB.rwdata ();
532
533 if (jobq != 'N')
534 m_right_sm.resize (nrow_q, n);
535
536 P *q = m_right_sm.rwdata ();
537
538 real_matrix alpha (n, 1);
539 real_matrix beta (n, 1);
540
541 OCTAVE_LOCAL_BUFFER(F77_INT, iwork, n);
542
543 if (! gsvd_initialized)
544 initialize_gsvd ();
545
546 F77_INT lwork;
547 if (have_DGGSVD3)
548 {
549 lwork = -1;
550 P work_tmp;
551
552 gsvd<T>::ggsvd (jobu, jobv, jobq, m, n, p, k, l,
553 tmp_dataA, m, tmp_dataB, p,
554 alpha, beta, u, nrow_u, v, nrow_v, q, nrow_q,
555 &work_tmp, lwork, iwork, info);
556
557 lwork = static_cast<F77_INT> (std::abs (work_tmp));
558 }
559 else
560 {
561 lwork = std::max ({3*n, m, p}) + n;
562 }
563 info = 0;
564
565 OCTAVE_LOCAL_BUFFER(P, work, lwork);
566
567 gsvd<T>::ggsvd (jobu, jobv, jobq, m, n, p, k, l,
568 tmp_dataA, m, tmp_dataB, p,
569 alpha, beta, u, nrow_u, v, nrow_v, q, nrow_q,
570 work, lwork, iwork, info);
571
572 if (info < 0)
573 (*current_liboctave_error_handler)
574 ("*ggsvd.f: argument %" OCTAVE_F77_INT_TYPE_FORMAT " illegal",
575 -info);
576
577 if (info > 0)
578 (*current_liboctave_error_handler)
579 ("*ggsvd.f: Jacobi-type procedure failed to converge");
580
581 F77_INT i, j;
582
583 if (gsvd_type != gsvd<T>::Type::sigma_only)
584 {
585 // Size according to LAPACK is k+l,k+l, but this needs
586 // to be nxn for Matlab compatibility.
587 T R (n, n, 0.0);
588 int astart = n-k-l;
589 if (m - k - l >= 0)
590 {
591 // R is stored in A(1:K+L,N-K-L+1:N)
592 for (i = 0; i < k+l; i++)
593 for (j = 0; j < k+l; j++)
594 R.xelem (i, j) = atmp.xelem (i, astart + j);
595 }
596 else
597 {
598 // (R11 R12 R13 ) is stored in A(1:M, N-K-L+1:N)
599 // ( 0 R22 R23 )
600 for (i = 0; i < m; i++)
601 for (j = 0; j < k+l; j++)
602 R.xelem (i, j) = atmp.xelem (i, astart + j);
603 // and R33 is stored in B(M-K+1:L,N+M-K-L+1:N)
604 for (i = m; i < k + l; i++)
605 for (j = n - l - k + m; j < n; j++)
606 R.xelem (i, j) = btmp.xelem (i - k, j);
607 }
608
609 // Output X = Q*R'
610 // FIXME: Is there a way to call BLAS multiply directly
611 // with flags so that R is transposed?
612 m_right_sm = m_right_sm * R.hermitian ();
613 }
614
615 // Fill in C and S
616 F77_INT rank;
617 bool fill_ptn;
618 if (m-k-l >= 0)
619 {
620 rank = l;
621 fill_ptn = true;
622 }
623 else
624 {
625 rank = m-k;
626 fill_ptn = false;
627 }
628
629 if (gsvd_type == gsvd<T>::Type::sigma_only)
630 {
631 // Return column vector with results
632 m_sigmaA.resize (k+l, 1);
633 m_sigmaB.resize (k+l, 1);
634
635 if (fill_ptn)
636 {
637 for (i = 0; i < k; i++)
638 {
639 m_sigmaA.xelem (i) = 1.0;
640 m_sigmaB.xelem (i) = 0.0;
641 }
642 for (i = k, j = k+l-1; i < k+l; i++, j--)
643 {
644 m_sigmaA.xelem (i) = alpha.xelem (i);
645 m_sigmaB.xelem (i) = beta.xelem (i);
646 }
647 }
648 else
649 {
650 for (i = 0; i < k; i++)
651 {
652 m_sigmaA.xelem (i) = 1.0;
653 m_sigmaB.xelem (i) = 0.0;
654 }
655 for (i = k; i < m; i++)
656 {
657 m_sigmaA.xelem (i) = alpha.xelem (i);
658 m_sigmaB.xelem (i) = beta.xelem (i);
659 }
660 for (i = m; i < k+l; i++)
661 {
662 m_sigmaA.xelem (i) = 0.0;
663 m_sigmaB.xelem (i) = 1.0;
664 }
665 }
666 }
667 else // returning all matrices
668 {
669 // Number of columns according to LAPACK is k+l, but this needs
670 // to be n for Matlab compatibility.
671 m_sigmaA.resize (m, n);
672 m_sigmaB.resize (p, n);
673
674 for (i = 0; i < k; i++)
675 m_sigmaA.xelem (i, i) = 1.0;
676
677 for (i = 0; i < rank; i++)
678 {
679 m_sigmaA.xelem (k+i, k+i) = alpha.xelem (k+i);
680 m_sigmaB.xelem (i, k+i) = beta.xelem (k+i);
681 }
682
683 if (! fill_ptn)
684 {
685 for (i = m; i < n; i++)
686 m_sigmaB.xelem (i-k, i) = 1.0;
687 }
688
689 }
690}
691
692// Instantiations needed in octave::math namespace.
693template class gsvd<Matrix>;
694template class gsvd<FloatMatrix>;
695template class gsvd<ComplexMatrix>;
696template class gsvd<FloatComplexMatrix>;
697
698OCTAVE_END_NAMESPACE(math)
699OCTAVE_END_NAMESPACE(octave)
Array::rwdata
T * rwdata()
Size of the specified dimension.
Definition Array-base.cc:1787
dynamic_library
Definition oct-shlib.h:42
dynamic_library::search
void * search(const std::string &nm, const name_mangler &mangler=name_mangler()) const
Definition oct-shlib.h:178
gsvd
Definition gsvd.h:37
gsvd::left_singular_matrix_B
T left_singular_matrix_B() const
Definition gsvd.cc:451
gsvd::left_singular_matrix_A
T left_singular_matrix_A() const
Definition gsvd.cc:440
gsvd::right_singular_matrix
T right_singular_matrix() const
Definition gsvd.cc:462
gsvd::gsvd
gsvd()
Definition gsvd.h:47
gsvd::Type
Type
Definition gsvd.h:41
OCTAVE_BEGIN_NAMESPACE
OCTAVE_BEGIN_NAMESPACE(octave) static octave_value daspk_fcn
F77_DBLE_CMPLX_ARG
#define F77_DBLE_CMPLX_ARG(x)
Definition f77-fcn.h:316
F77_CMPLX_ARG
#define F77_CMPLX_ARG(x)
Definition f77-fcn.h:310
F77_INT
octave_f77_int_type F77_INT
Definition f77-fcn.h:306
dggsvd_type
real_ggsvd_ptr< F77_DBLE >::type dggsvd_type
Definition gsvd.cc:225
STRINGIZE
#define STRINGIZE(x)
Definition gsvd.cc:54
sggsvd3_type
real_ggsvd3_ptr< F77_REAL >::type sggsvd3_type
Definition gsvd.cc:228
cggsvd3_type
comp_ggsvd3_ptr< F77_CMPLX, F77_REAL >::type cggsvd3_type
Definition gsvd.cc:232
zggsvd_type
comp_ggsvd_ptr< F77_DBLE_CMPLX, F77_DBLE >::type zggsvd_type
Definition gsvd.cc:229
cggsvd_type
comp_ggsvd_ptr< F77_CMPLX, F77_REAL >::type cggsvd_type
Definition gsvd.cc:231
zggsvd3_type
comp_ggsvd3_ptr< F77_DBLE_CMPLX, F77_DBLE >::type zggsvd3_type
Definition gsvd.cc:230
dggsvd3_type
real_ggsvd3_ptr< F77_DBLE >::type dggsvd3_type
Definition gsvd.cc:226
sggsvd_type
real_ggsvd_ptr< F77_REAL >::type sggsvd_type
Definition gsvd.cc:227
F77_CONST_CHAR_ARG_DECL
F77_RET_T F77_CONST_CHAR_ARG_DECL
Definition lo-arpack-proto.h:42
F77_CHAR_ARG_LEN_DECL
F77_RET_T const F77_INT F77_INT const F77_DBLE F77_DBLE const F77_INT F77_DBLE const F77_INT F77_INT F77_INT F77_DBLE F77_DBLE const F77_INT F77_INT &F77_CHAR_ARG_LEN_DECL F77_CHAR_ARG_LEN_DECL
Definition lo-arpack-proto.h:45
F77_RET_T
F77_RET_T(F77_CONST_CHAR_ARG_DECL, F77_CONST_CHAR_ARG_DECL, F77_CONST_CHAR_ARG_DECL, const F77_INT &, const F77_INT &, const F77_INT &, F77_INT &, F77_INT &, F77_DBLE *, const F77_INT &, F77_DBLE *, const F77_INT &, F77_DBLE *, F77_DBLE *, F77_DBLE *, const F77_INT &, F77_DBLE *, const F77_INT &, F77_DBLE *, const F77_INT &, F77_DBLE *, F77_INT *, F77_INT &F77_CHAR_ARG_LEN_DECL F77_CHAR_ARG_LEN_DECL F77_CHAR_ARG_LEN_DECL)
Complex
std::complex< double > Complex
Definition oct-cmplx.h:33
FloatComplex
std::complex< float > FloatComplex
Definition oct-cmplx.h:34
OCTAVE_LOCAL_BUFFER
#define OCTAVE_LOCAL_BUFFER(T, buf, size)
Definition oct-locbuf.h:44
F77_FUNC
F77_RET_T F77_FUNC(xerbla, XERBLA)(F77_CONST_CHAR_ARG_DEF(s_arg