Mercurial > hg > octave-lyh

comparison liboctave/fCmplxQR.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
166 } 166 }
167 } 167 }
168 168
169 FloatComplexQR::FloatComplexQR (const FloatComplexMatrix& q_arg, const FloatComplexMatrix& r_arg) 169 FloatComplexQR::FloatComplexQR (const FloatComplexMatrix& q_arg, const FloatComplexMatrix& r_arg)
170 { 170 {
171 if (q_arg.columns () != r_arg.rows ()) 171 octave_idx_type qr = q_arg.rows (), qc = q_arg.columns ();
172 { 172 octave_idx_type rr = r_arg.rows (), rc = r_arg.columns ();
173 (*current_liboctave_error_handler) ("QR dimensions mismatch"); 173 if (qc == rr && (qr == qc || (qr > qc && rr == rc)))
174 return; 174 {
175 } 175 q = q_arg;
176 176 r = r_arg;
177 this->q = q_arg; 177 }
178 this->r = r_arg; 178 else
179 (*current_liboctave_error_handler) ("QR dimensions mismatch");
180 }
181
182 QR::type
183 FloatComplexQR::get_type (void) const
184 {
185 QR::type retval;
186 if (!q.is_empty () && q.is_square ())
187 retval = QR::std;
188 else if (q.rows () > q.columns () && r.is_square ())
189 retval = QR::economy;
190 else
191 retval = QR::raw;
192 return retval;
179 } 193 }
180 194
181 #ifdef HAVE_QRUPDATE 195 #ifdef HAVE_QRUPDATE
182 196
183 void 197 void
194 OCTAVE_LOCAL_BUFFER (float, rw, k); 208 OCTAVE_LOCAL_BUFFER (float, rw, k);
195 F77_XFCN (cqr1up, CQR1UP, (m, n, k, q.fortran_vec (), m, r.fortran_vec (), k, 209 F77_XFCN (cqr1up, CQR1UP, (m, n, k, q.fortran_vec (), m, r.fortran_vec (), k,
196 utmp.fortran_vec (), vtmp.fortran_vec (), w, rw)); 210 utmp.fortran_vec (), vtmp.fortran_vec (), w, rw));
197 } 211 }
198 else 212 else
199 (*current_liboctave_error_handler) ("QR update dimensions mismatch"); 213 (*current_liboctave_error_handler) ("qrupdate: dimensions mismatch");
200 } 214 }
201 215
202 void 216 void
203 FloatComplexQR::update (const FloatComplexMatrix& u, const FloatComplexMatrix& v) 217 FloatComplexQR::update (const FloatComplexMatrix& u, const FloatComplexMatrix& v)
204 { 218 {
428 r.fortran_vec (), r.rows (), 442 r.fortran_vec (), r.rows (),
429 i + 1, j + 1, w, rw)); 443 i + 1, j + 1, w, rw));
430 } 444 }
431 } 445 }
432 446
447 #else
448
449 // Replacement update methods.
450
451 void
452 FloatComplexQR::update (const FloatComplexColumnVector& u, const FloatComplexColumnVector& v)
453 {
454 warn_qrupdate_once ();
455
456 octave_idx_type m = q.rows ();
457 octave_idx_type n = r.columns ();
458
459 if (u.length () == m && v.length () == n)
460 {
461 init(q*r + FloatComplexMatrix (u) * FloatComplexMatrix (v).hermitian (), get_type ());
462 }
463 else
464 (*current_liboctave_error_handler) ("qrupdate: dimensions mismatch");
465 }
466
467 void
468 FloatComplexQR::update (const FloatComplexMatrix& u, const FloatComplexMatrix& v)
469 {
470 warn_qrupdate_once ();
471
472 octave_idx_type m = q.rows ();
473 octave_idx_type n = r.columns ();
474
475 if (u.rows () == m && v.rows () == n && u.cols () == v.cols ())
476 {
477 init(q*r + u * v.hermitian (), get_type ());
478 }
479 else
480 (*current_liboctave_error_handler) ("qrupdate: dimensions mismatch");
481 }
482
483 static
484 FloatComplexMatrix insert_col (const FloatComplexMatrix& a, octave_idx_type i,
485 const FloatComplexColumnVector& x)
486 {
487 FloatComplexMatrix retval (a.rows (), a.columns () + 1);
488 retval.assign (idx_vector::colon, idx_vector (0, i),
489 a.index (idx_vector::colon, idx_vector (0, i)));
490 retval.assign (idx_vector::colon, idx_vector (i), x);
491 retval.assign (idx_vector::colon, idx_vector (i+1, retval.columns ()),
492 a.index (idx_vector::colon, idx_vector (i, a.columns ())));
493 return retval;
494 }
495
496 static
497 FloatComplexMatrix insert_row (const FloatComplexMatrix& a, octave_idx_type i,
498 const FloatComplexRowVector& x)
499 {
500 FloatComplexMatrix retval (a.rows () + 1, a.columns ());
501 retval.assign (idx_vector (0, i), idx_vector::colon,
502 a.index (idx_vector (0, i), idx_vector::colon));
503 retval.assign (idx_vector (i), idx_vector::colon, x);
504 retval.assign (idx_vector (i+1, retval.rows ()), idx_vector::colon,
505 a.index (idx_vector (i, a.rows ()), idx_vector::colon));
506 return retval;
507 }
508
509 static
510 FloatComplexMatrix delete_col (const FloatComplexMatrix& a, octave_idx_type i)
511 {
512 FloatComplexMatrix retval = a;
513 retval.delete_elements (1, idx_vector (i));
514 return retval;
515 }
516
517 static
518 FloatComplexMatrix delete_row (const FloatComplexMatrix& a, octave_idx_type i)
519 {
520 FloatComplexMatrix retval = a;
521 retval.delete_elements (0, idx_vector (i));
522 return retval;
523 }
524
525 static
526 FloatComplexMatrix shift_cols (const FloatComplexMatrix& a,
527 octave_idx_type i, octave_idx_type j)
528 {
529 octave_idx_type n = a.columns ();
530 Array<octave_idx_type> p (n);
531 for (octave_idx_type k = 0; k < n; k++) p(k) = k;
532 if (i < j)
533 {
534 for (octave_idx_type k = i; k < j; k++) p(k) = k+1;
535 p(j) = i;
536 }
537 else if (j < i)
538 {
539 p(j) = i;
540 for (octave_idx_type k = j+1; k < i+1; k++) p(k) = k-1;
541 }
542
543 return a.index (idx_vector::colon, idx_vector (p));
544 }
545
546 void
547 FloatComplexQR::insert_col (const FloatComplexColumnVector& u, octave_idx_type j)
548 {
549 warn_qrupdate_once ();
550
551 octave_idx_type m = q.rows ();
552 octave_idx_type n = r.columns ();
553
554 if (u.length () != m)
555 (*current_liboctave_error_handler) ("qrinsert: dimensions mismatch");
556 else if (j < 0 || j > n)
557 (*current_liboctave_error_handler) ("qrinsert: index out of range");
558 else
559 {
560 init (::insert_col (q*r, j, u), get_type ());
561 }
562 }
563
564 void
565 FloatComplexQR::insert_col (const FloatComplexMatrix& u, const Array<octave_idx_type>& j)
566 {
567 warn_qrupdate_once ();
568
569 octave_idx_type m = q.rows ();
570 octave_idx_type n = r.columns ();
571
572 Array<octave_idx_type> jsi;
573 Array<octave_idx_type> js = j.sort (jsi, ASCENDING);
574 octave_idx_type nj = js.length ();
575 bool dups = false;
576 for (octave_idx_type i = 0; i < nj - 1; i++)
577 dups = dups && js(i) == js(i+1);
578
579 if (dups)
580 (*current_liboctave_error_handler) ("qrinsert: duplicate index detected");
581 else if (u.length () != m || u.columns () != nj)
582 (*current_liboctave_error_handler) ("qrinsert: dimensions mismatch");
583 else if (nj > 0 && (js(0) < 0 || js(nj-1) > n))
584 (*current_liboctave_error_handler) ("qrinsert: index out of range");
585 else if (nj > 0)
586 {
587 FloatComplexMatrix a = q*r;
588 for (octave_idx_type i = 0; i < js.length (); i++)
589 a = ::insert_col (a, js(i), u.column (i));
590 init (a, get_type ());
591 }
592 }
593
594 void
595 FloatComplexQR::delete_col (octave_idx_type j)
596 {
597 warn_qrupdate_once ();
598
599 octave_idx_type m = q.rows ();
600 octave_idx_type n = r.columns ();
601
602 if (j < 0 || j > n-1)
603 (*current_liboctave_error_handler) ("qrdelete: index out of range");
604 else
605 {
606 init (::delete_col (q*r, j), get_type ());
607 }
608 }
609
610 void
611 FloatComplexQR::delete_col (const Array<octave_idx_type>& j)
612 {
613 warn_qrupdate_once ();
614
615 octave_idx_type m = q.rows ();
616 octave_idx_type n = r.columns ();
617
618 Array<octave_idx_type> jsi;
619 Array<octave_idx_type> js = j.sort (jsi, DESCENDING);
620 octave_idx_type nj = js.length ();
621 bool dups = false;
622 for (octave_idx_type i = 0; i < nj - 1; i++)
623 dups = dups && js(i) == js(i+1);
624
625 if (dups)
626 (*current_liboctave_error_handler) ("qrinsert: duplicate index detected");
627 else if (nj > 0 && (js(0) > n-1 || js(nj-1) < 0))
628 (*current_liboctave_error_handler) ("qrinsert: index out of range");
629 else if (nj > 0)
630 {
631 FloatComplexMatrix a = q*r;
632 for (octave_idx_type i = 0; i < js.length (); i++)
633 a = ::delete_col (a, js(i));
634 init (a, get_type ());
635 }
636 }
637
638 void
639 FloatComplexQR::insert_row (const FloatComplexRowVector& u, octave_idx_type j)
640 {
641 warn_qrupdate_once ();
642
643 octave_idx_type m = r.rows ();
644 octave_idx_type n = r.columns ();
645
646 if (! q.is_square () || u.length () != n)
647 (*current_liboctave_error_handler) ("qrinsert: dimensions mismatch");
648 else if (j < 0 || j > m)
649 (*current_liboctave_error_handler) ("qrinsert: index out of range");
650 else
651 {
652 init (::insert_row (q*r, j, u), get_type ());
653 }
654 }
655
656 void
657 FloatComplexQR::delete_row (octave_idx_type j)
658 {
659 warn_qrupdate_once ();
660
661 octave_idx_type m = r.rows ();
662 octave_idx_type n = r.columns ();
663
664 if (! q.is_square ())
665 (*current_liboctave_error_handler) ("qrdelete: dimensions mismatch");
666 else if (j < 0 || j > m-1)
667 (*current_liboctave_error_handler) ("qrdelete: index out of range");
668 else
669 {
670 init (::delete_row (q*r, j), get_type ());
671 }
672 }
673
674 void
675 FloatComplexQR::shift_cols (octave_idx_type i, octave_idx_type j)
676 {
677 warn_qrupdate_once ();
678
679 octave_idx_type m = q.rows ();
680 octave_idx_type n = r.columns ();
681
682 if (i < 0 || i > n-1 || j < 0 || j > n-1)
683 (*current_liboctave_error_handler) ("qrshift: index out of range");
684 else
685 {
686 init (::shift_cols (q*r, i, j), get_type ());
687 }
688 }
689
433 #endif 690 #endif
434 691
435 /* 692 /*
436 ;;; Local Variables: *** 693 ;;; Local Variables: ***
437 ;;; mode: C++ *** 694 ;;; mode: C++ ***