Mercurial > octave-nkf

comparison liboctave/fCMatrix.cc @ 8335:64cf956a109c

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
templatize & fix DET
author Jaroslav Hajek <highegg@gmail.com>
date Wed, 19 Nov 2008 11:23:07 +0100
parents 851803f7bb4d
children 9813c07ca946
comparison
equal deleted inserted replaced
8334:70dd33450061 8335:64cf956a109c
36 #include <sys/types.h> 36 #include <sys/types.h>
37 #endif 37 #endif
38 38
39 #include "Array-util.h" 39 #include "Array-util.h"
40 #include "fCMatrix.h" 40 #include "fCMatrix.h"
41 #include "fCmplxDET.h" 41 #include "DET.h"
42 #include "fCmplxSCHUR.h" 42 #include "fCmplxSCHUR.h"
43 #include "fCmplxSVD.h" 43 #include "fCmplxSVD.h"
44 #include "fCmplxCHOL.h" 44 #include "fCmplxCHOL.h"
45 #include "f77-fcn.h" 45 #include "f77-fcn.h"
46 #include "functor.h" 46 #include "functor.h"
1622 info = -1; 1622 info = -1;
1623 retval = FloatComplexDET (); 1623 retval = FloatComplexDET ();
1624 } 1624 }
1625 else 1625 else
1626 { 1626 {
1627 FloatComplex c = 1.0; 1627 retval = FloatComplexDET (1.0);
1628 int e = 0; 1628
1629
1630 for (octave_idx_type i = 0; i < nc; i++) 1629 for (octave_idx_type i = 0; i < nc; i++)
1631 { 1630 {
1632 if (ipvt(i) != (i+1)) 1631 FloatComplex c = atmp(i,i);
1633 c = -c; 1632 retval *= (ipvt(i) != (i+1)) ? -c : c;
1634 1633 }
1635 c *= atmp(i,i);
1636
1637 if (c == static_cast<float> (0.0))
1638 break;
1639
1640 while (std::abs(c) < 0.5)
1641 {
1642 c *= 2.0;
1643 e--;
1644 }
1645
1646 while (std::abs(c) >= 2.0)
1647 {
1648 c /= 2.0;
1649 e++;
1650 }
1651 }
1652
1653 retval = FloatComplexDET (c, e);
1654 } 1634 }
1655 } 1635 }
1656 } 1636 }
1657 1637
1658 return retval; 1638 return retval;