Mercurial > octave

comparison libinterp/dldfcn/qr.cc @ 22206:21684fa513ce

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Return C = Q'*B not Q when qr has two arguments (bug #41567, bug #46912) * qr.cc (Fqr): For full case, explicitly calculate Q. Add 8 tests.
author mfasi <mogrob.sanit@gmail.com>
date Thu, 18 Feb 2016 11:57:12 +1100
parents 278fc29b69ca
children 99454a60bf5e
comparison
equal deleted inserted replaced
22205:8a456af24e6b 22206:21684fa513ce
212 int arg_is_empty = empty_arg ("qr", arg.rows (), arg.columns ()); 212 int arg_is_empty = empty_arg ("qr", arg.rows (), arg.columns ());
213 213
214 if (arg_is_empty < 0) 214 if (arg_is_empty < 0)
215 return retval; 215 return retval;
216 216
217 bool economy = false;
218 bool is_cmplx = false;
219 bool b_mat = false;
220 int have_b = 0;
221
222 if (arg.is_complex_type ())
223 is_cmplx = true;
224 if (nargin > 1)
225 {
226 have_b = 1;
227 if (args(nargin-1).is_scalar_type ())
228 {
229 int val = args(nargin-1).int_value ();
230 if (val == 0)
231 {
232 economy = true;
233 have_b = (nargin > 2 ? 2 : 0);
234 }
235 }
236 else if (args(nargin-1).is_matrix_type ())
237 b_mat = true;
238 if (have_b > 0 && args(have_b).is_complex_type ())
239 is_cmplx = true;
240 }
241
217 if (arg.is_sparse_type ()) 242 if (arg.is_sparse_type ())
218 { 243 {
219 bool economy = false;
220 bool is_cmplx = false;
221 int have_b = 0;
222
223 if (arg.is_complex_type ())
224 is_cmplx = true;
225 if (nargin > 1)
226 {
227 have_b = 1;
228 if (args(nargin-1).is_scalar_type ())
229 {
230 int val = args(nargin-1).int_value ();
231 if (val == 0)
232 {
233 economy = true;
234 have_b = (nargin > 2 ? 2 : 0);
235 }
236 }
237 if (have_b > 0 && args(have_b).is_complex_type ())
238 is_cmplx = true;
239 }
240 244
241 if (is_cmplx) 245 if (is_cmplx)
242 { 246 {
243 sparse_qr<SparseComplexMatrix> q (arg.sparse_complex_matrix_value ()); 247 sparse_qr<SparseComplexMatrix> q (arg.sparse_complex_matrix_value ());
244 248
293 297
294 case 2: 298 case 2:
295 { 299 {
296 qr<FloatMatrix> fact (m, type); 300 qr<FloatMatrix> fact (m, type);
297 retval = ovl (fact.Q (), get_qr_r (fact)); 301 retval = ovl (fact.Q (), get_qr_r (fact));
302 if (b_mat)
303 {
304 if (is_cmplx)
305 retval(0) = fact.Q ().transpose ()
306 * args(1).float_complex_matrix_value ();
307 else
308 retval(0) = fact.Q ().transpose ()
309 * args(1).float_matrix_value ();
310 }
298 } 311 }
299 break; 312 break;
300 313
301 default: 314 default:
302 { 315 {
329 342
330 case 2: 343 case 2:
331 { 344 {
332 qr<FloatComplexMatrix> fact (m, type); 345 qr<FloatComplexMatrix> fact (m, type);
333 retval = ovl (fact.Q (), get_qr_r (fact)); 346 retval = ovl (fact.Q (), get_qr_r (fact));
347 if (b_mat)
348 retval (0) = conj (fact.Q ().transpose ())
349 * args(1).float_complex_matrix_value ();
334 } 350 }
335 break; 351 break;
336 352
337 default: 353 default:
338 { 354 {
366 382
367 case 2: 383 case 2:
368 { 384 {
369 qr<Matrix> fact (m, type); 385 qr<Matrix> fact (m, type);
370 retval = ovl (fact.Q (), get_qr_r (fact)); 386 retval = ovl (fact.Q (), get_qr_r (fact));
387 if (b_mat)
388 {
389 if (is_cmplx)
390 retval(0) = fact.Q ().transpose ()
391 * args(1).complex_matrix_value ();
392 else
393 retval(0) = fact.Q ().transpose ()
394 * args(1).matrix_value ();
395 }
371 } 396 }
372 break; 397 break;
373 398
374 default: 399 default:
375 { 400 {
401 426
402 case 2: 427 case 2:
403 { 428 {
404 qr<ComplexMatrix> fact (m, type); 429 qr<ComplexMatrix> fact (m, type);
405 retval = ovl (fact.Q (), get_qr_r (fact)); 430 retval = ovl (fact.Q (), get_qr_r (fact));
431 if (b_mat)
432 retval (0) = conj (fact.Q ().transpose ())
433 * args(1).complex_matrix_value ();
406 } 434 }
407 break; 435 break;
408 436
409 default: 437 default:
410 { 438 {
436 %! assert (qe * re, a, sqrt (eps)); 464 %! assert (qe * re, a, sqrt (eps));
437 465
438 %!test 466 %!test
439 %! a = [0, 2, 1; 2, 1, 2]; 467 %! a = [0, 2, 1; 2, 1, 2];
440 %! 468 %!
469 %! [q, r] = qr (a);
470 %! [qe, re] = qr (a, 0);
471 %!
472 %! assert (q * r, a, sqrt (eps));
473 %! assert (qe * re, a, sqrt (eps));
474
475 %!test
476 %! a = [0, 2, 1; 2, 1, 2];
477 %!
441 %! [q, r, p] = qr (a); # FIXME: not giving right dimensions. 478 %! [q, r, p] = qr (a); # FIXME: not giving right dimensions.
442 %! [qe, re, pe] = qr (a, 0); 479 %! [qe, re, pe] = qr (a, 0);
443 %! 480 %!
444 %! assert (q * r, a * p, sqrt (eps)); 481 %! assert (q * r, a * p, sqrt (eps));
445 %! assert (qe * re, a(:, pe), sqrt (eps)); 482 %! assert (qe * re, a(:, pe), sqrt (eps));
459 %! [q, r, p] = qr (a); 496 %! [q, r, p] = qr (a);
460 %! [qe, re, pe] = qr (a, 0); 497 %! [qe, re, pe] = qr (a, 0);
461 %! 498 %!
462 %! assert (q * r, a * p, sqrt (eps)); 499 %! assert (q * r, a * p, sqrt (eps));
463 %! assert (qe * re, a(:, pe), sqrt (eps)); 500 %! assert (qe * re, a(:, pe), sqrt (eps));
501
502 %!test
503 %! a = [0, 2, 1; 2, 1, 2; 3, 1, 2];
504 %! b = [1, 3, 2; 1, 1, 0; 3, 0, 2];
505 %!
506 %! [q, r] = qr (a);
507 %! [c, re] = qr (a, b);
508 %!
509 %! assert (r, re, sqrt (eps));
510 %! assert (q'*b, c, sqrt (eps));
511
512 %!test
513 %! a = [0, 2, i; 2, 1, 2; 3, 1, 2];
514 %! b = [1, 3, 2; 1, i, 0; 3, 0, 2];
515 %!
516 %! [q, r] = qr (a);
517 %! [c, re] = qr (a, b);
518 %!
519 %! assert (r, re, sqrt (eps));
520 %! assert (q'*b, c, sqrt (eps));
521
522 %!test
523 %! a = [0, 2, i; 2, 1, 2; 3, 1, 2];
524 %! b = [1, 3, 2; 1, 1, 0; 3, 0, 2];
525 %!
526 %! [q, r] = qr (a);
527 %! [c, re] = qr (a, b);
528 %!
529 %! assert (r, re, sqrt (eps));
530 %! assert (q'*b, c, sqrt (eps));
531
532 %!test
533 %! a = [0, 2, 1; 2, 1, 2; 3, 1, 2];
534 %! b = [1, 3, 2; 1, i, 0; 3, 0, 2];
535 %!
536 %! [q, r] = qr (a);
537 %! [c, re] = qr (a, b);
538 %!
539 %! assert (r, re, sqrt (eps));
540 %! assert (q'*b, c, sqrt (eps));
464 541
465 %!error qr () 542 %!error qr ()
466 %!error qr ([1, 2; 3, 4], 0, 2) 543 %!error qr ([1, 2; 3, 4], 0, 2)
467 544
468 %!function retval = __testqr (q, r, a, p) 545 %!function retval = __testqr (q, r, a, p)
573 %! [qe, re, pe] = qr (a, 0); 650 %! [qe, re, pe] = qr (a, 0);
574 %! 651 %!
575 %! assert (q * r, a * p, sqrt (eps ("single"))); 652 %! assert (q * r, a * p, sqrt (eps ("single")));
576 %! assert (qe * re, a(:, pe), sqrt (eps ("single"))); 653 %! assert (qe * re, a(:, pe), sqrt (eps ("single")));
577 654
655 %!test
656 %! a = single([0, 2, 1; 2, 1, 2; 3, 1, 2]);
657 %! b = single([1, 3, 2; 1, 1, 0; 3, 0, 2]);
658 %!
659 %! [q, r] = qr (a);
660 %! [c, re] = qr (a, b);
661 %!
662 %! assert (r, re, sqrt (eps ("single")));
663 %! assert (q'*b, c, sqrt (eps ("single")));
664
665 %!test
666 %! a = single([0, 2, i; 2, 1, 2; 3, 1, 2]);
667 %! b = single([1, 3, 2; 1, i, 0; 3, 0, 2]);
668 %!
669 %! [q, r] = qr (a);
670 %! [c, re] = qr (a, b);
671 %!
672 %! assert (r, re, sqrt (eps ("single")));
673 %! assert (q'*b, c, sqrt (eps ("single")));
674
675 %!test
676 %! a = single([0, 2, i; 2, 1, 2; 3, 1, 2]);
677 %! b = single([1, 3, 2; 1, 1, 0; 3, 0, 2]);
678 %!
679 %! [q, r] = qr (a);
680 %! [c, re] = qr (a, b);
681 %!
682 %! assert (r, re, sqrt (eps));
683 %! assert (q'*b, c, sqrt (eps));
684
685 %!test
686 %! a = single([0, 2, 1; 2, 1, 2; 3, 1, 2]);
687 %! b = single([1, 3, 2; 1, i, 0; 3, 0, 2]);
688 %!
689 %! [q, r] = qr (a);
690 %! [c, re] = qr (a, b);
691 %!
692 %! assert (r, re, sqrt (eps ("single")));
693 %! assert (q'*b, c, sqrt (eps ("single")));
694
578 %!error qr () 695 %!error qr ()
579 %!error qr ([1, 2; 3, 4], 0, 2) 696 %!error qr ([1, 2; 3, 4], 0, 2)
580 697
581 %!test 698 %!test
582 %! t = ones (24, 1); 699 %! t = ones (24, 1);