Mercurial > octave

comparison libinterp/corefcn/gsvd.cc @ 30300:4ee01c14fccd

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Rewrite gsvd function and return *correct* 3rd output (bug #60273). Correctly use the API interface to LAPACK which Octave got wrong for certain matrix combinations. Perform post-processing on LAPACK results to return Matlab-compatible outputs. * libinterp/corefcn/gsvd.cc (gsvd_type): Add second input (nargin) to be able to distinguish the type of gsvd to generate. * libinterp/corefcn/gsvd.cc (do_gsvd): Add another input (nargin) to be able to calculate the correct gsvd type. Update code for API change in liboctave which now returns a full Matrix for singular values rather than a DiagMatrix. Sort singular values in ascending order for compatibility with Matlab. * libinterp/corefcn/gsvd.cc (Fgsvd): Update documation to warn about rank-deficient input case. Add new input validation to check that number of columns of A and B match. Move special check for empty matrices from libinterp to liboctave. Update all BIST tests for new behavior of gsvd code in liboctave. * liboctave/numeric/gsvd.h (gsvd): Remove member variable R from gsvd class. * liboctave/numeric/gsvd.h (gsvd::singular_values_A, gsvd_singular_values_B): Change return type to T::real_matrix_type from T::real_diag_matrix_type. * liboctave/numeric/gsvd.h (gsvd::ggsvd): Change default gsvd::Type to "std" from "economy". * liboctave/numeric/gsvd.cc: Add #includes for <algorithm>, <undordered_map>, and "oct-locbuf.h". Change gsvd_fcn from std::map to std::unordered_map. * liboctave/numeric/gsvd.cc (initialize_gsvd): Throw an error if Octave is unable to connect to LAPACK library. * liboctave/numeric/gsvd.cc (ggsvd): Change function prototype to not use references to Octave Matrix datatypes, but rather point directly to underlying plain old datatype. For example, use "double *" rather than "Matrix&". Replace scratch memory instances created with Matrix() with OCTAVE_LOCAL_BUFFER which relies on C++ STL. * liboctave/numeric/gsvd.cc (left_singular_matrix_A, left_singular_matrix_B, right_singular_matrix): Use the fact that current_liboctave_error_handler never returns to eliminate else branch of if/else code. * liboctave/numeric/gsvd.cc (R_matrix): Delete unused function. * liboctave/numeric/gsvd.cc (gsvd): Add input validation for empty matrices. Simplify 'job' values processing which seems to have been fixed in LAPACK 3.0. Replace scratch memory created with std::vector with OCTAVE_LOCAL_BUFFER. Replace scratch memory created with Octave Matrix classes with OCTAVE_LOCAL_BUFFER. Use initializer list to std::max to simplify code. Re-code interface to LAPACK to extract 'R' matrix. Post-process to calculate third output of gsvd as X = Q*R'. Correct algorithm to extract singular values from LAPACK..
author Rik <rik@octave.org>
date Wed, 08 Sep 2021 16:08:03 -0700
parents 7d6709900da7
children cd63a97cb9be
comparison
equal deleted inserted replaced
30299:f306fe9bfa0d 30300:4ee01c14fccd
42 42
43 OCTAVE_NAMESPACE_BEGIN 43 OCTAVE_NAMESPACE_BEGIN
44 44
45 template <typename T> 45 template <typename T>
46 static typename math::gsvd<T>::Type 46 static typename math::gsvd<T>::Type
47 gsvd_type (int nargout) 47 gsvd_type (int nargout, int nargin)
48 { 48 {
49 return ((nargout == 0 || nargout == 1) 49 if (nargout == 0 || nargout == 1)
50 ? math::gsvd<T>::Type::sigma_only 50 return octave::math::gsvd<T>::Type::sigma_only;
51 : (nargout > 5) ? math::gsvd<T>::Type::std 51 else if (nargin < 3)
52 : math::gsvd<T>::Type::economy); 52 return octave::math::gsvd<T>::Type::std;
53 else
54 return octave::math::gsvd<T>::Type::economy;
53 } 55 }
54 56
55 // Named like this to avoid conflicts with the gsvd class. 57 // Named do_gsvd to avoid conflicts with the gsvd class itself.
56 template <typename T> 58 template <typename T>
57 static octave_value_list 59 static octave_value_list
58 do_gsvd (const T& A, const T& B, const octave_idx_type nargout, 60 do_gsvd (const T& A, const T& B,
61 const octave_idx_type nargout, const octave_idx_type nargin,
59 bool is_single = false) 62 bool is_single = false)
60 { 63 {
61 math::gsvd<T> result (A, B, gsvd_type<T> (nargout)); 64 math::gsvd<T> result (A, B, gsvd_type<T> (nargout, nargin));
62 65
63 octave_value_list retval (nargout); 66 octave_value_list retval (nargout);
64 if (nargout < 2) 67 if (nargout <= 1)
65 { 68 {
66 if (is_single) 69 if (is_single)
67 { 70 {
68 FloatDiagMatrix sigA = result.singular_values_A (); 71 FloatMatrix sigA = result.singular_values_A ();
69 FloatDiagMatrix sigB = result.singular_values_B (); 72 FloatMatrix sigB = result.singular_values_B ();
70 for (int i = sigA.rows () - 1; i >= 0; i--) 73 for (int i = sigA.rows () - 1; i >= 0; i--)
71 sigA.dgxelem(i) /= sigB.dgxelem(i); 74 sigA.xelem (i) /= sigB.xelem (i);
72 retval(0) = sigA.diag (); 75 retval(0) = sigA.sort ();
73 } 76 }
74 else 77 else
75 { 78 {
76 DiagMatrix sigA = result.singular_values_A (); 79 Matrix sigA = result.singular_values_A ();
77 DiagMatrix sigB = result.singular_values_B (); 80 Matrix sigB = result.singular_values_B ();
78 for (int i = sigA.rows () - 1; i >= 0; i--) 81 for (int i = sigA.rows () - 1; i >= 0; i--)
79 sigA.dgxelem(i) /= sigB.dgxelem(i); 82 sigA.xelem (i) /= sigB.xelem (i);
80 retval(0) = sigA.diag (); 83 retval(0) = sigA.sort ();
81 } 84 }
82 } 85 }
83 else 86 else
84 { 87 {
88 switch (nargout)
89 {
90 case 5:
91 retval(4) = result.singular_values_B ();
92 OCTAVE_FALLTHROUGH;
93
94 case 4:
95 retval(3) = result.singular_values_A ();
96 OCTAVE_FALLTHROUGH;
97
98 case 3:
99 retval(2) = result.right_singular_matrix ();
100 OCTAVE_FALLTHROUGH;
101 }
102 retval(1) = result.left_singular_matrix_B ();
85 retval(0) = result.left_singular_matrix_A (); 103 retval(0) = result.left_singular_matrix_A ();
86 retval(1) = result.left_singular_matrix_B ();
87 if (nargout > 2)
88 retval(2) = result.right_singular_matrix ();
89 if (nargout > 3)
90 retval(3) = result.singular_values_A ();
91 if (nargout > 4)
92 retval(4) = result.singular_values_B ();
93 if (nargout > 5)
94 retval(5) = result.R_matrix ();
95 } 104 }
105
96 return retval; 106 return retval;
97 } 107 }
98 108
99 DEFUN (gsvd, args, nargout, 109 DEFUN (gsvd, args, nargout,
100 doc: /* -*- texinfo -*- 110 doc: /* -*- texinfo -*-
143 "economy-sized" decomposition where the number of columns of @var{U}, @var{V} 153 "economy-sized" decomposition where the number of columns of @var{U}, @var{V}
144 and the number of rows of @var{C}, @var{S} is less than or equal to the number 154 and the number of rows of @var{C}, @var{S} is less than or equal to the number
145 of columns of @var{A}. This option is not yet implemented. 155 of columns of @var{A}. This option is not yet implemented.
146 156
147 Programming Note: the code is a wrapper to the corresponding @sc{lapack} dggsvd 157 Programming Note: the code is a wrapper to the corresponding @sc{lapack} dggsvd
148 and zggsvd routines. 158 and zggsvd routines. If matrices @var{A} and @var{B} are @emph{both} rank
149 159 deficient then @sc{lapack} will return an incorrect factorization. Programmers
160 should avoid this combination.
150 @seealso{svd} 161 @seealso{svd}
151 @end deftypefn */) 162 @end deftypefn */)
152 { 163 {
153 int nargin = args.length (); 164 int nargin = args.length ();
154 165
155 if (nargin < 2 || nargin > 3) 166 if (nargin < 2 || nargin > 3)
156 print_usage (); 167 print_usage ();
157 else if (nargin == 3) 168 else if (nargin == 3)
158 warning ("gsvd: economy-sized decomposition is not yet implemented, returning full decomposition"); 169 {
170 // FIXME: when "economy" is implemented delete this code
171 warning ("gsvd: economy-sized decomposition is not yet implemented, returning full decomposition");
172 nargin = 2;
173 }
159 174
160 octave_value_list retval; 175 octave_value_list retval;
161 176
162 octave_value argA = args(0); 177 octave_value argA = args(0);
163 octave_value argB = args(1); 178 octave_value argB = args(1);
164 179
165 octave_idx_type nr = argA.rows (); 180 if (argA.columns () != argB.columns ())
166 octave_idx_type nc = argA.columns (); 181 error ("gsvd: A and B must have the same number of columns");
167 182
168 octave_idx_type np = argB.columns (); 183 if (argA.is_single_type () || argB.is_single_type ())
169
170 // FIXME: This "special" case should be handled in the gsvd class, not here
171 if (nr == 0 || nc == 0)
172 { 184 {
173 retval = octave_value_list (nargout); 185 if (argA.isreal () && argB.isreal ())
174 if (nargout < 2) // S = gsvd (A, B) 186 {
175 { 187 FloatMatrix tmpA = argA.xfloat_matrix_value ("gsvd: A must be a real or complex matrix");
176 if (argA.is_single_type () || argB.is_single_type ()) 188 FloatMatrix tmpB = argB.xfloat_matrix_value ("gsvd: B must be a real or complex matrix");
177 retval(0) = FloatMatrix (0, 1); 189
178 else 190 if (tmpA.any_element_is_inf_or_nan ())
179 retval(0) = Matrix (0, 1); 191 error ("gsvd: A cannot have Inf or NaN values");
180 } 192 if (tmpB.any_element_is_inf_or_nan ())
181 else // [U, V, X, C, S, R] = gsvd (A, B) 193 error ("gsvd: B cannot have Inf or NaN values");
182 { 194
183 if (argA.is_single_type () || argB.is_single_type ()) 195 retval = do_gsvd (tmpA, tmpB, nargout, nargin, true);
184 { 196 }
185 retval(0) = float_identity_matrix (nc, nc); 197 else if (argA.iscomplex () || argB.iscomplex ())
186 retval(1) = float_identity_matrix (nc, nc); 198 {
187 if (nargout > 2) 199 FloatComplexMatrix ctmpA = argA.xfloat_complex_matrix_value ("gsvd: A must be a real or complex matrix");
188 retval(2) = float_identity_matrix (nr, nr); 200 FloatComplexMatrix ctmpB = argB.xfloat_complex_matrix_value ("gsvd: B must be a real or complex matrix");
189 if (nargout > 3) 201
190 retval(3) = FloatMatrix (nr, nc); 202 if (ctmpA.any_element_is_inf_or_nan ())
191 if (nargout > 4) 203 error ("gsvd: A cannot have Inf or NaN values");
192 retval(4) = float_identity_matrix (nr, nr); 204 if (ctmpB.any_element_is_inf_or_nan ())
193 if (nargout > 5) 205 error ("gsvd: B cannot have Inf or NaN values");
194 retval(5) = float_identity_matrix (nr, nr); 206
195 } 207 retval = do_gsvd (ctmpA, ctmpB, nargout, nargin, true);
196 else 208 }
197 { 209 else
198 retval(0) = identity_matrix (nc, nc); 210 error ("gsvd: A and B must be real or complex matrices");
199 retval(1) = identity_matrix (nc, nc);
200 if (nargout > 2)
201 retval(2) = identity_matrix (nr, nr);
202 if (nargout > 3)
203 retval(3) = Matrix (nr, nc);
204 if (nargout > 4)
205 retval(4) = identity_matrix (nr, nr);
206 if (nargout > 5)
207 retval(5) = identity_matrix (nr, nr);
208 }
209 }
210 } 211 }
211 else 212 else
212 { 213 {
213 if (nc != np) 214 if (argA.isreal () && argB.isreal ())
214 print_usage (); 215 {
215 216 Matrix tmpA = argA.xmatrix_value ("gsvd: A must be a real or complex matrix");
216 if (argA.is_single_type () || argB.is_single_type ()) 217 Matrix tmpB = argB.xmatrix_value ("gsvd: B must be a real or complex matrix");
217 { 218
218 if (argA.isreal () && argB.isreal ()) 219 if (tmpA.any_element_is_inf_or_nan ())
219 { 220 error ("gsvd: A cannot have Inf or NaN values");
220 FloatMatrix tmpA = argA.xfloat_matrix_value ("gsvd: A must be a real or complex matrix"); 221 if (tmpB.any_element_is_inf_or_nan ())
221 FloatMatrix tmpB = argB.xfloat_matrix_value ("gsvd: B must be a real or complex matrix"); 222 error ("gsvd: B cannot have Inf or NaN values");
222 223
223 if (tmpA.any_element_is_inf_or_nan ()) 224 retval = do_gsvd (tmpA, tmpB, nargout, nargin);
224 error ("gsvd: A cannot have Inf or NaN values"); 225 }
225 if (tmpB.any_element_is_inf_or_nan ()) 226 else if (argA.iscomplex () || argB.iscomplex ())
226 error ("gsvd: B cannot have Inf or NaN values"); 227 {
227 228 ComplexMatrix ctmpA = argA.xcomplex_matrix_value ("gsvd: A must be a real or complex matrix");
228 retval = do_gsvd (tmpA, tmpB, nargout, true); 229 ComplexMatrix ctmpB = argB.xcomplex_matrix_value ("gsvd: B must be a real or complex matrix");
229 } 230
230 else if (argA.iscomplex () || argB.iscomplex ()) 231 if (ctmpA.any_element_is_inf_or_nan ())
231 { 232 error ("gsvd: A cannot have Inf or NaN values");
232 FloatComplexMatrix ctmpA = argA.xfloat_complex_matrix_value ("gsvd: A must be a real or complex matrix"); 233 if (ctmpB.any_element_is_inf_or_nan ())
233 FloatComplexMatrix ctmpB = argB.xfloat_complex_matrix_value ("gsvd: B must be a real or complex matrix"); 234 error ("gsvd: B cannot have Inf or NaN values");
234 235
235 if (ctmpA.any_element_is_inf_or_nan ()) 236 retval = do_gsvd (ctmpA, ctmpB, nargout, nargin);
236 error ("gsvd: A cannot have Inf or NaN values");
237 if (ctmpB.any_element_is_inf_or_nan ())
238 error ("gsvd: B cannot have Inf or NaN values");
239
240 retval = do_gsvd (ctmpA, ctmpB, nargout, true);
241 }
242 else
243 error ("gsvd: A and B must be real or complex matrices");
244 } 237 }
245 else 238 else
246 { 239 error ("gsvd: A and B must be real or complex matrices");
247 if (argA.isreal () && argB.isreal ())
248 {
249 Matrix tmpA = argA.xmatrix_value ("gsvd: A must be a real or complex matrix");
250 Matrix tmpB = argB.xmatrix_value ("gsvd: B must be a real or complex matrix");
251
252 if (tmpA.any_element_is_inf_or_nan ())
253 error ("gsvd: A cannot have Inf or NaN values");
254 if (tmpB.any_element_is_inf_or_nan ())
255 error ("gsvd: B cannot have Inf or NaN values");
256
257 retval = do_gsvd (tmpA, tmpB, nargout);
258 }
259 else if (argA.iscomplex () || argB.iscomplex ())
260 {
261 ComplexMatrix ctmpA = argA.xcomplex_matrix_value ("gsvd: A must be a real or complex matrix");
262 ComplexMatrix ctmpB = argB.xcomplex_matrix_value ("gsvd: B must be a real or complex matrix");
263
264 if (ctmpA.any_element_is_inf_or_nan ())
265 error ("gsvd: A cannot have Inf or NaN values");
266 if (ctmpB.any_element_is_inf_or_nan ())
267 error ("gsvd: B cannot have Inf or NaN values");
268
269 retval = do_gsvd (ctmpA, ctmpB, nargout);
270 }
271 else
272 error ("gsvd: A and B must be real or complex matrices");
273 }
274 } 240 }
275 241
276 return retval; 242 return retval;
277 } 243 }
278 244
279 /* 245 /*
280 246
281 ## Basic test of decomposition 247 ## Basic tests of decomposition
282 %!test <48807> 248 %!test <60273>
283 %! A = reshape (1:15,5,3); 249 %! A = reshape (1:15,5,3);
284 %! B = magic (3); 250 %! B = magic (3);
285 %! [U,V,X,C,S] = gsvd (A,B); 251 %! [U,V,X,C,S] = gsvd (A,B);
252 %! assert (size (U), [5, 5]);
253 %! assert (size (V), [3, 3]);
254 %! assert (size (X), [3, 3]);
255 %! assert (size (C), [5, 3]);
256 %! assert (C(4:5, :), zeros (2,3));
257 %! assert (size (S), [3, 3]);
286 %! assert (U*C*X', A, 50*eps); 258 %! assert (U*C*X', A, 50*eps);
287 %! assert (V*S*X', B, 50*eps); 259 %! assert (V*S*X', B, 50*eps);
288 %! S0 = gsvd (A, B); 260 %! S0 = gsvd (A, B);
289 %! S1 = svd (A / B); 261 %! assert (size (S0), [3, 1]);
262 %! S1 = sort (svd (A / B));
290 %! assert (S0, S1, 10*eps); 263 %! assert (S0, S1, 10*eps);
291 264
265 %!test <60273>
266 %! A = reshape (1:15,3,5);
267 %! B = magic (5);
268 %! [U,V,X,C,S] = gsvd (A,B);
269 %! assert (size (U), [3, 3]);
270 %! assert (size (V), [5, 5]);
271 %! assert (size (X), [5, 5]);
272 %! assert (size (C), [3, 5]);
273 %! assert (C(:, 4:5), zeros (3,2));
274 %! assert (size (S), [5, 5]);
275 %! assert (U*C*X', A, 100*eps); # less accurate in this orientation
276 %! assert (V*S*X', B, 125*eps); # for some reason.
277 %! S0 = gsvd (A, B);
278 %! assert (size (S0), [5, 1]);
279 %! S0 = S0(3:end);
280 %! S1 = sort (svd (A / B));
281 %! assert (S0, S1, 20*eps);
282
292 ## a few tests for gsvd.m 283 ## a few tests for gsvd.m
293 %!shared A, A0, B, B0, U, V, C, S, X, R, D1, D2 284 %!shared A, A0, B, B0, U, V, C, S, X, old_state, restore_state
285 %! old_state = randn ("state");
286 %! restore_state = onCleanup (@() randn ("state", old_state));
287 %! randn ("state", 40); # initialize generator to make behavior reproducible
294 %! A0 = randn (5, 3); 288 %! A0 = randn (5, 3);
295 %! B0 = diag ([1 2 4]); 289 %! B0 = diag ([1 2 4]);
296 %! A = A0; 290 %! A = A0;
297 %! B = B0; 291 %! B = B0;
298 292
299 ## A (5x3) and B (3x3) are full rank 293 ## A (5x3) and B (3x3) are full rank
300 %!test <48807> 294 %!test <48807>
301 %! [U, V, X, C, S, R] = gsvd (A, B); 295 %! [U, V, X, C, S] = gsvd (A, B);
302 %! D1 = zeros (5, 3); D1(1:3, 1:3) = C; 296 %! assert (C'*C + S'*S, eye (3), 5*eps);
303 %! D2 = S; 297 %! assert (U*C*X', A, 10*eps);
304 %! assert (norm (diag (C).^2 + diag (S).^2 - ones (3, 1)) <= 1e-6); 298 %! assert (V*S*X', B, 20*eps);
305 %! assert (norm ((U'*A*X) - D1*R) <= 1e-6);
306 %! assert (norm ((V'*B*X) - D2*R) <= 1e-6);
307 299
308 ## A: 5x3 full rank, B: 3x3 rank deficient 300 ## A: 5x3 full rank, B: 3x3 rank deficient
309 %!test <48807> 301 %!test <48807>
310 %! B(2, 2) = 0; 302 %! B(2, 2) = 0;
311 %! [U, V, X, C, S, R] = gsvd (A, B); 303 %! [U, V, X, C, S] = gsvd (A, B);
312 %! D1 = zeros (5, 3); D1(1, 1) = 1; D1(2:3, 2:3) = C; 304 %! assert (C'*C + S'*S, eye (3), 5*eps);
313 %! D2 = [zeros(2, 1) S; zeros(1, 3)]; 305 %! assert (U*C*X', A, 10*eps);
314 %! assert (norm (diag (C).^2 + diag (S).^2 - ones (2, 1)) <= 1e-6); 306 %! assert (V*S*X', B, 20*eps);
315 %! assert (norm ((U'*A*X) - D1*R) <= 1e-6);
316 %! assert (norm ((V'*B*X) - D2*R) <= 1e-6);
317 307
318 ## A: 5x3 rank deficient, B: 3x3 full rank 308 ## A: 5x3 rank deficient, B: 3x3 full rank
319 %!test <48807> 309 %!test <48807>
320 %! B = B0; 310 %! B = B0;
321 %! A(:, 3) = 2*A(:, 1) - A(:, 2); 311 %! A(:, 3) = 2*A(:, 1) - A(:, 2);
322 %! [U, V, X, C, S, R] = gsvd (A, B); 312 %! [U, V, X, C, S] = gsvd (A, B);
323 %! D1 = zeros (5, 3); D1(1:3, 1:3) = C; 313 %! assert (C'*C + S'*S, eye (3), 5*eps);
324 %! D2 = S; 314 %! assert (U*C*X', A, 10*eps);
325 %! assert (norm (diag (C).^2 + diag (S).^2 - ones (3, 1)) <= 1e-6); 315 %! assert (V*S*X', B, 20*eps);
326 %! assert (norm ((U'*A*X) - D1*R) <= 1e-6);
327 %! assert (norm ((V'*B*X) - D2*R) <= 1e-6);
328 316
329 ## A and B are both rank deficient 317 ## A and B are both rank deficient
330 %!test <48807> 318 ## FIXME: LAPACK seems to be completely broken for this case
319 %!#test <48807>
331 %! B(:, 3) = 2*B(:, 1) - B(:, 2); 320 %! B(:, 3) = 2*B(:, 1) - B(:, 2);
332 %! [U, V, X, C, S, R] = gsvd (A, B); 321 %! [U, V, X, C, S] = gsvd (A, B);
333 %! D1 = zeros (5, 2); D1(1:2, 1:2) = C; 322 %! assert (C'*C + S'*S, eye (3), 5*eps);
334 %! D2 = [S; zeros(1, 2)]; 323 %! assert (U*C*X', A, 10*eps);
335 %! assert (norm (diag (C).^2 + diag (S).^2 - ones (2, 1)) <= 1e-6); 324 %! assert (V*S*X', B, 20*eps);
336 %! assert (norm ((U'*A*X) - D1*[zeros(2, 1) R]) <= 1e-6);
337 %! assert (norm ((V'*B*X) - D2*[zeros(2, 1) R]) <= 1e-6);
338 325
339 ## A (now 3x5) and B (now 5x5) are full rank 326 ## A (now 3x5) and B (now 5x5) are full rank
340 %!test <48807> 327 %!test <48807>
341 %! A = A0.'; 328 %! A = A0.';
342 %! B0 = diag ([1 2 4 8 16]); 329 %! B0 = diag ([1 2 4 8 16]);
343 %! B = B0; 330 %! B = B0;
344 %! [U, V, X, C, S, R] = gsvd (A, B); 331 %! [U, V, X, C, S] = gsvd (A, B);
345 %! D1 = [C zeros(3,2)]; 332 %! assert (C'*C + S'*S, eye (5), 5*eps);
346 %! D2 = [S zeros(3,2); zeros(2, 3) eye(2)]; 333 %! assert (U*C*X', A, 10*eps);
347 %! assert (norm (diag (C).^2 + diag (S).^2 - ones (3, 1)) <= 1e-6); 334 %! assert (V*S*X', B, 40*eps);
348 %! assert (norm ((U'*A*X) - D1*R) <= 1e-6);
349 %! assert (norm ((V'*B*X) - D2*R) <= 1e-6);
350 335
351 ## A: 3x5 full rank, B: 5x5 rank deficient 336 ## A: 3x5 full rank, B: 5x5 rank deficient
352 %!test <48807> 337 %!test <48807>
353 %! B(2, 2) = 0; 338 %! B(2, 2) = 0;
354 %! [U, V, X, C, S, R] = gsvd (A, B); 339 %! [U, V, X, C, S] = gsvd (A, B);
355 %! D1 = zeros (3, 5); D1(1, 1) = 1; D1(2:3, 2:3) = C; 340 %! assert (C'*C + S'*S, eye (5), 5*eps);
356 %! D2 = zeros (5, 5); D2(1:2, 2:3) = S; D2(3:4, 4:5) = eye (2); 341 %! assert (U*C*X', A, 10*eps);
357 %! assert (norm (diag (C).^2 + diag (S).^2 - ones (2, 1)) <= 1e-6); 342 %! assert (V*S*X', B, 40*eps);
358 %! assert (norm ((U'*A*X) - D1*R) <= 1e-6);
359 %! assert (norm ((V'*B*X) - D2*R) <= 1e-6);
360 343
361 ## A: 3x5 rank deficient, B: 5x5 full rank 344 ## A: 3x5 rank deficient, B: 5x5 full rank
362 %!test <48807> 345 %!test <48807>
363 %! B = B0; 346 %! B = B0;
364 %! A(3, :) = 2*A(1, :) - A(2, :); 347 %! A(3, :) = 2*A(1, :) - A(2, :);
365 %! [U, V, X, C, S, R] = gsvd (A, B); 348 %! [U, V, X, C, S] = gsvd (A, B);
366 %! D1 = zeros (3, 5); D1(1:3, 1:3) = C; 349 %! assert (C'*C + S'*S, eye (5), 5*eps);
367 %! D2 = zeros (5, 5); D2(1:3, 1:3) = S; D2(4:5, 4:5) = eye (2); 350 %! assert (U*C*X', A, 10*eps);
368 %! assert (norm (diag (C).^2 + diag (S).^2 - ones (3, 1)) <= 1e-6); 351 %! assert (V*S*X', B, 40*eps);
369 %! assert (norm ((U'*A*X) - D1*R) <= 1e-6);
370 %! assert (norm ((V'*B*X) - D2*R) <= 1e-6);
371 352
372 ## A and B are both rank deficient 353 ## A and B are both rank deficient
373 %!test <48807> 354 ## FIXME: LAPACK seems to be completely broken for this case
355 %!#test <48807>
374 %! A = A0.'; B = B0.'; 356 %! A = A0.'; B = B0.';
375 %! A(:, 3) = 2*A(:, 1) - A(:, 2); 357 %! A(:, 3) = 2*A(:, 1) - A(:, 2);
376 %! B(:, 3) = 2*B(:, 1) - B(:, 2); 358 %! B(:, 3) = 2*B(:, 1) - B(:, 2);
377 %! [U, V, X, C, S, R]=gsvd (A, B); 359 %! [U, V, X, C, S] = gsvd (A, B);
378 %! D1 = zeros (3, 4); D1(1:3, 1:3) = C; 360 %! assert (C'*C + S'*S, eye (3), 5*eps);
379 %! D2 = eye (4); D2(1:3, 1:3) = S; D2(5,:) = 0; 361 %! assert (U*C*X', A, 10*eps);
380 %! assert (norm (diag (C).^2 + diag (S).^2 - ones (3, 1)) <= 1e-6); 362 %! assert (V*S*X', B, 20*eps);
381 %! assert (norm ((U'*A*X) - D1*[zeros(4, 1) R]) <= 1e-6);
382 %! assert (norm ((V'*B*X) - D2*[zeros(4, 1) R]) <= 1e-6);
383 363
384 ## A: 5x3 complex full rank, B: 3x3 complex full rank 364 ## A: 5x3 complex full rank, B: 3x3 complex full rank
385 %!test <48807> 365 %!test <48807>
386 %! A0 = A0 + j*randn (5, 3); 366 %! A0 = A0 + j*randn (5, 3);
387 %! B0 = diag ([1 2 4]) + j*diag ([4 -2 -1]); 367 %! B0 = diag ([1 2 4]) + j*diag ([4 -2 -1]);
388 %! A = A0; 368 %! A = A0;
389 %! B = B0; 369 %! B = B0;
390 %! [U, V, X, C, S, R] = gsvd (A, B); 370 %! [U, V, X, C, S] = gsvd (A, B);
391 %! D1 = zeros (5, 3); D1(1:3, 1:3) = C; 371 %! assert (C'*C + S'*S, eye (3), 5*eps);
392 %! D2 = S; 372 %! assert (U*C*X', A, 10*eps);
393 %! assert (norm (diag (C).^2 + diag (S).^2 - ones (3, 1)) <= 1e-6); 373 %! assert (V*S*X', B, 25*eps);
394 %! assert (norm ((U'*A*X) - D1*R) <= 1e-6);
395 %! assert (norm ((V'*B*X) - D2*R) <= 1e-6);
396 374
397 ## A: 5x3 complex full rank, B: 3x3 complex rank deficient 375 ## A: 5x3 complex full rank, B: 3x3 complex rank deficient
398 %!test <48807> 376 %!test <48807>
399 %! B(2, 2) = 0; 377 %! B(2, 2) = 0;
400 %! [U, V, X, C, S, R] = gsvd (A, B); 378 %! [U, V, X, C, S] = gsvd (A, B);
401 %! D1 = zeros (5, 3); D1(1, 1) = 1; D1(2:3, 2:3) = C; 379 %! assert (C'*C + S'*S, eye (3), 5*eps);
402 %! D2 = [zeros(2, 1) S; zeros(1, 3)]; 380 %! assert (U*C*X', A, 10*eps);
403 %! assert (norm (diag (C).^2 + diag (S).^2 - ones (2, 1)) <= 1e-6); 381 %! assert (V*S*X', B, 25*eps);
404 %! assert (norm ((U'*A*X) - D1*R) <= 1e-6);
405 %! assert (norm ((V'*B*X) - D2*R) <= 1e-6);
406 382
407 ## A: 5x3 complex rank deficient, B: 3x3 complex full rank 383 ## A: 5x3 complex rank deficient, B: 3x3 complex full rank
408 %!test <48807> 384 %!test <48807>
409 %! B = B0; 385 %! B = B0;
410 %! A(:, 3) = 2*A(:, 1) - A(:, 2); 386 %! A(:, 3) = 2*A(:, 1) - A(:, 2);
411 %! [U, V, X, C, S, R] = gsvd (A, B); 387 %! [U, V, X, C, S] = gsvd (A, B);
412 %! D1 = zeros (5, 3); D1(1:3, 1:3) = C; 388 %! assert (C'*C + S'*S, eye (3), 5*eps);
413 %! D2 = S; 389 %! assert (U*C*X', A, 10*eps);
414 %! assert (norm (diag (C).^2 + diag (S).^2 - ones (3, 1)) <= 1e-6); 390 %! assert (V*S*X', B, 25*eps);
415 %! assert (norm ((U'*A*X) - D1*R) <= 1e-6);
416 %! assert (norm ((V'*B*X) - D2*R) <= 1e-6);
417 391
418 ## A (5x3) and B (3x3) are both complex rank deficient 392 ## A (5x3) and B (3x3) are both complex rank deficient
419 %!test <48807> 393 ## FIXME: LAPACK seems to be completely broken for this case
394 %!#test <48807>
420 %! B(:, 3) = 2*B(:, 1) - B(:, 2); 395 %! B(:, 3) = 2*B(:, 1) - B(:, 2);
421 %! [U, V, X, C, S, R] = gsvd (A, B); 396 %! [U, V, X, C, S] = gsvd (A, B);
422 %! D1 = zeros (5, 2); D1(1:2, 1:2) = C; 397 %! assert (C'*C + S'*S, eye (3), 5*eps);
423 %! D2 = [S; zeros(1, 2)]; 398 %! assert (U*C*X', A, 10*eps);
424 %! assert (norm (diag (C).^2 + diag (S).^2 - ones (2, 1)) <= 1e-6); 399 %! assert (V*S*X', B, 20*eps);
425 %! assert (norm ((U'*A*X) - D1*[zeros(2, 1) R]) <= 1e-6);
426 %! assert (norm ((V'*B*X) - D2*[zeros(2, 1) R]) <= 1e-6);
427 400
428 ## A (now 3x5) complex and B (now 5x5) complex are full rank 401 ## A (now 3x5) complex and B (now 5x5) complex are full rank
429 ## now, A is 3x5 402 ## now, A is 3x5
430 %!test <48807> 403 %!test <48807>
431 %! A = A0.'; 404 %! A = A0.';
432 %! B0 = diag ([1 2 4 8 16]) + j*diag ([-5 4 -3 2 -1]); 405 %! B0 = diag ([1 2 4 8 16]) + j*diag ([-5 4 -3 2 -1]);
433 %! B = B0; 406 %! B = B0;
434 %! [U, V, X, C, S, R] = gsvd (A, B); 407 %! [U, V, X, C, S] = gsvd (A, B);
435 %! D1 = [C zeros(3,2)]; 408 %! assert (C'*C + S'*S, eye (5), 5*eps);
436 %! D2 = [S zeros(3,2); zeros(2, 3) eye(2)]; 409 %! assert (U*C*X', A, 25*eps);
437 %! assert (norm (diag (C).^2 + diag (S).^2 - ones (3, 1)) <= 1e-6); 410 %! assert (V*S*X', B, 85*eps);
438 %! assert (norm ((U'*A*X) - D1*R) <= 1e-6);
439 %! assert (norm ((V'*B*X) - D2*R) <= 1e-6);
440 411
441 ## A: 3x5 complex full rank, B: 5x5 complex rank deficient 412 ## A: 3x5 complex full rank, B: 5x5 complex rank deficient
442 %!test <48807> 413 %!test <48807>
443 %! B(2, 2) = 0; 414 %! B(2, 2) = 0;
444 %! [U, V, X, C, S, R] = gsvd (A, B); 415 %! [U, V, X, C, S] = gsvd (A, B);
445 %! D1 = zeros (3, 5); D1(1, 1) = 1; D1(2:3, 2:3) = C; 416 %! assert (C'*C + S'*S, eye (5), 5*eps);
446 %! D2 = zeros (5,5); D2(1:2, 2:3) = S; D2(3:4, 4:5) = eye (2); 417 %! assert (U*C*X', A, 10*eps);
447 %! assert (norm (diag (C).^2 + diag (S).^2 - ones (2, 1)) <= 1e-6); 418 %! assert (V*S*X', B, 85*eps);
448 %! assert (norm ((U'*A*X) - D1*R) <= 1e-6);
449 %! assert (norm ((V'*B*X) - D2*R) <= 1e-6);
450 419
451 ## A: 3x5 complex rank deficient, B: 5x5 complex full rank 420 ## A: 3x5 complex rank deficient, B: 5x5 complex full rank
452 %!test <48807> 421 %!test <48807>
453 %! B = B0; 422 %! B = B0;
454 %! A(3, :) = 2*A(1, :) - A(2, :); 423 %! A(3, :) = 2*A(1, :) - A(2, :);
455 %! [U, V, X, C, S, R] = gsvd (A, B); 424 %! [U, V, X, C, S] = gsvd (A, B);
456 %! D1 = zeros (3, 5); D1(1:3, 1:3) = C; 425 %! assert (C'*C + S'*S, eye (5), 5*eps);
457 %! D2 = zeros (5,5); D2(1:3, 1:3) = S; D2(4:5, 4:5) = eye (2); 426 %! assert (U*C*X', A, 10*eps);
458 %! assert (norm (diag (C).^2 + diag (S).^2 - ones (3, 1)) <= 1e-6); 427 %! assert (V*S*X', B, 85*eps);
459 %! assert (norm ((U'*A*X) - D1*R) <= 1e-6);
460 %! assert (norm ((V'*B*X) - D2*R) <= 1e-6);
461 428
462 ## A and B are both complex rank deficient 429 ## A and B are both complex rank deficient
463 %!test <48807> 430 ## FIXME: LAPACK seems to be completely broken for this case
431 %!#test <48807>
464 %! A = A0.'; 432 %! A = A0.';
465 %! B = B0.'; 433 %! B = B0.';
466 %! A(:, 3) = 2*A(:, 1) - A(:, 2); 434 %! A(:, 3) = 2*A(:, 1) - A(:, 2);
467 %! B(:, 3) = 2*B(:, 1) - B(:, 2); 435 %! B(:, 3) = 2*B(:, 1) - B(:, 2);
468 %! [U, V, X, C, S, R] = gsvd (A, B); 436 %! [U, V, X, C, S] = gsvd (A, B);
469 %! D1 = zeros (3, 4); D1(1:3, 1:3) = C; 437 %! assert (C'*C + S'*S, eye (5), 5*eps);
470 %! D2 = eye (4); D2(1:3, 1:3) = S; D2(5,:) = 0; 438 %! assert (U*C*X', A, 10*eps);
471 %! assert (norm (diag (C).^2 + diag (S).^2 - ones (3, 1)) <= 1e-6); 439 %! assert (V*S*X', B, 85*eps);
472 %! assert (norm ((U'*A*X) - D1*[zeros(4, 1) R]) <= 1e-6);
473 %! assert (norm ((V'*B*X) - D2*[zeros(4, 1) R]) <= 1e-6);
474 440
475 ## Test that single inputs produce single outputs 441 ## Test that single inputs produce single outputs
476 %!test 442 %!test
477 %! s = gsvd (single (ones (0,1)), B); 443 %! s = gsvd (single (eye (5)), B);
478 %! assert (class (s), "single"); 444 %! assert (class (s), "single");
479 %! s = gsvd (single (ones (1,0)), B); 445 %! [U,V,X,C,S] = gsvd (single (eye(5)), B);
480 %! assert (class (s), "single");
481 %! s = gsvd (single (ones (1,0)), B);
482 %! [U,V,X,C,S,R] = gsvd (single ([]), B);
483 %! assert (class (U), "single"); 446 %! assert (class (U), "single");
484 %! assert (class (V), "single"); 447 %! assert (class (V), "single");
485 %! assert (class (X), "single"); 448 %! assert (class (X), "single");
486 %! assert (class (C), "single"); 449 %! assert (class (C), "single");
487 %! assert (class (S), "single"); 450 %! assert (class (S), "single");
488 %! assert (class (R), "single");
489 %! 451 %!
490 %! s = gsvd (single (A), B); 452 %! s = gsvd (A, single (eye (5)));
491 %! assert (class (s), "single"); 453 %! assert (class (s), "single");
492 %! [U,V,X,C,S,R] = gsvd (single (A), B); 454 %! [U,V,X,C,S] = gsvd (A, single (eye (5)));
493 %! assert (class (U), "single"); 455 %! assert (class (U), "single");
494 %! assert (class (V), "single"); 456 %! assert (class (V), "single");
495 %! assert (class (X), "single"); 457 %! assert (class (X), "single");
496 %! assert (class (C), "single"); 458 %! assert (class (C), "single");
497 %! assert (class (S), "single"); 459 %! assert (class (S), "single");
498 %! assert (class (R), "single"); 460
461 ## Test input validation
462 %!error <Invalid call> gsvd ()
463 %!error <Invalid call> gsvd (1)
464 %!error <Invalid call> gsvd (1,2,3,4)
465 %!warning <economy-sized decomposition is not yet implemented> gsvd (1,2,0);
466 %!error <A and B must have the same number of columns> gsvd (1,[1, 2])
467 ## Test input validation for single (real and complex) inputs.
468 %!error <A cannot have Inf or NaN values> gsvd (Inf, single (2))
469 %!error <A cannot have Inf or NaN values> gsvd (NaN, single (2))
470 %!error <B cannot have Inf or NaN values> gsvd (single (1), Inf)
471 %!error <B cannot have Inf or NaN values> gsvd (single (1), NaN)
472 %!error <A must be a real or complex matrix> gsvd ({1}, single (2i))
473 %!error <B must be a real or complex matrix> gsvd (single (i), {2})
474 %!error <A cannot have Inf or NaN values> gsvd (Inf, single (2i))
475 %!error <A cannot have Inf or NaN values> gsvd (NaN, single (2i))
476 %!error <B cannot have Inf or NaN values> gsvd (single (i), Inf)
477 %!error <B cannot have Inf or NaN values> gsvd (single (i), NaN)
478 ## Test input validation for single, but not real or complex, inputs.
479 %!error <A and B must be real or complex matrices> gsvd ({1}, single (2))
480 %!error <A and B must be real or complex matrices> gsvd (single (1), {2})
481 ## Test input validation for double (real and complex) inputs.
482 %!error <A cannot have Inf or NaN values> gsvd (Inf, 2)
483 %!error <A cannot have Inf or NaN values> gsvd (NaN, 2)
484 %!error <B cannot have Inf or NaN values> gsvd (1, Inf)
485 %!error <B cannot have Inf or NaN values> gsvd (1, NaN)
486 %!error <A must be a real or complex matrix> gsvd ({1}, 2i)
487 %!error <B must be a real or complex matrix> gsvd (i, {2})
488 %!error <A cannot have Inf or NaN values> gsvd (Inf, 2i)
489 %!error <A cannot have Inf or NaN values> gsvd (NaN, 2i)
490 %!error <B cannot have Inf or NaN values> gsvd (i, Inf)
491 %!error <B cannot have Inf or NaN values> gsvd (i, NaN)
492 ## Test input validation for double, but not real or complex, inputs.
493 %!error <A and B must be real or complex matrices> gsvd ({1}, double (2))
494 %!error <A and B must be real or complex matrices> gsvd (double (1), {2})
495 ## Test input validation in liboctave/numeric/gsvd.cc
496 %!error <A and B cannot be empty matrices> gsvd (zeros (0,1), 1)
497 %!error <A and B cannot be empty matrices> gsvd (1, zeros (0,1))
499 498
500 */ 499 */
501 500
502 OCTAVE_NAMESPACE_END 501 OCTAVE_NAMESPACE_END