Mercurial > octave

comparison liboctave/numeric/gsvd.cc @ 31607:aac27ad79be6 stable

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
maint: Re-indent code after switch to using namespace macros. * build-env.h, build-env.in.cc, Cell.h, __betainc__.cc, __eigs__.cc, __ftp__.cc, __ichol__.cc, __ilu__.cc, __isprimelarge__.cc, __magick_read__.cc, __pchip_deriv__.cc, amd.cc, base-text-renderer.cc, base-text-renderer.h, besselj.cc, bitfcns.cc, bsxfun.cc, c-file-ptr-stream.h, call-stack.cc, call-stack.h, ccolamd.cc, cellfun.cc, chol.cc, colamd.cc, dasrt.cc, data.cc, debug.cc, defaults.cc, defaults.h, det.cc, display.cc, display.h, dlmread.cc, dynamic-ld.cc, dynamic-ld.h, ellipj.cc, environment.cc, environment.h, error.cc, error.h, errwarn.h, event-manager.cc, event-manager.h, event-queue.cc, event-queue.h, fcn-info.cc, fcn-info.h, fft.cc, fft2.cc, file-io.cc, filter.cc, find.cc, ft-text-renderer.cc, ft-text-renderer.h, gcd.cc, gl-render.cc, gl-render.h, gl2ps-print.cc, gl2ps-print.h, graphics-toolkit.cc, graphics-toolkit.h, graphics.cc, gsvd.cc, gtk-manager.cc, gtk-manager.h, help.cc, help.h, hook-fcn.cc, hook-fcn.h, input.cc, input.h, interpreter-private.cc, interpreter-private.h, interpreter.cc, interpreter.h, inv.cc, jsondecode.cc, jsonencode.cc, latex-text-renderer.cc, latex-text-renderer.h, load-path.cc, load-path.h, load-save.cc, load-save.h, lookup.cc, ls-hdf5.cc, ls-mat4.cc, ls-mat5.cc, lsode.cc, lu.cc, mappers.cc, matrix_type.cc, max.cc, mex.cc, mexproto.h, mxarray.h, mxtypes.in.h, oct-errno.in.cc, oct-hdf5-types.cc, oct-hist.cc, oct-hist.h, oct-map.cc, oct-map.h, oct-opengl.h, oct-prcstrm.h, oct-process.cc, oct-process.h, oct-stdstrm.h, oct-stream.cc, oct-stream.h, oct-strstrm.h, octave-default-image.h, ordqz.cc, ordschur.cc, pager.cc, pager.h, pinv.cc, pow2.cc, pr-output.cc, psi.cc, qr.cc, quadcc.cc, rand.cc, regexp.cc, settings.cc, settings.h, sighandlers.cc, sighandlers.h, sparse-xpow.cc, sqrtm.cc, stack-frame.cc, stack-frame.h, stream-euler.cc, strfns.cc, svd.cc, syminfo.cc, syminfo.h, symrcm.cc, symrec.cc, symrec.h, symscope.cc, symscope.h, symtab.cc, symtab.h, sysdep.cc, sysdep.h, text-engine.cc, text-engine.h, text-renderer.cc, text-renderer.h, time.cc, toplev.cc, typecast.cc, url-handle-manager.cc, url-handle-manager.h, urlwrite.cc, utils.cc, utils.h, variables.cc, variables.h, xdiv.cc, __delaunayn__.cc, __init_fltk__.cc, __init_gnuplot__.cc, __ode15__.cc, __voronoi__.cc, audioread.cc, convhulln.cc, gzip.cc, cdef-class.cc, cdef-class.h, cdef-fwd.h, cdef-manager.cc, cdef-manager.h, cdef-method.cc, cdef-method.h, cdef-object.cc, cdef-object.h, cdef-package.cc, cdef-package.h, cdef-property.cc, cdef-property.h, cdef-utils.cc, cdef-utils.h, ov-base-diag.cc, ov-base-int.cc, ov-base-mat.cc, ov-base-mat.h, ov-base-scalar.cc, ov-base.cc, ov-base.h, ov-bool-mat.cc, ov-bool-mat.h, ov-bool-sparse.cc, ov-bool.cc, ov-builtin.h, ov-cell.cc, ov-ch-mat.cc, ov-class.cc, ov-class.h, ov-classdef.cc, ov-classdef.h, ov-complex.cc, ov-cx-diag.cc, ov-cx-mat.cc, ov-cx-sparse.cc, ov-dld-fcn.cc, ov-dld-fcn.h, ov-fcn-handle.cc, ov-fcn-handle.h, ov-fcn.h, ov-float.cc, ov-flt-complex.cc, ov-flt-cx-diag.cc, ov-flt-cx-mat.cc, ov-flt-re-diag.cc, ov-flt-re-mat.cc, ov-flt-re-mat.h, ov-intx.h, ov-java.cc, ov-lazy-idx.cc, ov-legacy-range.cc, ov-magic-int.cc, ov-mex-fcn.cc, ov-mex-fcn.h, ov-null-mat.cc, ov-perm.cc, ov-range.cc, ov-re-diag.cc, ov-re-mat.cc, ov-re-mat.h, ov-re-sparse.cc, ov-scalar.cc, ov-str-mat.cc, ov-struct.cc, ov-typeinfo.cc, ov-typeinfo.h, ov-usr-fcn.cc, ov-usr-fcn.h, ov.cc, ov.h, ovl.h, octave.cc, octave.h, op-b-sbm.cc, op-bm-sbm.cc, op-cs-scm.cc, op-fm-fcm.cc, op-fs-fcm.cc, op-s-scm.cc, op-scm-cs.cc, op-scm-s.cc, op-sm-cs.cc, ops.h, anon-fcn-validator.cc, anon-fcn-validator.h, bp-table.cc, bp-table.h, comment-list.cc, comment-list.h, filepos.h, lex.h, oct-lvalue.cc, oct-lvalue.h, parse.h, profiler.cc, profiler.h, pt-anon-scopes.cc, pt-anon-scopes.h, pt-arg-list.cc, pt-arg-list.h, pt-args-block.cc, pt-args-block.h, pt-array-list.cc, pt-array-list.h, pt-assign.cc, pt-assign.h, pt-binop.cc, pt-binop.h, pt-bp.cc, pt-bp.h, pt-cbinop.cc, pt-cbinop.h, pt-cell.cc, pt-cell.h, pt-check.cc, pt-check.h, pt-classdef.cc, pt-classdef.h, pt-cmd.h, pt-colon.cc, pt-colon.h, pt-const.cc, pt-const.h, pt-decl.cc, pt-decl.h, pt-eval.cc, pt-eval.h, pt-except.cc, pt-except.h, pt-exp.cc, pt-exp.h, pt-fcn-handle.cc, pt-fcn-handle.h, pt-id.cc, pt-id.h, pt-idx.cc, pt-idx.h, pt-jump.h, pt-loop.cc, pt-loop.h, pt-mat.cc, pt-mat.h, pt-misc.cc, pt-misc.h, pt-pr-code.cc, pt-pr-code.h, pt-select.cc, pt-select.h, pt-spmd.cc, pt-spmd.h, pt-stmt.cc, pt-stmt.h, pt-tm-const.cc, pt-tm-const.h, pt-unop.cc, pt-unop.h, pt-walk.cc, pt-walk.h, pt.cc, pt.h, token.cc, token.h, Range.cc, Range.h, idx-vector.cc, idx-vector.h, range-fwd.h, CollocWt.cc, CollocWt.h, aepbalance.cc, aepbalance.h, chol.cc, chol.h, gepbalance.cc, gepbalance.h, gsvd.cc, gsvd.h, hess.cc, hess.h, lo-mappers.cc, lo-mappers.h, lo-specfun.cc, lo-specfun.h, lu.cc, lu.h, oct-convn.cc, oct-convn.h, oct-fftw.cc, oct-fftw.h, oct-norm.cc, oct-norm.h, oct-rand.cc, oct-rand.h, oct-spparms.cc, oct-spparms.h, qr.cc, qr.h, qrp.cc, qrp.h, randgamma.cc, randgamma.h, randmtzig.cc, randmtzig.h, randpoisson.cc, randpoisson.h, schur.cc, schur.h, sparse-chol.cc, sparse-chol.h, sparse-lu.cc, sparse-lu.h, sparse-qr.cc, sparse-qr.h, svd.cc, svd.h, child-list.cc, child-list.h, dir-ops.cc, dir-ops.h, file-ops.cc, file-ops.h, file-stat.cc, file-stat.h, lo-sysdep.cc, lo-sysdep.h, lo-sysinfo.cc, lo-sysinfo.h, mach-info.cc, mach-info.h, oct-env.cc, oct-env.h, oct-group.cc, oct-group.h, oct-password.cc, oct-password.h, oct-syscalls.cc, oct-syscalls.h, oct-time.cc, oct-time.h, oct-uname.cc, oct-uname.h, action-container.cc, action-container.h, base-list.h, cmd-edit.cc, cmd-edit.h, cmd-hist.cc, cmd-hist.h, f77-fcn.h, file-info.cc, file-info.h, lo-array-errwarn.cc, lo-array-errwarn.h, lo-hash.cc, lo-hash.h, lo-ieee.h, lo-regexp.cc, lo-regexp.h, lo-utils.cc, lo-utils.h, oct-base64.cc, oct-base64.h, oct-glob.cc, oct-glob.h, oct-inttypes.h, oct-mutex.cc, oct-mutex.h, oct-refcount.h, oct-shlib.cc, oct-shlib.h, oct-sparse.cc, oct-sparse.h, oct-string.h, octave-preserve-stream-state.h, pathsearch.cc, pathsearch.h, quit.cc, quit.h, unwind-prot.cc, unwind-prot.h, url-transfer.cc, url-transfer.h: Re-indent code after switch to using namespace macros.
author Rik <rik@octave.org>
date Thu, 01 Dec 2022 18:02:15 -0800
parents e88a07dec498
children 597f3ee61a48
comparison
equal deleted inserted replaced
31605:e88a07dec498 31607:aac27ad79be6
42 #include "oct-locbuf.h" 42 #include "oct-locbuf.h"
43 #include "oct-shlib.h" 43 #include "oct-shlib.h"
44 44
45 OCTAVE_BEGIN_NAMESPACE(octave) 45 OCTAVE_BEGIN_NAMESPACE(octave)
46 46
47 static std::unordered_map<std::string, void *> gsvd_fcn; 47 static std::unordered_map<std::string, void *> gsvd_fcn;
48 48
49 static bool have_DGGSVD3 = false; 49 static bool have_DGGSVD3 = false;
50 static bool gsvd_initialized = false; 50 static bool gsvd_initialized = false;
51 51
52 /* Hack to stringize results of F77_FUNC macro. */ 52 /* Hack to stringize results of F77_FUNC macro. */
53 #define xSTRINGIZE(x) #x 53 #define xSTRINGIZE(x) #x
54 #define STRINGIZE(x) xSTRINGIZE(x) 54 #define STRINGIZE(x) xSTRINGIZE(x)
55 55
56 static void initialize_gsvd (void) 56 static void initialize_gsvd (void)
57 { 57 {
58 if (gsvd_initialized) 58 if (gsvd_initialized)
59 return; 59 return;
60 60
61 dynamic_library libs (""); 61 dynamic_library libs ("");
62 if (! libs) 62 if (! libs)
63 (*current_liboctave_error_handler) 63 (*current_liboctave_error_handler)
64 ("gsvd: unable to query LAPACK library"); 64 ("gsvd: unable to query LAPACK library");
65 65
66 have_DGGSVD3 = (libs.search (STRINGIZE (F77_FUNC (dggsvd3, DGGSVD3))) 66 have_DGGSVD3 = (libs.search (STRINGIZE (F77_FUNC (dggsvd3, DGGSVD3)))
67 != nullptr); 67 != nullptr);
68 68
69 if (have_DGGSVD3) 69 if (have_DGGSVD3)
70 { 70 {
71 gsvd_fcn["dg"] = libs.search (STRINGIZE (F77_FUNC (dggsvd3, DGGSVD3))); 71 gsvd_fcn["dg"] = libs.search (STRINGIZE (F77_FUNC (dggsvd3, DGGSVD3)));
72 gsvd_fcn["sg"] = libs.search (STRINGIZE (F77_FUNC (sggsvd3, SGGSVD3))); 72 gsvd_fcn["sg"] = libs.search (STRINGIZE (F77_FUNC (sggsvd3, SGGSVD3)));
73 gsvd_fcn["zg"] = libs.search (STRINGIZE (F77_FUNC (zggsvd3, ZGGSVD3))); 73 gsvd_fcn["zg"] = libs.search (STRINGIZE (F77_FUNC (zggsvd3, ZGGSVD3)));
74 gsvd_fcn["cg"] = libs.search (STRINGIZE (F77_FUNC (cggsvd3, CGGSVD3))); 74 gsvd_fcn["cg"] = libs.search (STRINGIZE (F77_FUNC (cggsvd3, CGGSVD3)));
75 } 75 }
76 else 76 else
77 { 77 {
78 gsvd_fcn["dg"] = libs.search (STRINGIZE (F77_FUNC (dggsvd, DGGSVD))); 78 gsvd_fcn["dg"] = libs.search (STRINGIZE (F77_FUNC (dggsvd, DGGSVD)));
79 gsvd_fcn["sg"] = libs.search (STRINGIZE (F77_FUNC (sggsvd, SGGSVD))); 79 gsvd_fcn["sg"] = libs.search (STRINGIZE (F77_FUNC (sggsvd, SGGSVD)));
80 gsvd_fcn["zg"] = libs.search (STRINGIZE (F77_FUNC (zggsvd, ZGGSVD))); 80 gsvd_fcn["zg"] = libs.search (STRINGIZE (F77_FUNC (zggsvd, ZGGSVD)));
81 gsvd_fcn["cg"] = libs.search (STRINGIZE (F77_FUNC (cggsvd, CGGSVD))); 81 gsvd_fcn["cg"] = libs.search (STRINGIZE (F77_FUNC (cggsvd, CGGSVD)));
82 } 82 }
83 83
84 gsvd_initialized = true; 84 gsvd_initialized = true;
85 } 85 }
86 86
87 /* Clean up macro namespace as soon as we are done using them */ 87 /* Clean up macro namespace as soon as we are done using them */
88 #undef xSTRINGIZE 88 #undef xSTRINGIZE
89 #undef STRINGIZE 89 #undef STRINGIZE
90 90
91 template<class T1> 91 template<class T1>
92 struct real_ggsvd_ptr 92 struct real_ggsvd_ptr
93 { 93 {
94 typedef F77_RET_T (*type) 94 typedef F77_RET_T (*type)
95 (F77_CONST_CHAR_ARG_DECL, // JOBU 95 (F77_CONST_CHAR_ARG_DECL, // JOBU
96 F77_CONST_CHAR_ARG_DECL, // JOBV 96 F77_CONST_CHAR_ARG_DECL, // JOBV
97 F77_CONST_CHAR_ARG_DECL, // JOBQ 97 F77_CONST_CHAR_ARG_DECL, // JOBQ
98 const F77_INT&, // M 98 const F77_INT&, // M
99 const F77_INT&, // N 99 const F77_INT&, // N
100 const F77_INT&, // P 100 const F77_INT&, // P
101 F77_INT&, // K 101 F77_INT&, // K
102 F77_INT&, // L 102 F77_INT&, // L
103 T1 *, // A(LDA,N) 103 T1 *, // A(LDA,N)
104 const F77_INT&, // LDA 104 const F77_INT&, // LDA
105 T1 *, // B(LDB,N) 105 T1 *, // B(LDB,N)
106 const F77_INT&, // LDB 106 const F77_INT&, // LDB
107 T1 *, // ALPHA(N) 107 T1 *, // ALPHA(N)
108 T1 *, // BETA(N) 108 T1 *, // BETA(N)
109 T1 *, // U(LDU,M) 109 T1 *, // U(LDU,M)
110 const F77_INT&, // LDU 110 const F77_INT&, // LDU
111 T1 *, // V(LDV,P) 111 T1 *, // V(LDV,P)
112 const F77_INT&, // LDV 112 const F77_INT&, // LDV
113 T1 *, // Q(LDQ,N) 113 T1 *, // Q(LDQ,N)
114 const F77_INT&, // LDQ 114 const F77_INT&, // LDQ
115 T1 *, // WORK 115 T1 *, // WORK
116 F77_INT *, // IWORK(N) 116 F77_INT *, // IWORK(N)
117 F77_INT& // INFO 117 F77_INT& // INFO
118 F77_CHAR_ARG_LEN_DECL 118 F77_CHAR_ARG_LEN_DECL
119 F77_CHAR_ARG_LEN_DECL 119 F77_CHAR_ARG_LEN_DECL
120 F77_CHAR_ARG_LEN_DECL); 120 F77_CHAR_ARG_LEN_DECL);
121 }; 121 };
122 122
123 template<class T1> 123 template<class T1>
124 struct real_ggsvd3_ptr 124 struct real_ggsvd3_ptr
125 { 125 {
126 typedef F77_RET_T (*type) 126 typedef F77_RET_T (*type)
127 (F77_CONST_CHAR_ARG_DECL, // JOBU 127 (F77_CONST_CHAR_ARG_DECL, // JOBU
128 F77_CONST_CHAR_ARG_DECL, // JOBV 128 F77_CONST_CHAR_ARG_DECL, // JOBV
129 F77_CONST_CHAR_ARG_DECL, // JOBQ 129 F77_CONST_CHAR_ARG_DECL, // JOBQ
130 const F77_INT&, // M 130 const F77_INT&, // M
131 const F77_INT&, // N 131 const F77_INT&, // N
132 const F77_INT&, // P 132 const F77_INT&, // P
133 F77_INT&, // K 133 F77_INT&, // K
134 F77_INT&, // L 134 F77_INT&, // L
135 T1 *, // A(LDA,N) 135 T1 *, // A(LDA,N)
136 const F77_INT&, // LDA 136 const F77_INT&, // LDA
137 T1 *, // B(LDB,N) 137 T1 *, // B(LDB,N)
138 const F77_INT&, // LDB 138 const F77_INT&, // LDB
139 T1 *, // ALPHA(N) 139 T1 *, // ALPHA(N)
140 T1 *, // BETA(N) 140 T1 *, // BETA(N)
141 T1 *, // U(LDU,M) 141 T1 *, // U(LDU,M)
142 const F77_INT&, // LDU 142 const F77_INT&, // LDU
143 T1 *, // V(LDV,P) 143 T1 *, // V(LDV,P)
144 const F77_INT&, // LDV 144 const F77_INT&, // LDV
145 T1 *, // Q(LDQ,N) 145 T1 *, // Q(LDQ,N)
146 const F77_INT&, // LDQ 146 const F77_INT&, // LDQ
147 T1 *, // WORK 147 T1 *, // WORK
148 const F77_INT&, // LWORK 148 const F77_INT&, // LWORK
149 F77_INT *, // IWORK(N) 149 F77_INT *, // IWORK(N)
150 F77_INT& // INFO 150 F77_INT& // INFO
151 F77_CHAR_ARG_LEN_DECL 151 F77_CHAR_ARG_LEN_DECL
152 F77_CHAR_ARG_LEN_DECL 152 F77_CHAR_ARG_LEN_DECL
153 F77_CHAR_ARG_LEN_DECL); 153 F77_CHAR_ARG_LEN_DECL);
154 }; 154 };
155 155
156 template<class T1, class T2> 156 template<class T1, class T2>
157 struct comp_ggsvd_ptr 157 struct comp_ggsvd_ptr
158 { 158 {
159 typedef F77_RET_T (*type) 159 typedef F77_RET_T (*type)
160 (F77_CONST_CHAR_ARG_DECL, // JOBU 160 (F77_CONST_CHAR_ARG_DECL, // JOBU
161 F77_CONST_CHAR_ARG_DECL, // JOBV 161 F77_CONST_CHAR_ARG_DECL, // JOBV
162 F77_CONST_CHAR_ARG_DECL, // JOBQ 162 F77_CONST_CHAR_ARG_DECL, // JOBQ
163 const F77_INT&, // M 163 const F77_INT&, // M
164 const F77_INT&, // N 164 const F77_INT&, // N
165 const F77_INT&, // P 165 const F77_INT&, // P
166 F77_INT&, // K 166 F77_INT&, // K
167 F77_INT&, // L 167 F77_INT&, // L
168 T1 *, // A(LDA,N) 168 T1 *, // A(LDA,N)
169 const F77_INT&, // LDA 169 const F77_INT&, // LDA
170 T1 *, // B(LDB,N) 170 T1 *, // B(LDB,N)
171 const F77_INT&, // LDB 171 const F77_INT&, // LDB
172 T2 *, // ALPHA(N) 172 T2 *, // ALPHA(N)
173 T2 *, // BETA(N) 173 T2 *, // BETA(N)
174 T1 *, // U(LDU,M) 174 T1 *, // U(LDU,M)
175 const F77_INT&, // LDU 175 const F77_INT&, // LDU
176 T1 *, // V(LDV,P) 176 T1 *, // V(LDV,P)
177 const F77_INT&, // LDV 177 const F77_INT&, // LDV
178 T1 *, // Q(LDQ,N) 178 T1 *, // Q(LDQ,N)
179 const F77_INT&, // LDQ 179 const F77_INT&, // LDQ
180 T1 *, // WORK 180 T1 *, // WORK
181 T2 *, // RWORK 181 T2 *, // RWORK
182 F77_INT *, // IWORK(N) 182 F77_INT *, // IWORK(N)
183 F77_INT& // INFO 183 F77_INT& // INFO
184 F77_CHAR_ARG_LEN_DECL 184 F77_CHAR_ARG_LEN_DECL
185 F77_CHAR_ARG_LEN_DECL 185 F77_CHAR_ARG_LEN_DECL
186 F77_CHAR_ARG_LEN_DECL); 186 F77_CHAR_ARG_LEN_DECL);
187 }; 187 };
188 188
189 template<class T1, class T2> 189 template<class T1, class T2>
190 struct comp_ggsvd3_ptr 190 struct comp_ggsvd3_ptr
191 { 191 {
192 typedef F77_RET_T (*type) 192 typedef F77_RET_T (*type)
193 (F77_CONST_CHAR_ARG_DECL, // JOBU 193 (F77_CONST_CHAR_ARG_DECL, // JOBU
194 F77_CONST_CHAR_ARG_DECL, // JOBV 194 F77_CONST_CHAR_ARG_DECL, // JOBV
195 F77_CONST_CHAR_ARG_DECL, // JOBQ 195 F77_CONST_CHAR_ARG_DECL, // JOBQ
196 const F77_INT&, // M 196 const F77_INT&, // M
197 const F77_INT&, // N 197 const F77_INT&, // N
198 const F77_INT&, // P 198 const F77_INT&, // P
199 F77_INT&, // K 199 F77_INT&, // K
200 F77_INT&, // L 200 F77_INT&, // L
201 T1 *, // A(LDA,N) 201 T1 *, // A(LDA,N)
202 const F77_INT&, // LDA 202 const F77_INT&, // LDA
203 T1 *, // B(LDB,N) 203 T1 *, // B(LDB,N)
204 const F77_INT&, // LDB 204 const F77_INT&, // LDB
205 T2 *, // ALPHA(N) 205 T2 *, // ALPHA(N)
206 T2 *, // BETA(N) 206 T2 *, // BETA(N)
207 T1 *, // U(LDU,M) 207 T1 *, // U(LDU,M)
208 const F77_INT&, // LDU 208 const F77_INT&, // LDU
209 T1 *, // V(LDV,P) 209 T1 *, // V(LDV,P)
210 const F77_INT&, // LDV 210 const F77_INT&, // LDV
211 T1 *, // Q(LDQ,N) 211 T1 *, // Q(LDQ,N)
212 const F77_INT&, // LDQ 212 const F77_INT&, // LDQ
213 T1 *, // WORK 213 T1 *, // WORK
214 const F77_INT&, // LWORK 214 const F77_INT&, // LWORK
215 T2 *, // RWORK 215 T2 *, // RWORK
216 F77_INT *, // IWORK(N) 216 F77_INT *, // IWORK(N)
217 F77_INT& // INFO 217 F77_INT& // INFO
218 F77_CHAR_ARG_LEN_DECL 218 F77_CHAR_ARG_LEN_DECL
219 F77_CHAR_ARG_LEN_DECL 219 F77_CHAR_ARG_LEN_DECL
220 F77_CHAR_ARG_LEN_DECL); 220 F77_CHAR_ARG_LEN_DECL);
221 }; 221 };
222 222
223 // template specializations 223 // template specializations
224 typedef real_ggsvd_ptr<F77_DBLE>::type dggsvd_type; 224 typedef real_ggsvd_ptr<F77_DBLE>::type dggsvd_type;
225 typedef real_ggsvd3_ptr<F77_DBLE>::type dggsvd3_type; 225 typedef real_ggsvd3_ptr<F77_DBLE>::type dggsvd3_type;
226 typedef real_ggsvd_ptr<F77_REAL>::type sggsvd_type; 226 typedef real_ggsvd_ptr<F77_REAL>::type sggsvd_type;
227 typedef real_ggsvd3_ptr<F77_REAL>::type sggsvd3_type; 227 typedef real_ggsvd3_ptr<F77_REAL>::type sggsvd3_type;
228 typedef comp_ggsvd_ptr<F77_DBLE_CMPLX, F77_DBLE>::type zggsvd_type; 228 typedef comp_ggsvd_ptr<F77_DBLE_CMPLX, F77_DBLE>::type zggsvd_type;
229 typedef comp_ggsvd3_ptr<F77_DBLE_CMPLX, F77_DBLE>::type zggsvd3_type; 229 typedef comp_ggsvd3_ptr<F77_DBLE_CMPLX, F77_DBLE>::type zggsvd3_type;
230 typedef comp_ggsvd_ptr<F77_CMPLX, F77_REAL>::type cggsvd_type; 230 typedef comp_ggsvd_ptr<F77_CMPLX, F77_REAL>::type cggsvd_type;
231 typedef comp_ggsvd3_ptr<F77_CMPLX, F77_REAL>::type cggsvd3_type; 231 typedef comp_ggsvd3_ptr<F77_CMPLX, F77_REAL>::type cggsvd3_type;
232 232
233 OCTAVE_BEGIN_NAMESPACE(math) 233 OCTAVE_BEGIN_NAMESPACE(math)
234 234
235 template <> 235 template <>
236 void 236 void
237 gsvd<Matrix>::ggsvd (char& jobu, char& jobv, char& jobq, F77_INT m, 237 gsvd<Matrix>::ggsvd (char& jobu, char& jobv, char& jobq, F77_INT m,
238 F77_INT n, F77_INT p, F77_INT& k, F77_INT& l, 238 F77_INT n, F77_INT p, F77_INT& k, F77_INT& l,
239 double *tmp_dataA, F77_INT m1, double *tmp_dataB, 239 double *tmp_dataA, F77_INT m1, double *tmp_dataB,
240 F77_INT p1, Matrix& alpha, Matrix& beta, double *u, 240 F77_INT p1, Matrix& alpha, Matrix& beta, double *u,
241 F77_INT nrow_u, double *v, F77_INT nrow_v, double *q, 241 F77_INT nrow_u, double *v, F77_INT nrow_v, double *q,
242 F77_INT nrow_q, double *work, F77_INT lwork, 242 F77_INT nrow_q, double *work, F77_INT lwork,
243 F77_INT *iwork, F77_INT& info) 243 F77_INT *iwork, F77_INT& info)
244 { 244 {
245 if (! gsvd_initialized) 245 if (! gsvd_initialized)
246 initialize_gsvd (); 246 initialize_gsvd ();
247 247
248 if (have_DGGSVD3) 248 if (have_DGGSVD3)
249 {
250 dggsvd3_type f_ptr = reinterpret_cast<dggsvd3_type> (gsvd_fcn["dg"]);
251 f_ptr (F77_CONST_CHAR_ARG2 (&jobu, 1),
252 F77_CONST_CHAR_ARG2 (&jobv, 1),
253 F77_CONST_CHAR_ARG2 (&jobq, 1),
254 m, n, p, k, l, tmp_dataA, m1, tmp_dataB, p1,
255 alpha.fortran_vec (), beta.fortran_vec (),
256 u, nrow_u, v, nrow_v, q, nrow_q,
257 work, lwork, iwork, info
258 F77_CHAR_ARG_LEN (1)
259 F77_CHAR_ARG_LEN (1)
260 F77_CHAR_ARG_LEN (1));
261 }
262 else
263 {
264 dggsvd_type f_ptr = reinterpret_cast<dggsvd_type> (gsvd_fcn["dg"]);
265 f_ptr (F77_CONST_CHAR_ARG2 (&jobu, 1),
266 F77_CONST_CHAR_ARG2 (&jobv, 1),
267 F77_CONST_CHAR_ARG2 (&jobq, 1),
268 m, n, p, k, l, tmp_dataA, m1, tmp_dataB, p1,
269 alpha.fortran_vec (), beta.fortran_vec (),
270 u, nrow_u, v, nrow_v, q, nrow_q,
271 work, iwork, info
272 F77_CHAR_ARG_LEN (1)
273 F77_CHAR_ARG_LEN (1)
274 F77_CHAR_ARG_LEN (1));
275 }
276 }
277
278 template <>
279 void
280 gsvd<FloatMatrix>::ggsvd (char& jobu, char& jobv, char& jobq, F77_INT m,
281 F77_INT n, F77_INT p, F77_INT& k, F77_INT& l,
282 float *tmp_dataA, F77_INT m1, float *tmp_dataB,
283 F77_INT p1, FloatMatrix& alpha,
284 FloatMatrix& beta, float *u, F77_INT nrow_u,
285 float *v, F77_INT nrow_v, float *q,
286 F77_INT nrow_q, float *work,
287 F77_INT lwork, F77_INT *iwork, F77_INT& info)
288 {
289 if (! gsvd_initialized)
290 initialize_gsvd ();
291
292 if (have_DGGSVD3)
293 {
294 sggsvd3_type f_ptr = reinterpret_cast<sggsvd3_type> (gsvd_fcn["sg"]);
295 f_ptr (F77_CONST_CHAR_ARG2 (&jobu, 1),
296 F77_CONST_CHAR_ARG2 (&jobv, 1),
297 F77_CONST_CHAR_ARG2 (&jobq, 1),
298 m, n, p, k, l, tmp_dataA, m1, tmp_dataB, p1,
299 alpha.fortran_vec (), beta.fortran_vec (),
300 u, nrow_u, v, nrow_v, q, nrow_q,
301 work, lwork, iwork, info
302 F77_CHAR_ARG_LEN (1)
303 F77_CHAR_ARG_LEN (1)
304 F77_CHAR_ARG_LEN (1));
305 }
306 else
307 {
308 sggsvd_type f_ptr = reinterpret_cast<sggsvd_type> (gsvd_fcn["sg"]);
309 f_ptr (F77_CONST_CHAR_ARG2 (&jobu, 1),
310 F77_CONST_CHAR_ARG2 (&jobv, 1),
311 F77_CONST_CHAR_ARG2 (&jobq, 1),
312 m, n, p, k, l, tmp_dataA, m1, tmp_dataB, p1,
313 alpha.fortran_vec (), beta.fortran_vec (),
314 u, nrow_u, v, nrow_v, q, nrow_q,
315 work, iwork, info
316 F77_CHAR_ARG_LEN (1)
317 F77_CHAR_ARG_LEN (1)
318 F77_CHAR_ARG_LEN (1));
319 }
320 }
321
322 template <>
323 void
324 gsvd<ComplexMatrix>::ggsvd (char& jobu, char& jobv, char& jobq,
325 F77_INT m, F77_INT n, F77_INT p, F77_INT& k,
326 F77_INT& l, Complex *tmp_dataA, F77_INT m1,
327 Complex *tmp_dataB, F77_INT p1, Matrix& alpha,
328 Matrix& beta, Complex *u, F77_INT nrow_u,
329 Complex *v, F77_INT nrow_v, Complex *q,
330 F77_INT nrow_q, Complex *work,
331 F77_INT lwork, F77_INT *iwork, F77_INT& info)
332 {
333 if (! gsvd_initialized)
334 initialize_gsvd ();
335
336 OCTAVE_LOCAL_BUFFER(double, rwork, 2*n);
337
338 if (have_DGGSVD3)
339 {
340 zggsvd3_type f_ptr = reinterpret_cast<zggsvd3_type> (gsvd_fcn["zg"]);
341 f_ptr (F77_CONST_CHAR_ARG2 (&jobu, 1),
342 F77_CONST_CHAR_ARG2 (&jobv, 1),
343 F77_CONST_CHAR_ARG2 (&jobq, 1),
344 m, n, p, k, l,
345 F77_DBLE_CMPLX_ARG (tmp_dataA), m1,
346 F77_DBLE_CMPLX_ARG (tmp_dataB), p1,
347 alpha.fortran_vec (), beta.fortran_vec (),
348 F77_DBLE_CMPLX_ARG (u), nrow_u,
349 F77_DBLE_CMPLX_ARG (v), nrow_v,
350 F77_DBLE_CMPLX_ARG (q), nrow_q,
351 F77_DBLE_CMPLX_ARG (work),
352 lwork, rwork, iwork, info
353 F77_CHAR_ARG_LEN (1)
354 F77_CHAR_ARG_LEN (1)
355 F77_CHAR_ARG_LEN (1));
356 }
357 else
358 {
359 zggsvd_type f_ptr = reinterpret_cast<zggsvd_type> (gsvd_fcn["zg"]);
360 f_ptr (F77_CONST_CHAR_ARG2 (&jobu, 1),
361 F77_CONST_CHAR_ARG2 (&jobv, 1),
362 F77_CONST_CHAR_ARG2 (&jobq, 1),
363 m, n, p, k, l,
364 F77_DBLE_CMPLX_ARG (tmp_dataA), m1,
365 F77_DBLE_CMPLX_ARG (tmp_dataB), p1,
366 alpha.fortran_vec (), beta.fortran_vec (),
367 F77_DBLE_CMPLX_ARG (u), nrow_u,
368 F77_DBLE_CMPLX_ARG (v), nrow_v,
369 F77_DBLE_CMPLX_ARG (q), nrow_q,
370 F77_DBLE_CMPLX_ARG (work),
371 rwork, iwork, info
372 F77_CHAR_ARG_LEN (1)
373 F77_CHAR_ARG_LEN (1)
374 F77_CHAR_ARG_LEN (1));
375 }
376 }
377
378 template <>
379 void
380 gsvd<FloatComplexMatrix>::ggsvd (char& jobu, char& jobv, char& jobq,
381 F77_INT m, F77_INT n, F77_INT p,
382 F77_INT& k, F77_INT& l,
383 FloatComplex *tmp_dataA, F77_INT m1,
384 FloatComplex *tmp_dataB, F77_INT p1,
385 FloatMatrix& alpha, FloatMatrix& beta,
386 FloatComplex *u, F77_INT nrow_u,
387 FloatComplex *v, F77_INT nrow_v,
388 FloatComplex *q, F77_INT nrow_q,
389 FloatComplex *work, F77_INT lwork,
390 F77_INT *iwork, F77_INT& info)
391 {
392 if (! gsvd_initialized)
393 initialize_gsvd ();
394
395 OCTAVE_LOCAL_BUFFER(float, rwork, 2*n);
396
397 if (have_DGGSVD3)
398 {
399 cggsvd3_type f_ptr = reinterpret_cast<cggsvd3_type> (gsvd_fcn["cg"]);
400 f_ptr (F77_CONST_CHAR_ARG2 (&jobu, 1),
401 F77_CONST_CHAR_ARG2 (&jobv, 1),
402 F77_CONST_CHAR_ARG2 (&jobq, 1),
403 m, n, p, k, l,
404 F77_CMPLX_ARG (tmp_dataA), m1,
405 F77_CMPLX_ARG (tmp_dataB), p1,
406 alpha.fortran_vec (), beta.fortran_vec (),
407 F77_CMPLX_ARG (u), nrow_u,
408 F77_CMPLX_ARG (v), nrow_v,
409 F77_CMPLX_ARG (q), nrow_q,
410 F77_CMPLX_ARG (work), lwork,
411 rwork, iwork, info
412 F77_CHAR_ARG_LEN (1)
413 F77_CHAR_ARG_LEN (1)
414 F77_CHAR_ARG_LEN (1));
415 }
416 else
417 {
418 cggsvd_type f_ptr = reinterpret_cast<cggsvd_type> (gsvd_fcn["cg"]);
419 f_ptr (F77_CONST_CHAR_ARG2 (&jobu, 1),
420 F77_CONST_CHAR_ARG2 (&jobv, 1),
421 F77_CONST_CHAR_ARG2 (&jobq, 1),
422 m, n, p, k, l,
423 F77_CMPLX_ARG (tmp_dataA), m1,
424 F77_CMPLX_ARG (tmp_dataB), p1,
425 alpha.fortran_vec (), beta.fortran_vec (),
426 F77_CMPLX_ARG (u), nrow_u,
427 F77_CMPLX_ARG (v), nrow_v,
428 F77_CMPLX_ARG (q), nrow_q,
429 F77_CMPLX_ARG (work),
430 rwork, iwork, info
431 F77_CHAR_ARG_LEN (1)
432 F77_CHAR_ARG_LEN (1)
433 F77_CHAR_ARG_LEN (1));
434 }
435 }
436
437 template <typename T>
438 T
439 gsvd<T>::left_singular_matrix_A (void) const
440 {
441 if (m_type == gsvd::Type::sigma_only)
442 (*current_liboctave_error_handler)
443 ("gsvd: U not computed because type == gsvd::sigma_only");
444
445 return m_left_smA;
446 }
447
448 template <typename T>
449 T
450 gsvd<T>::left_singular_matrix_B (void) const
451 {
452 if (m_type == gsvd::Type::sigma_only)
453 (*current_liboctave_error_handler)
454 ("gsvd: V not computed because type == gsvd::sigma_only");
455
456 return m_left_smB;
457 }
458
459 template <typename T>
460 T
461 gsvd<T>::right_singular_matrix (void) const
462 {
463 if (m_type == gsvd::Type::sigma_only)
464 (*current_liboctave_error_handler)
465 ("gsvd: X not computed because type == gsvd::sigma_only");
466
467 return m_right_sm;
468 }
469
470 template <typename T>
471 gsvd<T>::gsvd (const T& a, const T& b, gsvd::Type gsvd_type)
472 {
473 if (a.isempty () || b.isempty ())
474 (*current_liboctave_error_handler)
475 ("gsvd: A and B cannot be empty matrices");
476
477 F77_INT info;
478
479 F77_INT m = to_f77_int (a.rows ());
480 F77_INT n = to_f77_int (a.cols ());
481 F77_INT p = to_f77_int (b.rows ());
482
483 T atmp = a;
484 P *tmp_dataA = atmp.fortran_vec ();
485
486 T btmp = b;
487 P *tmp_dataB = btmp.fortran_vec ();
488
489 char jobu = 'U';
490 char jobv = 'V';
491 char jobq = 'Q';
492
493 F77_INT nrow_u = m;
494 F77_INT nrow_v = p;
495 F77_INT nrow_q = n;
496
497 F77_INT k, l;
498
499 switch (gsvd_type)
500 {
501 case gsvd<T>::Type::sigma_only:
502
503 // FIXME: In LAPACK 3.0, problem below seems to be fixed,
504 // so we now set jobX = 'N'.
505 //
506 // For calculating sigma_only, both jobu and jobv should be 'N', but
507 // there seems to be a bug in dgesvd from LAPACK V2.0. To
508 // demonstrate the bug, set both jobu and jobv to 'N' and find the
509 // singular values of [eye(3), eye(3)]. The result is
510 // [-sqrt(2), -sqrt(2), -sqrt(2)].
511
512 jobu = jobv = jobq = 'N';
513 nrow_u = nrow_v = nrow_q = 1;
514 break;
515
516 default:
517 break;
518 }
519
520 m_type = gsvd_type;
521
522 if (jobu != 'N')
523 m_left_smA.resize (nrow_u, m);
524
525 P *u = m_left_smA.fortran_vec ();
526
527 if (jobv != 'N')
528 m_left_smB.resize (nrow_v, p);
529
530 P *v = m_left_smB.fortran_vec ();
531
532 if (jobq != 'N')
533 m_right_sm.resize (nrow_q, n);
534
535 P *q = m_right_sm.fortran_vec ();
536
537 real_matrix alpha (n, 1);
538 real_matrix beta (n, 1);
539
540 OCTAVE_LOCAL_BUFFER(F77_INT, iwork, n);
541
542 if (! gsvd_initialized)
543 initialize_gsvd ();
544
545 F77_INT lwork;
546 if (have_DGGSVD3)
547 {
548 lwork = -1;
549 P work_tmp;
550
551 gsvd<T>::ggsvd (jobu, jobv, jobq, m, n, p, k, l,
552 tmp_dataA, m, tmp_dataB, p,
553 alpha, beta, u, nrow_u, v, nrow_v, q, nrow_q,
554 &work_tmp, lwork, iwork, info);
555
556 lwork = static_cast<F77_INT> (std::abs (work_tmp));
557 }
558 else
559 {
560 lwork = std::max ({3*n, m, p}) + n;
561 }
562 info = 0;
563
564 OCTAVE_LOCAL_BUFFER(P, work, lwork);
565
566 gsvd<T>::ggsvd (jobu, jobv, jobq, m, n, p, k, l,
567 tmp_dataA, m, tmp_dataB, p,
568 alpha, beta, u, nrow_u, v, nrow_v, q, nrow_q,
569 work, lwork, iwork, info);
570
571 if (info < 0)
572 (*current_liboctave_error_handler)
573 ("*ggsvd.f: argument %" OCTAVE_F77_INT_TYPE_FORMAT " illegal",
574 -info);
575
576 if (info > 0)
577 (*current_liboctave_error_handler)
578 ("*ggsvd.f: Jacobi-type procedure failed to converge");
579
580 F77_INT i, j;
581
582 if (gsvd_type != gsvd<T>::Type::sigma_only)
583 {
584 // Size according to LAPACK is k+l,k+l, but this needs
585 // to be nxn for Matlab compatibility.
586 T R (n, n, 0.0);
587 int astart = n-k-l;
588 if (m - k - l >= 0)
249 { 589 {
250 dggsvd3_type f_ptr = reinterpret_cast<dggsvd3_type> (gsvd_fcn["dg"]); 590 // R is stored in A(1:K+L,N-K-L+1:N)
251 f_ptr (F77_CONST_CHAR_ARG2 (&jobu, 1), 591 for (i = 0; i < k+l; i++)
252 F77_CONST_CHAR_ARG2 (&jobv, 1), 592 for (j = 0; j < k+l; j++)
253 F77_CONST_CHAR_ARG2 (&jobq, 1), 593 R.xelem (i, j) = atmp.xelem (i, astart + j);
254 m, n, p, k, l, tmp_dataA, m1, tmp_dataB, p1,
255 alpha.fortran_vec (), beta.fortran_vec (),
256 u, nrow_u, v, nrow_v, q, nrow_q,
257 work, lwork, iwork, info
258 F77_CHAR_ARG_LEN (1)
259 F77_CHAR_ARG_LEN (1)
260 F77_CHAR_ARG_LEN (1));
261 } 594 }
262 else 595 else
263 { 596 {
264 dggsvd_type f_ptr = reinterpret_cast<dggsvd_type> (gsvd_fcn["dg"]); 597 // (R11 R12 R13 ) is stored in A(1:M, N-K-L+1:N)
265 f_ptr (F77_CONST_CHAR_ARG2 (&jobu, 1), 598 // ( 0 R22 R23 )
266 F77_CONST_CHAR_ARG2 (&jobv, 1), 599 for (i = 0; i < m; i++)
267 F77_CONST_CHAR_ARG2 (&jobq, 1), 600 for (j = 0; j < k+l; j++)
268 m, n, p, k, l, tmp_dataA, m1, tmp_dataB, p1, 601 R.xelem (i, j) = atmp.xelem (i, astart + j);
269 alpha.fortran_vec (), beta.fortran_vec (), 602 // and R33 is stored in B(M-K+1:L,N+M-K-L+1:N)
270 u, nrow_u, v, nrow_v, q, nrow_q, 603 for (i = m; i < k + l; i++)
271 work, iwork, info 604 for (j = n - l - k + m; j < n; j++)
272 F77_CHAR_ARG_LEN (1) 605 R.xelem (i, j) = btmp.xelem (i - k, j);
273 F77_CHAR_ARG_LEN (1)
274 F77_CHAR_ARG_LEN (1));
275 } 606 }
276 } 607
277 608 // Output X = Q*R'
278 template <> 609 // FIXME: Is there a way to call BLAS multiply directly
279 void 610 // with flags so that R is transposed?
280 gsvd<FloatMatrix>::ggsvd (char& jobu, char& jobv, char& jobq, F77_INT m, 611 m_right_sm = m_right_sm * R.hermitian ();
281 F77_INT n, F77_INT p, F77_INT& k, F77_INT& l, 612 }
282 float *tmp_dataA, F77_INT m1, float *tmp_dataB, 613
283 F77_INT p1, FloatMatrix& alpha, 614 // Fill in C and S
284 FloatMatrix& beta, float *u, F77_INT nrow_u, 615 F77_INT rank;
285 float *v, F77_INT nrow_v, float *q, 616 bool fill_ptn;
286 F77_INT nrow_q, float *work, 617 if (m-k-l >= 0)
287 F77_INT lwork, F77_INT *iwork, F77_INT& info) 618 {
288 { 619 rank = l;
289 if (! gsvd_initialized) 620 fill_ptn = true;
290 initialize_gsvd (); 621 }
291 622 else
292 if (have_DGGSVD3) 623 {
624 rank = m-k;
625 fill_ptn = false;
626 }
627
628 if (gsvd_type == gsvd<T>::Type::sigma_only)
629 {
630 // Return column vector with results
631 m_sigmaA.resize (k+l, 1);
632 m_sigmaB.resize (k+l, 1);
633
634 if (fill_ptn)
293 { 635 {
294 sggsvd3_type f_ptr = reinterpret_cast<sggsvd3_type> (gsvd_fcn["sg"]); 636 for (i = 0; i < k; i++)
295 f_ptr (F77_CONST_CHAR_ARG2 (&jobu, 1), 637 {
296 F77_CONST_CHAR_ARG2 (&jobv, 1), 638 m_sigmaA.xelem (i) = 1.0;
297 F77_CONST_CHAR_ARG2 (&jobq, 1), 639 m_sigmaB.xelem (i) = 0.0;
298 m, n, p, k, l, tmp_dataA, m1, tmp_dataB, p1, 640 }
299 alpha.fortran_vec (), beta.fortran_vec (), 641 for (i = k, j = k+l-1; i < k+l; i++, j--)
300 u, nrow_u, v, nrow_v, q, nrow_q, 642 {
301 work, lwork, iwork, info 643 m_sigmaA.xelem (i) = alpha.xelem (i);
302 F77_CHAR_ARG_LEN (1) 644 m_sigmaB.xelem (i) = beta.xelem (i);
303 F77_CHAR_ARG_LEN (1) 645 }
304 F77_CHAR_ARG_LEN (1));
305 } 646 }
306 else 647 else
307 { 648 {
308 sggsvd_type f_ptr = reinterpret_cast<sggsvd_type> (gsvd_fcn["sg"]); 649 for (i = 0; i < k; i++)
309 f_ptr (F77_CONST_CHAR_ARG2 (&jobu, 1),
310 F77_CONST_CHAR_ARG2 (&jobv, 1),
311 F77_CONST_CHAR_ARG2 (&jobq, 1),
312 m, n, p, k, l, tmp_dataA, m1, tmp_dataB, p1,
313 alpha.fortran_vec (), beta.fortran_vec (),
314 u, nrow_u, v, nrow_v, q, nrow_q,
315 work, iwork, info
316 F77_CHAR_ARG_LEN (1)
317 F77_CHAR_ARG_LEN (1)
318 F77_CHAR_ARG_LEN (1));
319 }
320 }
321
322 template <>
323 void
324 gsvd<ComplexMatrix>::ggsvd (char& jobu, char& jobv, char& jobq,
325 F77_INT m, F77_INT n, F77_INT p, F77_INT& k,
326 F77_INT& l, Complex *tmp_dataA, F77_INT m1,
327 Complex *tmp_dataB, F77_INT p1, Matrix& alpha,
328 Matrix& beta, Complex *u, F77_INT nrow_u,
329 Complex *v, F77_INT nrow_v, Complex *q,
330 F77_INT nrow_q, Complex *work,
331 F77_INT lwork, F77_INT *iwork, F77_INT& info)
332 {
333 if (! gsvd_initialized)
334 initialize_gsvd ();
335
336 OCTAVE_LOCAL_BUFFER(double, rwork, 2*n);
337
338 if (have_DGGSVD3)
339 {
340 zggsvd3_type f_ptr = reinterpret_cast<zggsvd3_type> (gsvd_fcn["zg"]);
341 f_ptr (F77_CONST_CHAR_ARG2 (&jobu, 1),
342 F77_CONST_CHAR_ARG2 (&jobv, 1),
343 F77_CONST_CHAR_ARG2 (&jobq, 1),
344 m, n, p, k, l,
345 F77_DBLE_CMPLX_ARG (tmp_dataA), m1,
346 F77_DBLE_CMPLX_ARG (tmp_dataB), p1,
347 alpha.fortran_vec (), beta.fortran_vec (),
348 F77_DBLE_CMPLX_ARG (u), nrow_u,
349 F77_DBLE_CMPLX_ARG (v), nrow_v,
350 F77_DBLE_CMPLX_ARG (q), nrow_q,
351 F77_DBLE_CMPLX_ARG (work),
352 lwork, rwork, iwork, info
353 F77_CHAR_ARG_LEN (1)
354 F77_CHAR_ARG_LEN (1)
355 F77_CHAR_ARG_LEN (1));
356 }
357 else
358 {
359 zggsvd_type f_ptr = reinterpret_cast<zggsvd_type> (gsvd_fcn["zg"]);
360 f_ptr (F77_CONST_CHAR_ARG2 (&jobu, 1),
361 F77_CONST_CHAR_ARG2 (&jobv, 1),
362 F77_CONST_CHAR_ARG2 (&jobq, 1),
363 m, n, p, k, l,
364 F77_DBLE_CMPLX_ARG (tmp_dataA), m1,
365 F77_DBLE_CMPLX_ARG (tmp_dataB), p1,
366 alpha.fortran_vec (), beta.fortran_vec (),
367 F77_DBLE_CMPLX_ARG (u), nrow_u,
368 F77_DBLE_CMPLX_ARG (v), nrow_v,
369 F77_DBLE_CMPLX_ARG (q), nrow_q,
370 F77_DBLE_CMPLX_ARG (work),
371 rwork, iwork, info
372 F77_CHAR_ARG_LEN (1)
373 F77_CHAR_ARG_LEN (1)
374 F77_CHAR_ARG_LEN (1));
375 }
376 }
377
378 template <>
379 void
380 gsvd<FloatComplexMatrix>::ggsvd (char& jobu, char& jobv, char& jobq,
381 F77_INT m, F77_INT n, F77_INT p,
382 F77_INT& k, F77_INT& l,
383 FloatComplex *tmp_dataA, F77_INT m1,
384 FloatComplex *tmp_dataB, F77_INT p1,
385 FloatMatrix& alpha, FloatMatrix& beta,
386 FloatComplex *u, F77_INT nrow_u,
387 FloatComplex *v, F77_INT nrow_v,
388 FloatComplex *q, F77_INT nrow_q,
389 FloatComplex *work, F77_INT lwork,
390 F77_INT *iwork, F77_INT& info)
391 {
392 if (! gsvd_initialized)
393 initialize_gsvd ();
394
395 OCTAVE_LOCAL_BUFFER(float, rwork, 2*n);
396
397 if (have_DGGSVD3)
398 {
399 cggsvd3_type f_ptr = reinterpret_cast<cggsvd3_type> (gsvd_fcn["cg"]);
400 f_ptr (F77_CONST_CHAR_ARG2 (&jobu, 1),
401 F77_CONST_CHAR_ARG2 (&jobv, 1),
402 F77_CONST_CHAR_ARG2 (&jobq, 1),
403 m, n, p, k, l,
404 F77_CMPLX_ARG (tmp_dataA), m1,
405 F77_CMPLX_ARG (tmp_dataB), p1,
406 alpha.fortran_vec (), beta.fortran_vec (),
407 F77_CMPLX_ARG (u), nrow_u,
408 F77_CMPLX_ARG (v), nrow_v,
409 F77_CMPLX_ARG (q), nrow_q,
410 F77_CMPLX_ARG (work), lwork,
411 rwork, iwork, info
412 F77_CHAR_ARG_LEN (1)
413 F77_CHAR_ARG_LEN (1)
414 F77_CHAR_ARG_LEN (1));
415 }
416 else
417 {
418 cggsvd_type f_ptr = reinterpret_cast<cggsvd_type> (gsvd_fcn["cg"]);
419 f_ptr (F77_CONST_CHAR_ARG2 (&jobu, 1),
420 F77_CONST_CHAR_ARG2 (&jobv, 1),
421 F77_CONST_CHAR_ARG2 (&jobq, 1),
422 m, n, p, k, l,
423 F77_CMPLX_ARG (tmp_dataA), m1,
424 F77_CMPLX_ARG (tmp_dataB), p1,
425 alpha.fortran_vec (), beta.fortran_vec (),
426 F77_CMPLX_ARG (u), nrow_u,
427 F77_CMPLX_ARG (v), nrow_v,
428 F77_CMPLX_ARG (q), nrow_q,
429 F77_CMPLX_ARG (work),
430 rwork, iwork, info
431 F77_CHAR_ARG_LEN (1)
432 F77_CHAR_ARG_LEN (1)
433 F77_CHAR_ARG_LEN (1));
434 }
435 }
436
437 template <typename T>
438 T
439 gsvd<T>::left_singular_matrix_A (void) const
440 {
441 if (m_type == gsvd::Type::sigma_only)
442 (*current_liboctave_error_handler)
443 ("gsvd: U not computed because type == gsvd::sigma_only");
444
445 return m_left_smA;
446 }
447
448 template <typename T>
449 T
450 gsvd<T>::left_singular_matrix_B (void) const
451 {
452 if (m_type == gsvd::Type::sigma_only)
453 (*current_liboctave_error_handler)
454 ("gsvd: V not computed because type == gsvd::sigma_only");
455
456 return m_left_smB;
457 }
458
459 template <typename T>
460 T
461 gsvd<T>::right_singular_matrix (void) const
462 {
463 if (m_type == gsvd::Type::sigma_only)
464 (*current_liboctave_error_handler)
465 ("gsvd: X not computed because type == gsvd::sigma_only");
466
467 return m_right_sm;
468 }
469
470 template <typename T>
471 gsvd<T>::gsvd (const T& a, const T& b, gsvd::Type gsvd_type)
472 {
473 if (a.isempty () || b.isempty ())
474 (*current_liboctave_error_handler)
475 ("gsvd: A and B cannot be empty matrices");
476
477 F77_INT info;
478
479 F77_INT m = to_f77_int (a.rows ());
480 F77_INT n = to_f77_int (a.cols ());
481 F77_INT p = to_f77_int (b.rows ());
482
483 T atmp = a;
484 P *tmp_dataA = atmp.fortran_vec ();
485
486 T btmp = b;
487 P *tmp_dataB = btmp.fortran_vec ();
488
489 char jobu = 'U';
490 char jobv = 'V';
491 char jobq = 'Q';
492
493 F77_INT nrow_u = m;
494 F77_INT nrow_v = p;
495 F77_INT nrow_q = n;
496
497 F77_INT k, l;
498
499 switch (gsvd_type)
500 {
501 case gsvd<T>::Type::sigma_only:
502
503 // FIXME: In LAPACK 3.0, problem below seems to be fixed,
504 // so we now set jobX = 'N'.
505 //
506 // For calculating sigma_only, both jobu and jobv should be 'N', but
507 // there seems to be a bug in dgesvd from LAPACK V2.0. To
508 // demonstrate the bug, set both jobu and jobv to 'N' and find the
509 // singular values of [eye(3), eye(3)]. The result is
510 // [-sqrt(2), -sqrt(2), -sqrt(2)].
511
512 jobu = jobv = jobq = 'N';
513 nrow_u = nrow_v = nrow_q = 1;
514 break;
515
516 default:
517 break;
518 }
519
520 m_type = gsvd_type;
521
522 if (jobu != 'N')
523 m_left_smA.resize (nrow_u, m);
524
525 P *u = m_left_smA.fortran_vec ();
526
527 if (jobv != 'N')
528 m_left_smB.resize (nrow_v, p);
529
530 P *v = m_left_smB.fortran_vec ();
531
532 if (jobq != 'N')
533 m_right_sm.resize (nrow_q, n);
534
535 P *q = m_right_sm.fortran_vec ();
536
537 real_matrix alpha (n, 1);
538 real_matrix beta (n, 1);
539
540 OCTAVE_LOCAL_BUFFER(F77_INT, iwork, n);
541
542 if (! gsvd_initialized)
543 initialize_gsvd ();
544
545 F77_INT lwork;
546 if (have_DGGSVD3)
547 {
548 lwork = -1;
549 P work_tmp;
550
551 gsvd<T>::ggsvd (jobu, jobv, jobq, m, n, p, k, l,
552 tmp_dataA, m, tmp_dataB, p,
553 alpha, beta, u, nrow_u, v, nrow_v, q, nrow_q,
554 &work_tmp, lwork, iwork, info);
555
556 lwork = static_cast<F77_INT> (std::abs (work_tmp));
557 }
558 else
559 {
560 lwork = std::max ({3*n, m, p}) + n;
561 }
562 info = 0;
563
564 OCTAVE_LOCAL_BUFFER(P, work, lwork);
565
566 gsvd<T>::ggsvd (jobu, jobv, jobq, m, n, p, k, l,
567 tmp_dataA, m, tmp_dataB, p,
568 alpha, beta, u, nrow_u, v, nrow_v, q, nrow_q,
569 work, lwork, iwork, info);
570
571 if (info < 0)
572 (*current_liboctave_error_handler)
573 ("*ggsvd.f: argument %" OCTAVE_F77_INT_TYPE_FORMAT " illegal",
574 -info);
575
576 if (info > 0)
577 (*current_liboctave_error_handler)
578 ("*ggsvd.f: Jacobi-type procedure failed to converge");
579
580 F77_INT i, j;
581
582 if (gsvd_type != gsvd<T>::Type::sigma_only)
583 {
584 // Size according to LAPACK is k+l,k+l, but this needs
585 // to be nxn for Matlab compatibility.
586 T R (n, n, 0.0);
587 int astart = n-k-l;
588 if (m - k - l >= 0)
589 { 650 {
590 // R is stored in A(1:K+L,N-K-L+1:N) 651 m_sigmaA.xelem (i) = 1.0;
591 for (i = 0; i < k+l; i++) 652 m_sigmaB.xelem (i) = 0.0;
592 for (j = 0; j < k+l; j++)
593 R.xelem (i, j) = atmp.xelem (i, astart + j);
594 } 653 }
595 else 654 for (i = k; i < m; i++)
596 { 655 {
597 // (R11 R12 R13 ) is stored in A(1:M, N-K-L+1:N) 656 m_sigmaA.xelem (i) = alpha.xelem (i);
598 // ( 0 R22 R23 ) 657 m_sigmaB.xelem (i) = beta.xelem (i);
599 for (i = 0; i < m; i++)
600 for (j = 0; j < k+l; j++)
601 R.xelem (i, j) = atmp.xelem (i, astart + j);
602 // and R33 is stored in B(M-K+1:L,N+M-K-L+1:N)
603 for (i = m; i < k + l; i++)
604 for (j = n - l - k + m; j < n; j++)
605 R.xelem (i, j) = btmp.xelem (i - k, j);
606 } 658 }
607 659 for (i = m; i < k+l; i++)
608 // Output X = Q*R'
609 // FIXME: Is there a way to call BLAS multiply directly
610 // with flags so that R is transposed?
611 m_right_sm = m_right_sm * R.hermitian ();
612 }
613
614 // Fill in C and S
615 F77_INT rank;
616 bool fill_ptn;
617 if (m-k-l >= 0)
618 {
619 rank = l;
620 fill_ptn = true;
621 }
622 else
623 {
624 rank = m-k;
625 fill_ptn = false;
626 }
627
628 if (gsvd_type == gsvd<T>::Type::sigma_only)
629 {
630 // Return column vector with results
631 m_sigmaA.resize (k+l, 1);
632 m_sigmaB.resize (k+l, 1);
633
634 if (fill_ptn)
635 { 660 {
636 for (i = 0; i < k; i++) 661 m_sigmaA.xelem (i) = 0.0;
637 { 662 m_sigmaB.xelem (i) = 1.0;
638 m_sigmaA.xelem (i) = 1.0;
639 m_sigmaB.xelem (i) = 0.0;
640 }
641 for (i = k, j = k+l-1; i < k+l; i++, j--)
642 {
643 m_sigmaA.xelem (i) = alpha.xelem (i);
644 m_sigmaB.xelem (i) = beta.xelem (i);
645 }
646 }
647 else
648 {
649 for (i = 0; i < k; i++)
650 {
651 m_sigmaA.xelem (i) = 1.0;
652 m_sigmaB.xelem (i) = 0.0;
653 }
654 for (i = k; i < m; i++)
655 {
656 m_sigmaA.xelem (i) = alpha.xelem (i);
657 m_sigmaB.xelem (i) = beta.xelem (i);
658 }
659 for (i = m; i < k+l; i++)
660 {
661 m_sigmaA.xelem (i) = 0.0;
662 m_sigmaB.xelem (i) = 1.0;
663 }
664 } 663 }
665 } 664 }
666 else // returning all matrices 665 }
666 else // returning all matrices
667 {
668 // Number of columns according to LAPACK is k+l, but this needs
669 // to be n for Matlab compatibility.
670 m_sigmaA.resize (m, n);
671 m_sigmaB.resize (p, n);
672
673 for (i = 0; i < k; i++)
674 m_sigmaA.xelem (i, i) = 1.0;
675
676 for (i = 0; i < rank; i++)
667 { 677 {
668 // Number of columns according to LAPACK is k+l, but this needs 678 m_sigmaA.xelem (k+i, k+i) = alpha.xelem (k+i);
669 // to be n for Matlab compatibility. 679 m_sigmaB.xelem (i, k+i) = beta.xelem (k+i);
670 m_sigmaA.resize (m, n);
671 m_sigmaB.resize (p, n);
672
673 for (i = 0; i < k; i++)
674 m_sigmaA.xelem (i, i) = 1.0;
675
676 for (i = 0; i < rank; i++)
677 {
678 m_sigmaA.xelem (k+i, k+i) = alpha.xelem (k+i);
679 m_sigmaB.xelem (i, k+i) = beta.xelem (k+i);
680 }
681
682 if (! fill_ptn)
683 {
684 for (i = m; i < n; i++)
685 m_sigmaB.xelem (i-k, i) = 1.0;
686 }
687
688 } 680 }
689 } 681
690 682 if (! fill_ptn)
691 // Instantiations needed in octave::math namespace. 683 {
692 template class gsvd<Matrix>; 684 for (i = m; i < n; i++)
693 template class gsvd<FloatMatrix>; 685 m_sigmaB.xelem (i-k, i) = 1.0;
694 template class gsvd<ComplexMatrix>; 686 }
695 template class gsvd<FloatComplexMatrix>; 687
688 }
689 }
690
691 // Instantiations needed in octave::math namespace.
692 template class gsvd<Matrix>;
693 template class gsvd<FloatMatrix>;
694 template class gsvd<ComplexMatrix>;
695 template class gsvd<FloatComplexMatrix>;
696 696
697 OCTAVE_END_NAMESPACE(math) 697 OCTAVE_END_NAMESPACE(math)
698 OCTAVE_END_NAMESPACE(octave) 698 OCTAVE_END_NAMESPACE(octave)