Mercurial > hg > octave-nkf

comparison liboctave/dbleQR.cc @ 8562:a6edd5c23cb5

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
use replacement methods if qrupdate is not available
author Jaroslav Hajek <highegg@gmail.com>
date Thu, 22 Jan 2009 11:10:47 +0100
parents d66c9b6e506a
children c86718093c1b
comparison
equal deleted inserted replaced
8561:66165de2cc42 8562:a6edd5c23cb5
157 } 157 }
158 } 158 }
159 159
160 QR::QR (const Matrix& q_arg, const Matrix& r_arg) 160 QR::QR (const Matrix& q_arg, const Matrix& r_arg)
161 { 161 {
162 if (q_arg.columns () != r_arg.rows ()) 162 octave_idx_type qr = q_arg.rows (), qc = q_arg.columns ();
163 { 163 octave_idx_type rr = r_arg.rows (), rc = r_arg.columns ();
164 (*current_liboctave_error_handler) ("QR dimensions mismatch"); 164 if (qc == rr && (qr == qc || (qr > qc && rr == rc)))
165 return; 165 {
166 } 166 q = q_arg;
167 167 r = r_arg;
168 this->q = q_arg; 168 }
169 this->r = r_arg; 169 else
170 (*current_liboctave_error_handler) ("QR dimensions mismatch");
171 }
172
173 QR::type
174 QR::get_type (void) const
175 {
176 QR::type retval;
177 if (!q.is_empty () && q.is_square ())
178 retval = QR::std;
179 else if (q.rows () > q.columns () && r.is_square ())
180 retval = QR::economy;
181 else
182 retval = QR::raw;
183 return retval;
170 } 184 }
171 185
172 #ifdef HAVE_QRUPDATE 186 #ifdef HAVE_QRUPDATE
173 187
174 void 188 void
184 OCTAVE_LOCAL_BUFFER (double, w, 2*k); 198 OCTAVE_LOCAL_BUFFER (double, w, 2*k);
185 F77_XFCN (dqr1up, DQR1UP, (m, n, k, q.fortran_vec (), m, r.fortran_vec (), k, 199 F77_XFCN (dqr1up, DQR1UP, (m, n, k, q.fortran_vec (), m, r.fortran_vec (), k,
186 utmp.fortran_vec (), vtmp.fortran_vec (), w)); 200 utmp.fortran_vec (), vtmp.fortran_vec (), w));
187 } 201 }
188 else 202 else
189 (*current_liboctave_error_handler) ("QR update dimensions mismatch"); 203 (*current_liboctave_error_handler) ("qrupdate: dimensions mismatch");
190 } 204 }
191 205
192 void 206 void
193 QR::update (const Matrix& u, const Matrix& v) 207 QR::update (const Matrix& u, const Matrix& v)
194 { 208 {
416 r.fortran_vec (), r.rows (), 430 r.fortran_vec (), r.rows (),
417 i + 1, j + 1, w)); 431 i + 1, j + 1, w));
418 } 432 }
419 } 433 }
420 434
435 #else
436
437 // Replacement update methods.
438
439 void
440 QR::update (const ColumnVector& u, const ColumnVector& v)
441 {
442 warn_qrupdate_once ();
443
444 octave_idx_type m = q.rows ();
445 octave_idx_type n = r.columns ();
446
447 if (u.length () == m && v.length () == n)
448 {
449 init(q*r + Matrix (u) * Matrix (v).transpose (), get_type ());
450 }
451 else
452 (*current_liboctave_error_handler) ("qrupdate: dimensions mismatch");
453 }
454
455 void
456 QR::update (const Matrix& u, const Matrix& v)
457 {
458 warn_qrupdate_once ();
459
460 octave_idx_type m = q.rows ();
461 octave_idx_type n = r.columns ();
462
463 if (u.rows () == m && v.rows () == n && u.cols () == v.cols ())
464 {
465 init(q*r + u * v.transpose (), get_type ());
466 }
467 else
468 (*current_liboctave_error_handler) ("qrupdate: dimensions mismatch");
469 }
470
471 static
472 Matrix insert_col (const Matrix& a, octave_idx_type i,
473 const ColumnVector& x)
474 {
475 Matrix retval (a.rows (), a.columns () + 1);
476 retval.assign (idx_vector::colon, idx_vector (0, i),
477 a.index (idx_vector::colon, idx_vector (0, i)));
478 retval.assign (idx_vector::colon, idx_vector (i), x);
479 retval.assign (idx_vector::colon, idx_vector (i+1, retval.columns ()),
480 a.index (idx_vector::colon, idx_vector (i, a.columns ())));
481 return retval;
482 }
483
484 static
485 Matrix insert_row (const Matrix& a, octave_idx_type i,
486 const RowVector& x)
487 {
488 Matrix retval (a.rows () + 1, a.columns ());
489 retval.assign (idx_vector (0, i), idx_vector::colon,
490 a.index (idx_vector (0, i), idx_vector::colon));
491 retval.assign (idx_vector (i), idx_vector::colon, x);
492 retval.assign (idx_vector (i+1, retval.rows ()), idx_vector::colon,
493 a.index (idx_vector (i, a.rows ()), idx_vector::colon));
494 return retval;
495 }
496
497 static
498 Matrix delete_col (const Matrix& a, octave_idx_type i)
499 {
500 Matrix retval = a;
501 retval.delete_elements (1, idx_vector (i));
502 return retval;
503 }
504
505 static
506 Matrix delete_row (const Matrix& a, octave_idx_type i)
507 {
508 Matrix retval = a;
509 retval.delete_elements (0, idx_vector (i));
510 return retval;
511 }
512
513 static
514 Matrix shift_cols (const Matrix& a,
515 octave_idx_type i, octave_idx_type j)
516 {
517 octave_idx_type n = a.columns ();
518 Array<octave_idx_type> p (n);
519 for (octave_idx_type k = 0; k < n; k++) p(k) = k;
520 if (i < j)
521 {
522 for (octave_idx_type k = i; k < j; k++) p(k) = k+1;
523 p(j) = i;
524 }
525 else if (j < i)
526 {
527 p(j) = i;
528 for (octave_idx_type k = j+1; k < i+1; k++) p(k) = k-1;
529 }
530
531 return a.index (idx_vector::colon, idx_vector (p));
532 }
533
534 void
535 QR::insert_col (const ColumnVector& u, octave_idx_type j)
536 {
537 warn_qrupdate_once ();
538
539 octave_idx_type m = q.rows ();
540 octave_idx_type n = r.columns ();
541
542 if (u.length () != m)
543 (*current_liboctave_error_handler) ("qrinsert: dimensions mismatch");
544 else if (j < 0 || j > n)
545 (*current_liboctave_error_handler) ("qrinsert: index out of range");
546 else
547 {
548 init (::insert_col (q*r, j, u), get_type ());
549 }
550 }
551
552 void
553 QR::insert_col (const Matrix& u, const Array<octave_idx_type>& j)
554 {
555 warn_qrupdate_once ();
556
557 octave_idx_type m = q.rows ();
558 octave_idx_type n = r.columns ();
559
560 Array<octave_idx_type> jsi;
561 Array<octave_idx_type> js = j.sort (jsi, ASCENDING);
562 octave_idx_type nj = js.length ();
563 bool dups = false;
564 for (octave_idx_type i = 0; i < nj - 1; i++)
565 dups = dups && js(i) == js(i+1);
566
567 if (dups)
568 (*current_liboctave_error_handler) ("qrinsert: duplicate index detected");
569 else if (u.length () != m || u.columns () != nj)
570 (*current_liboctave_error_handler) ("qrinsert: dimensions mismatch");
571 else if (nj > 0 && (js(0) < 0 || js(nj-1) > n))
572 (*current_liboctave_error_handler) ("qrinsert: index out of range");
573 else if (nj > 0)
574 {
575 Matrix a = q*r;
576 for (octave_idx_type i = 0; i < js.length (); i++)
577 a = ::insert_col (a, js(i), u.column (i));
578 init (a, get_type ());
579 }
580 }
581
582 void
583 QR::delete_col (octave_idx_type j)
584 {
585 warn_qrupdate_once ();
586
587 octave_idx_type m = q.rows ();
588 octave_idx_type n = r.columns ();
589
590 if (j < 0 || j > n-1)
591 (*current_liboctave_error_handler) ("qrdelete: index out of range");
592 else
593 {
594 init (::delete_col (q*r, j), get_type ());
595 }
596 }
597
598 void
599 QR::delete_col (const Array<octave_idx_type>& j)
600 {
601 warn_qrupdate_once ();
602
603 octave_idx_type m = q.rows ();
604 octave_idx_type n = r.columns ();
605
606 Array<octave_idx_type> jsi;
607 Array<octave_idx_type> js = j.sort (jsi, DESCENDING);
608 octave_idx_type nj = js.length ();
609 bool dups = false;
610 for (octave_idx_type i = 0; i < nj - 1; i++)
611 dups = dups && js(i) == js(i+1);
612
613 if (dups)
614 (*current_liboctave_error_handler) ("qrinsert: duplicate index detected");
615 else if (nj > 0 && (js(0) > n-1 || js(nj-1) < 0))
616 (*current_liboctave_error_handler) ("qrinsert: index out of range");
617 else if (nj > 0)
618 {
619 Matrix a = q*r;
620 for (octave_idx_type i = 0; i < js.length (); i++)
621 a = ::delete_col (a, js(i));
622 init (a, get_type ());
623 }
624 }
625
626 void
627 QR::insert_row (const RowVector& u, octave_idx_type j)
628 {
629 warn_qrupdate_once ();
630
631 octave_idx_type m = r.rows ();
632 octave_idx_type n = r.columns ();
633
634 if (! q.is_square () || u.length () != n)
635 (*current_liboctave_error_handler) ("qrinsert: dimensions mismatch");
636 else if (j < 0 || j > m)
637 (*current_liboctave_error_handler) ("qrinsert: index out of range");
638 else
639 {
640 init (::insert_row (q*r, j, u), get_type ());
641 }
642 }
643
644 void
645 QR::delete_row (octave_idx_type j)
646 {
647 octave_idx_type m = r.rows ();
648 octave_idx_type n = r.columns ();
649
650 if (! q.is_square ())
651 (*current_liboctave_error_handler) ("qrdelete: dimensions mismatch");
652 else if (j < 0 || j > m-1)
653 (*current_liboctave_error_handler) ("qrdelete: index out of range");
654 else
655 {
656 init (::delete_row (q*r, j), get_type ());
657 }
658 }
659
660 void
661 QR::shift_cols (octave_idx_type i, octave_idx_type j)
662 {
663 warn_qrupdate_once ();
664
665 octave_idx_type m = q.rows ();
666 octave_idx_type n = r.columns ();
667
668 if (i < 0 || i > n-1 || j < 0 || j > n-1)
669 (*current_liboctave_error_handler) ("qrshift: index out of range");
670 else
671 {
672 init (::shift_cols (q*r, i, j), get_type ());
673 }
674 }
675
676 void warn_qrupdate_once (void)
677 {
678 static bool warned = false;
679 if (! warned)
680 {
681 (*current_liboctave_warning_handler)
682 ("In this version of Octave, QR & Cholesky updating routines\n"
683 "simply update the matrix and recalculate factorizations.\n"
684 "To use fast algorithms, link Octave with the qrupdate library.\n"
685 "See <http://sourceforge.net/projects/qrupdate>.\n");
686 warned = true;
687 }
688 }
689
421 #endif 690 #endif
422 691
423 /* 692 /*
424 ;;; Local Variables: *** 693 ;;; Local Variables: ***
425 ;;; mode: C++ *** 694 ;;; mode: C++ ***