Mercurial > hg > octave-lyh

comparison liboctave/CMatrix.cc @ 8336:9813c07ca946

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
make det take advantage of matrix type
author Jaroslav Hajek <highegg@gmail.com>
date Wed, 19 Nov 2008 15:26:39 +0100
parents 64cf956a109c
children e02242c54c49
comparison
equal deleted inserted replaced
8335:64cf956a109c 8336:9813c07ca946
54 #include "mx-cm-dm.h" 54 #include "mx-cm-dm.h"
55 #include "mx-dm-cm.h" 55 #include "mx-dm-cm.h"
56 #include "mx-cm-s.h" 56 #include "mx-cm-s.h"
57 #include "mx-inlines.cc" 57 #include "mx-inlines.cc"
58 #include "oct-cmplx.h" 58 #include "oct-cmplx.h"
59 #include "oct-norm.h"
59 60
60 #if defined (HAVE_FFTW3) 61 #if defined (HAVE_FFTW3)
61 #include "oct-fftw.h" 62 #include "oct-fftw.h"
62 #endif 63 #endif
63 64
1572 } 1573 }
1573 1574
1574 ComplexDET 1575 ComplexDET
1575 ComplexMatrix::determinant (octave_idx_type& info, double& rcon, int calc_cond) const 1576 ComplexMatrix::determinant (octave_idx_type& info, double& rcon, int calc_cond) const
1576 { 1577 {
1577 ComplexDET retval; 1578 MatrixType mattype (*this);
1578 1579 return determinant (mattype, info, rcon, calc_cond);
1579 octave_idx_type nr = rows (); 1580 }
1580 octave_idx_type nc = cols (); 1581
1581 1582 ComplexDET
1582 if (nr == 0 || nc == 0) 1583 ComplexMatrix::determinant (MatrixType& mattype,
1583 { 1584 octave_idx_type& info, double& rcon, int calc_cond) const
1584 retval = ComplexDET (1.0, 0); 1585 {
1585 } 1586 ComplexDET retval (1.0);
1587
1588 octave_idx_type nr = rows ();
1589 octave_idx_type nc = cols ();
1590
1591 if (nr != nc)
1592 (*current_liboctave_error_handler) ("matrix must be square");
1586 else 1593 else
1587 { 1594 {
1588 Array<octave_idx_type> ipvt (nr); 1595 int typ = mattype.type ();
1589 octave_idx_type *pipvt = ipvt.fortran_vec (); 1596
1590 1597 if (typ == MatrixType::Lower || typ == MatrixType::Upper)
1591 ComplexMatrix atmp = *this; 1598 {
1592 Complex *tmp_data = atmp.fortran_vec (); 1599 for (octave_idx_type i = 0; i < nc; i++)
1593 1600 retval *= elem (i,i);
1594 info = 0; 1601 }
1595 1602 else if (typ == MatrixType::Hermitian)
1596 // Calculate the norm of the matrix, for later use. 1603 {
1597 double anorm = 0; 1604 ComplexMatrix atmp = *this;
1598 if (calc_cond) 1605 Complex *tmp_data = atmp.fortran_vec ();
1599 anorm = atmp.abs().sum().row(static_cast<octave_idx_type>(0)).max(); 1606
1600 1607 info = 0;
1601 F77_XFCN (zgetrf, ZGETRF, (nr, nc, tmp_data, nr, pipvt, info)); 1608 double anorm = 0;
1602 1609 if (calc_cond) anorm = xnorm (*this, 1);
1603 // Throw-away extra info LAPACK gives so as to not change output. 1610
1604 rcon = 0.0; 1611
1605 if (info != 0) 1612 char job = 'L';
1606 { 1613 F77_XFCN (zpotrf, ZPOTRF, (F77_CONST_CHAR_ARG2 (&job, 1), nr,
1607 info = -1; 1614 tmp_data, nr, info
1608 retval = ComplexDET (); 1615 F77_CHAR_ARG_LEN (1)));
1609 } 1616
1610 else 1617 if (info != 0)
1611 { 1618 {
1612 if (calc_cond) 1619 rcon = 0.0;
1613 { 1620 mattype.mark_as_unsymmetric ();
1614 // Now calc the condition number for non-singular matrix. 1621 typ = MatrixType::Full;
1615 char job = '1'; 1622 }
1616 Array<Complex> z (2*nr); 1623 else
1617 Complex *pz = z.fortran_vec (); 1624 {
1618 Array<double> rz (2*nr); 1625 Array<Complex> z (2 * nc);
1619 double *prz = rz.fortran_vec (); 1626 Complex *pz = z.fortran_vec ();
1620 1627 Array<double> rz (nc);
1621 F77_XFCN (zgecon, ZGECON, (F77_CONST_CHAR_ARG2 (&job, 1), 1628 double *prz = rz.fortran_vec ();
1622 nc, tmp_data, nr, anorm, 1629
1623 rcon, pz, prz, info 1630 F77_XFCN (zpocon, ZPOCON, (F77_CONST_CHAR_ARG2 (&job, 1),
1624 F77_CHAR_ARG_LEN (1))); 1631 nr, tmp_data, nr, anorm,
1625 } 1632 rcon, pz, prz, info
1626 1633 F77_CHAR_ARG_LEN (1)));
1627 if (info != 0) 1634
1628 { 1635 if (info != 0)
1629 info = -1; 1636 rcon = 0.0;
1630 retval = ComplexDET (); 1637
1631 } 1638 for (octave_idx_type i = 0; i < nc; i++)
1632 else 1639 retval *= elem (i,i);
1633 { 1640
1634 retval = ComplexDET (1.0); 1641 retval = retval.square ();
1635 1642 }
1636 for (octave_idx_type i = 0; i < nc; i++) 1643 }
1637 { 1644 else if (typ != MatrixType::Full)
1638 Complex c = atmp(i,i); 1645 (*current_liboctave_error_handler) ("det: invalid dense matrix type");
1639 retval *= (ipvt(i) != (i+1)) ? -c : c; 1646
1647 if (typ == MatrixType::Full)
1648 {
1649 Array<octave_idx_type> ipvt (nr);
1650 octave_idx_type *pipvt = ipvt.fortran_vec ();
1651
1652 ComplexMatrix atmp = *this;
1653 Complex *tmp_data = atmp.fortran_vec ();
1654
1655 info = 0;
1656
1657 // Calculate the norm of the matrix, for later use.
1658 double anorm = 0;
1659 if (calc_cond) anorm = xnorm (*this, 1);
1660
1661 F77_XFCN (zgetrf, ZGETRF, (nr, nr, tmp_data, nr, pipvt, info));
1662
1663 // Throw-away extra info LAPACK gives so as to not change output.
1664 rcon = 0.0;
1665 if (info != 0)
1666 {
1667 info = -1;
1668 retval = ComplexDET ();
1669 }
1670 else
1671 {
1672 if (calc_cond)
1673 {
1674 // Now calc the condition number for non-singular matrix.
1675 char job = '1';
1676 Array<Complex> z (2 * nc);
1677 Complex *pz = z.fortran_vec ();
1678 Array<double> rz (2 * nc);
1679 double *prz = rz.fortran_vec ();
1680
1681 F77_XFCN (zgecon, ZGECON, (F77_CONST_CHAR_ARG2 (&job, 1),
1682 nc, tmp_data, nr, anorm,
1683 rcon, pz, prz, info
1684 F77_CHAR_ARG_LEN (1)));
1640 } 1685 }
1641 } 1686
1642 } 1687 if (info != 0)
1643 } 1688 {
1644 1689 info = -1;
1690 retval = ComplexDET ();
1691 }
1692 else
1693 {
1694 for (octave_idx_type i = 0; i < nc; i++)
1695 {
1696 Complex c = atmp(i,i);
1697 retval *= (ipvt(i) != (i+1)) ? -c : c;
1698 }
1699 }
1700 }
1701 }
1702 }
1703
1645 return retval; 1704 return retval;
1646 } 1705 }
1647 1706
1648 double 1707 double
1649 ComplexMatrix::rcond (void) const 1708 ComplexMatrix::rcond (void) const