Mercurial > octave

comparison liboctave/numeric/chol.cc @ 21270:230e186e292d

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
make building without qrupdate work again * dMatrix.h (Matrix::hermitian): New function. * fMatrix.h (FloatMatrix::hermitian): New function. * liboctave/numeric/chol.cc: Fix function declarations, definition of zero, and names of imag and conj functions.
author John W. Eaton <jwe@octave.org>
date Tue, 16 Feb 2016 11:13:55 -0500
parents 3c8a3d35661a
children d3b265a83adc
comparison
equal deleted inserted replaced
21269:3c8a3d35661a 21270:230e186e292d
25 # include <config.h> 25 # include <config.h>
26 #endif 26 #endif
27 27
28 #include <vector> 28 #include <vector>
29 29
30
31 #include "CColVector.h" 30 #include "CColVector.h"
32 #include "CMatrix.h" 31 #include "CMatrix.h"
32 #include "CRowVector.h"
33 #include "chol.h" 33 #include "chol.h"
34 #include "dColVector.h" 34 #include "dColVector.h"
35 #include "dMatrix.h" 35 #include "dMatrix.h"
36 #include "dRowVector.h"
36 #include "f77-fcn.h" 37 #include "f77-fcn.h"
37 #include "fCColVector.h" 38 #include "fCColVector.h"
38 #include "fCMatrix.h" 39 #include "fCMatrix.h"
40 #include "fCRowVector.h"
39 #include "fColVector.h" 41 #include "fColVector.h"
40 #include "fMatrix.h" 42 #include "fMatrix.h"
43 #include "fRowVector.h"
41 #include "lo-error.h" 44 #include "lo-error.h"
42 #include "oct-locbuf.h" 45 #include "oct-locbuf.h"
43 #include "oct-norm.h" 46 #include "oct-norm.h"
44 47
45 #if ! defined (HAVE_QRUPDATE) 48 #if ! defined (HAVE_QRUPDATE)
447 450
448 #if ! defined (HAVE_QRUPDATE) 451 #if ! defined (HAVE_QRUPDATE)
449 452
450 template <typename T> 453 template <typename T>
451 void 454 void
452 chol<T>::update (const T::VT& u) 455 chol<T>::update (const VT& u)
453 { 456 {
454 warn_qrupdate_once (); 457 warn_qrupdate_once ();
455 458
456 octave_idx_type n = chol_mat.rows (); 459 octave_idx_type n = chol_mat.rows ();
457 460
464 467
465 template <typename T> 468 template <typename T>
466 static bool 469 static bool
467 singular (const T& a) 470 singular (const T& a)
468 { 471 {
469 static typename T::element_type zero (); 472 static typename T::element_type zero (0);
470 for (octave_idx_type i = 0; i < a.rows (); i++) 473 for (octave_idx_type i = 0; i < a.rows (); i++)
471 if (a(i,i) == zero) return true; 474 if (a(i,i) == zero) return true;
472 return false; 475 return false;
473 } 476 }
474 477
475 template <typename T> 478 template <typename T>
476 octave_idx_type 479 octave_idx_type
477 chol<T>::downdate (const T::VT& u) 480 chol<T>::downdate (const VT& u)
478 { 481 {
479 warn_qrupdate_once (); 482 warn_qrupdate_once ();
480 483
481 octave_idx_type info = -1; 484 octave_idx_type info = -1;
482 485
497 return info; 500 return info;
498 } 501 }
499 502
500 template <typename T> 503 template <typename T>
501 octave_idx_type 504 octave_idx_type
502 chol<T>::insert_sym (const T::VT& u, octave_idx_type j) 505 chol<T>::insert_sym (const VT& u, octave_idx_type j)
503 { 506 {
504 static typename T::element_type zero (); 507 static typename T::element_type zero (0);
505 508
506 warn_qrupdate_once (); 509 warn_qrupdate_once ();
507 510
508 octave_idx_type info = -1; 511 octave_idx_type info = -1;
509 512
514 if (j < 0 || j > n) 517 if (j < 0 || j > n)
515 (*current_liboctave_error_handler) ("cholinsert: index out of range"); 518 (*current_liboctave_error_handler) ("cholinsert: index out of range");
516 519
517 if (singular (chol_mat)) 520 if (singular (chol_mat))
518 info = 2; 521 info = 2;
519 else if (ximag (u(j)) != zero) 522 else if (imag (u(j)) != zero)
520 info = 3; 523 info = 3;
521 else 524 else
522 { 525 {
523 T a = chol_mat.hermitian () * chol_mat; 526 T a = chol_mat.hermitian () * chol_mat;
524 T a1 (n+1, n+1); 527 T a1 (n+1, n+1);
526 for (octave_idx_type l = 0; l < n+1; l++) 529 for (octave_idx_type l = 0; l < n+1; l++)
527 { 530 {
528 if (l == j) 531 if (l == j)
529 a1(k, l) = u(k); 532 a1(k, l) = u(k);
530 else if (k == j) 533 else if (k == j)
531 a1(k, l) = xconj (u(l)); 534 a1(k, l) = conj (u(l));
532 else 535 else
533 a1(k, l) = a(k < j ? k : k-1, l < j ? l : l-1); 536 a1(k, l) = a(k < j ? k : k-1, l < j ? l : l-1);
534 } 537 }
535 info = init (a1, true, false); 538 info = init (a1, true, false);
536 if (info) info = 1; 539 if (info) info = 1;