Mercurial > octave-nkf

comparison liboctave/CMatrix.cc @ 5785:6b9cec830d72

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
[project @ 2006-05-03 19:32:46 by dbateman]
author dbateman
date Wed, 03 May 2006 19:32:48 +0000
parents ace8d8d26933
children cdef72fcd206
comparison
equal deleted inserted replaced
5784:70f7659d0fb9 5785:6b9cec830d72
109 F77_RET_T 109 F77_RET_T
110 F77_FUNC (zgelss, ZGELSS) (const octave_idx_type&, const octave_idx_type&, const octave_idx_type&, 110 F77_FUNC (zgelss, ZGELSS) (const octave_idx_type&, const octave_idx_type&, const octave_idx_type&,
111 Complex*, const octave_idx_type&, Complex*, 111 Complex*, const octave_idx_type&, Complex*,
112 const octave_idx_type&, double*, double&, octave_idx_type&, 112 const octave_idx_type&, double*, double&, octave_idx_type&,
113 Complex*, const octave_idx_type&, double*, octave_idx_type&); 113 Complex*, const octave_idx_type&, double*, octave_idx_type&);
114
115 F77_RET_T
116 F77_FUNC (zpotrf, ZPOTRF) (F77_CONST_CHAR_ARG_DECL, const octave_idx_type&,
117 Complex*, const octave_idx_type&,
118 octave_idx_type& F77_CHAR_ARG_LEN_DECL);
119
120 F77_RET_T
121 F77_FUNC (zpocon, ZPOCON) (F77_CONST_CHAR_ARG_DECL, const octave_idx_type&,
122 Complex*, const octave_idx_type&, const double&,
123 double&, Complex*, double*,
124 octave_idx_type& F77_CHAR_ARG_LEN_DECL);
125
126 F77_RET_T
127 F77_FUNC (zpotrs, ZPOTRS) (F77_CONST_CHAR_ARG_DECL, const octave_idx_type&,
128 const octave_idx_type&, const Complex*,
129 const octave_idx_type&, Complex*,
130 const octave_idx_type&, octave_idx_type&
131 F77_CHAR_ARG_LEN_DECL);
132
133 F77_RET_T
134 F77_FUNC (ztrcon, ZTRCON) (F77_CONST_CHAR_ARG_DECL, F77_CONST_CHAR_ARG_DECL,
135 F77_CONST_CHAR_ARG_DECL, const octave_idx_type&,
136 const Complex*, const octave_idx_type&, double&,
137 Complex*, double*, octave_idx_type&
138 F77_CHAR_ARG_LEN_DECL
139 F77_CHAR_ARG_LEN_DECL
140 F77_CHAR_ARG_LEN_DECL);
141
142 F77_RET_T
143 F77_FUNC (ztrtrs, ZTRTRS) (F77_CONST_CHAR_ARG_DECL, F77_CONST_CHAR_ARG_DECL,
144 F77_CONST_CHAR_ARG_DECL, const octave_idx_type&,
145 const octave_idx_type&, const Complex*,
146 const octave_idx_type&, Complex*,
147 const octave_idx_type&, octave_idx_type&
148 F77_CHAR_ARG_LEN_DECL
149 F77_CHAR_ARG_LEN_DECL
150 F77_CHAR_ARG_LEN_DECL);
114 151
115 // Note that the original complex fft routines were not written for 152 // Note that the original complex fft routines were not written for
116 // double complex arguments. They have been modified by adding an 153 // double complex arguments. They have been modified by adding an
117 // implicit double precision (a-h,o-z) statement at the beginning of 154 // implicit double precision (a-h,o-z) statement at the beginning of
118 // each subroutine. 155 // each subroutine.
1489 1526
1490 return retval; 1527 return retval;
1491 } 1528 }
1492 1529
1493 ComplexMatrix 1530 ComplexMatrix
1494 ComplexMatrix::solve (const Matrix& b) const 1531 ComplexMatrix::utsolve (MatrixType &mattype, const ComplexMatrix& b,
1495 { 1532 octave_idx_type& info, double& rcond,
1496 octave_idx_type info; 1533 solve_singularity_handler sing_handler,
1497 double rcond; 1534 bool calc_cond) const
1498 return solve (b, info, rcond, 0); 1535 {
1499 } 1536 ComplexMatrix retval;
1500 1537
1501 ComplexMatrix 1538 octave_idx_type nr = rows ();
1502 ComplexMatrix::solve (const Matrix& b, octave_idx_type& info) const 1539 octave_idx_type nc = cols ();
1503 { 1540
1504 double rcond; 1541 if (nr == 0 || nc == 0 || nr != b.rows ())
1505 return solve (b, info, rcond, 0); 1542 (*current_liboctave_error_handler)
1506 } 1543 ("matrix dimension mismatch solution of linear equations");
1507 1544 else
1508 ComplexMatrix 1545 {
1509 ComplexMatrix::solve (const Matrix& b, octave_idx_type& info, double& rcond) const 1546 volatile int typ = mattype.type ();
1510 { 1547
1511 return solve (b, info, rcond, 0); 1548 if (typ == MatrixType::Permuted_Upper ||
1512 } 1549 typ == MatrixType::Upper)
1513 1550 {
1514 ComplexMatrix 1551 octave_idx_type b_nc = b.cols ();
1515 ComplexMatrix::solve (const Matrix& b, octave_idx_type& info, double& rcond, 1552 rcond = 1.;
1516 solve_singularity_handler sing_handler) const 1553 info = 0;
1517 { 1554
1518 ComplexMatrix tmp (b); 1555 if (typ == MatrixType::Permuted_Upper)
1519 return solve (tmp, info, rcond, sing_handler); 1556 {
1520 } 1557 (*current_liboctave_error_handler)
1521 1558 ("Permuted triangular matrix not implemented");
1522 ComplexMatrix 1559 }
1523 ComplexMatrix::solve (const ComplexMatrix& b) const 1560 else
1524 { 1561 {
1525 octave_idx_type info; 1562 const Complex *tmp_data = fortran_vec ();
1526 double rcond; 1563
1527 return solve (b, info, rcond, 0); 1564 if (calc_cond)
1528 } 1565 {
1529 1566 char norm = '1';
1530 ComplexMatrix 1567 char uplo = 'U';
1531 ComplexMatrix::solve (const ComplexMatrix& b, octave_idx_type& info) const 1568 char dia = 'N';
1532 { 1569
1533 double rcond; 1570 Array<Complex> z (2 * nc);
1534 return solve (b, info, rcond, 0); 1571 Complex *pz = z.fortran_vec ();
1535 } 1572 Array<double> rz (nc);
1536 1573 double *prz = rz.fortran_vec ();
1537 ComplexMatrix 1574
1538 ComplexMatrix::solve (const ComplexMatrix& b, octave_idx_type& info, double& rcond) const 1575 F77_XFCN (ztrcon, ZTRCON, (F77_CONST_CHAR_ARG2 (&norm, 1),
1539 { 1576 F77_CONST_CHAR_ARG2 (&uplo, 1),
1540 return solve (b, info, rcond, 0); 1577 F77_CONST_CHAR_ARG2 (&dia, 1),
1541 } 1578 nr, tmp_data, nr, rcond,
1542 1579 pz, prz, info
1543 ComplexMatrix 1580 F77_CHAR_ARG_LEN (1)
1544 ComplexMatrix::solve (const ComplexMatrix& b, octave_idx_type& info, double& rcond, 1581 F77_CHAR_ARG_LEN (1)
1545 solve_singularity_handler sing_handler) const 1582 F77_CHAR_ARG_LEN (1)));
1583
1584 if (f77_exception_encountered)
1585 (*current_liboctave_error_handler)
1586 ("unrecoverable error in ztrcon");
1587
1588 if (info != 0)
1589 info = -2;
1590
1591 volatile double rcond_plus_one = rcond + 1.0;
1592
1593 if (rcond_plus_one == 1.0 || xisnan (rcond))
1594 {
1595 info = -2;
1596
1597 if (sing_handler)
1598 sing_handler (rcond);
1599 else
1600 (*current_liboctave_error_handler)
1601 ("matrix singular to machine precision, rcond = %g",
1602 rcond);
1603 }
1604 }
1605
1606 if (info == 0)
1607 {
1608 retval = b;
1609 Complex *result = retval.fortran_vec ();
1610
1611 char uplo = 'U';
1612 char trans = 'N';
1613 char dia = 'N';
1614
1615 F77_XFCN (ztrtrs, ZTRTRS, (F77_CONST_CHAR_ARG2 (&uplo, 1),
1616 F77_CONST_CHAR_ARG2 (&trans, 1),
1617 F77_CONST_CHAR_ARG2 (&dia, 1),
1618 nr, b_nc, tmp_data, nr,
1619 result, nr, info
1620 F77_CHAR_ARG_LEN (1)
1621 F77_CHAR_ARG_LEN (1)
1622 F77_CHAR_ARG_LEN (1)));
1623
1624 if (f77_exception_encountered)
1625 (*current_liboctave_error_handler)
1626 ("unrecoverable error in dtrtrs");
1627 }
1628 }
1629 }
1630 else
1631 (*current_liboctave_error_handler) ("incorrect matrix type");
1632 }
1633
1634 return retval;
1635 }
1636
1637 ComplexMatrix
1638 ComplexMatrix::ltsolve (MatrixType &mattype, const ComplexMatrix& b,
1639 octave_idx_type& info, double& rcond,
1640 solve_singularity_handler sing_handler,
1641 bool calc_cond) const
1642 {
1643 ComplexMatrix retval;
1644
1645 octave_idx_type nr = rows ();
1646 octave_idx_type nc = cols ();
1647
1648 if (nr == 0 || nc == 0 || nr != b.rows ())
1649 (*current_liboctave_error_handler)
1650 ("matrix dimension mismatch solution of linear equations");
1651 else
1652 {
1653 volatile int typ = mattype.type ();
1654
1655 if (typ == MatrixType::Permuted_Lower ||
1656 typ == MatrixType::Lower)
1657 {
1658 octave_idx_type b_nc = b.cols ();
1659 rcond = 1.;
1660 info = 0;
1661
1662 if (typ == MatrixType::Permuted_Lower)
1663 {
1664 (*current_liboctave_error_handler)
1665 ("Permuted triangular matrix not implemented");
1666 }
1667 else
1668 {
1669 const Complex *tmp_data = fortran_vec ();
1670
1671 if (calc_cond)
1672 {
1673 char norm = '1';
1674 char uplo = 'L';
1675 char dia = 'N';
1676
1677 Array<Complex> z (2 * nc);
1678 Complex *pz = z.fortran_vec ();
1679 Array<double> rz (nc);
1680 double *prz = rz.fortran_vec ();
1681
1682 F77_XFCN (ztrcon, ZTRCON, (F77_CONST_CHAR_ARG2 (&norm, 1),
1683 F77_CONST_CHAR_ARG2 (&uplo, 1),
1684 F77_CONST_CHAR_ARG2 (&dia, 1),
1685 nr, tmp_data, nr, rcond,
1686 pz, prz, info
1687 F77_CHAR_ARG_LEN (1)
1688 F77_CHAR_ARG_LEN (1)
1689 F77_CHAR_ARG_LEN (1)));
1690
1691 if (f77_exception_encountered)
1692 (*current_liboctave_error_handler)
1693 ("unrecoverable error in ztrcon");
1694
1695 if (info != 0)
1696 info = -2;
1697
1698 volatile double rcond_plus_one = rcond + 1.0;
1699
1700 if (rcond_plus_one == 1.0 || xisnan (rcond))
1701 {
1702 info = -2;
1703
1704 if (sing_handler)
1705 sing_handler (rcond);
1706 else
1707 (*current_liboctave_error_handler)
1708 ("matrix singular to machine precision, rcond = %g",
1709 rcond);
1710 }
1711 }
1712
1713 if (info == 0)
1714 {
1715 retval = b;
1716 Complex *result = retval.fortran_vec ();
1717
1718 char uplo = 'L';
1719 char trans = 'N';
1720 char dia = 'N';
1721
1722 F77_XFCN (ztrtrs, ZTRTRS, (F77_CONST_CHAR_ARG2 (&uplo, 1),
1723 F77_CONST_CHAR_ARG2 (&trans, 1),
1724 F77_CONST_CHAR_ARG2 (&dia, 1),
1725 nr, b_nc, tmp_data, nr,
1726 result, nr, info
1727 F77_CHAR_ARG_LEN (1)
1728 F77_CHAR_ARG_LEN (1)
1729 F77_CHAR_ARG_LEN (1)));
1730
1731 if (f77_exception_encountered)
1732 (*current_liboctave_error_handler)
1733 ("unrecoverable error in dtrtrs");
1734 }
1735 }
1736 }
1737 else
1738 (*current_liboctave_error_handler) ("incorrect matrix type");
1739 }
1740
1741 return retval;
1742 }
1743
1744 ComplexMatrix
1745 ComplexMatrix::fsolve (MatrixType &mattype, const ComplexMatrix& b,
1746 octave_idx_type& info, double& rcond,
1747 solve_singularity_handler sing_handler,
1748 bool calc_cond) const
1546 { 1749 {
1547 ComplexMatrix retval; 1750 ComplexMatrix retval;
1548 1751
1549 octave_idx_type nr = rows (); 1752 octave_idx_type nr = rows ();
1550 octave_idx_type nc = cols (); 1753 octave_idx_type nc = cols ();
1552 if (nr == 0 || nc == 0 || nr != nc || nr != b.rows ()) 1755 if (nr == 0 || nc == 0 || nr != nc || nr != b.rows ())
1553 (*current_liboctave_error_handler) 1756 (*current_liboctave_error_handler)
1554 ("matrix dimension mismatch in solution of linear equations"); 1757 ("matrix dimension mismatch in solution of linear equations");
1555 else 1758 else
1556 { 1759 {
1557 info = 0; 1760 volatile int typ = mattype.type ();
1558 1761
1559 Array<octave_idx_type> ipvt (nr); 1762 // Calculate the norm of the matrix, for later use.
1560 octave_idx_type *pipvt = ipvt.fortran_vec (); 1763 double anorm = -1.;
1561 1764
1562 ComplexMatrix atmp = *this; 1765 if (typ == MatrixType::Hermitian)
1563 Complex *tmp_data = atmp.fortran_vec ();
1564
1565 Array<Complex> z (2 * nc);
1566 Complex *pz = z.fortran_vec ();
1567 Array<double> rz (2 * nc);
1568 double *prz = rz.fortran_vec ();
1569
1570 // Calculate the norm of the matrix, for later use.
1571 double anorm = atmp.abs().sum().row(static_cast<octave_idx_type>(0)).max();
1572
1573 F77_XFCN (zgetrf, ZGETRF, (nr, nr, tmp_data, nr, pipvt, info));
1574
1575 if (f77_exception_encountered)
1576 (*current_liboctave_error_handler) ("unrecoverable error in zgetrf");
1577 else
1578 { 1766 {
1579 // Throw-away extra info LAPACK gives so as to not change output. 1767 info = 0;
1580 rcond = 0.0; 1768 char job = 'L';
1581 if (info != 0) 1769 ComplexMatrix atmp = *this;
1582 { 1770 Complex *tmp_data = atmp.fortran_vec ();
1583 info = -2; 1771 anorm = atmp.abs().sum().row(static_cast<octave_idx_type>(0)).max();
1584 1772
1585 if (sing_handler) 1773 F77_XFCN (zpotrf, ZPOTRF, (F77_CONST_CHAR_ARG2 (&job, 1), nr,
1586 sing_handler (rcond); 1774 tmp_data, nr, info
1587 else 1775 F77_CHAR_ARG_LEN (1)));
1588 (*current_liboctave_error_handler) 1776
1589 ("matrix singular to machine precision"); 1777 if (f77_exception_encountered)
1590 1778 (*current_liboctave_error_handler)
1591 } 1779 ("unrecoverable error in zpotrf");
1592 else 1780 else
1593 { 1781 {
1594 // Now calculate the condition number for non-singular matrix. 1782 // Throw-away extra info LAPACK gives so as to not change output.
1595 char job = '1'; 1783 rcond = 0.0;
1596 F77_XFCN (zgecon, ZGECON, (F77_CONST_CHAR_ARG2 (&job, 1),
1597 nc, tmp_data, nr, anorm,
1598 rcond, pz, prz, info
1599 F77_CHAR_ARG_LEN (1)));
1600
1601 if (f77_exception_encountered)
1602 (*current_liboctave_error_handler)
1603 ("unrecoverable error in zgecon");
1604
1605 if (info != 0) 1784 if (info != 0)
1606 info = -2;
1607
1608 volatile double rcond_plus_one = rcond + 1.0;
1609
1610 if (rcond_plus_one == 1.0 || xisnan (rcond))
1611 { 1785 {
1786 info = -2;
1787
1788 mattype.mark_as_unsymmetric ();
1789 typ = MatrixType::Full;
1790 }
1791 else
1792 {
1793 if (calc_cond)
1794 {
1795 Array<Complex> z (2 * nc);
1796 Complex *pz = z.fortran_vec ();
1797 Array<double> rz (nc);
1798 double *prz = rz.fortran_vec ();
1799
1800 F77_XFCN (zpocon, ZPOCON, (F77_CONST_CHAR_ARG2 (&job, 1),
1801 nr, tmp_data, nr, anorm,
1802 rcond, pz, prz, info
1803 F77_CHAR_ARG_LEN (1)));
1804
1805 if (f77_exception_encountered)
1806 (*current_liboctave_error_handler)
1807 ("unrecoverable error in zpocon");
1808
1809 if (info != 0)
1810 info = -2;
1811
1812 volatile double rcond_plus_one = rcond + 1.0;
1813
1814 if (rcond_plus_one == 1.0 || xisnan (rcond))
1815 {
1816 info = -2;
1817
1818 if (sing_handler)
1819 sing_handler (rcond);
1820 else
1821 (*current_liboctave_error_handler)
1822 ("matrix singular to machine precision, rcond = %g",
1823 rcond);
1824 }
1825 }
1826
1827 if (info == 0)
1828 {
1829 retval = b;
1830 Complex *result = retval.fortran_vec ();
1831
1832 octave_idx_type b_nc = b.cols ();
1833
1834 F77_XFCN (zpotrs, ZPOTRS, (F77_CONST_CHAR_ARG2 (&job, 1),
1835 nr, b_nc, tmp_data, nr,
1836 result, b.rows(), info
1837 F77_CHAR_ARG_LEN (1)));
1838
1839 if (f77_exception_encountered)
1840 (*current_liboctave_error_handler)
1841 ("unrecoverable error in zpotrs");
1842 }
1843 else
1844 {
1845 mattype.mark_as_unsymmetric ();
1846 typ = MatrixType::Full;
1847 }
1848 }
1849 }
1850 }
1851
1852 if (typ == MatrixType::Full)
1853 {
1854 info = 0;
1855
1856 Array<octave_idx_type> ipvt (nr);
1857 octave_idx_type *pipvt = ipvt.fortran_vec ();
1858
1859 ComplexMatrix atmp = *this;
1860 Complex *tmp_data = atmp.fortran_vec ();
1861
1862 Array<Complex> z (2 * nc);
1863 Complex *pz = z.fortran_vec ();
1864 Array<double> rz (2 * nc);
1865 double *prz = rz.fortran_vec ();
1866
1867 // Calculate the norm of the matrix, for later use.
1868 if (anorm < 0.)
1869 anorm = atmp.abs().sum().row(static_cast<octave_idx_type>(0)).max();
1870
1871 F77_XFCN (zgetrf, ZGETRF, (nr, nr, tmp_data, nr, pipvt, info));
1872
1873 if (f77_exception_encountered)
1874 (*current_liboctave_error_handler)
1875 ("unrecoverable error in zgetrf");
1876 else
1877 {
1878 // Throw-away extra info LAPACK gives so as to not change output.
1879 rcond = 0.0;
1880 if (info != 0)
1881 {
1612 info = -2; 1882 info = -2;
1613 1883
1614 if (sing_handler) 1884 if (sing_handler)
1615 sing_handler (rcond); 1885 sing_handler (rcond);
1616 else 1886 else
1617 (*current_liboctave_error_handler) 1887 (*current_liboctave_error_handler)
1618 ("matrix singular to machine precision, rcond = %g", 1888 ("matrix singular to machine precision");
1619 rcond); 1889
1620 } 1890 mattype.mark_as_rectangular ();
1621 else 1891 }
1892 else
1622 { 1893 {
1623 retval = b; 1894 if (calc_cond)
1624 Complex *result = retval.fortran_vec (); 1895 {
1625 1896 // Now calculate the condition number for
1626 octave_idx_type b_nc = b.cols (); 1897 // non-singular matrix.
1627 1898 char job = '1';
1628 job = 'N'; 1899 F77_XFCN (zgecon, ZGECON, (F77_CONST_CHAR_ARG2 (&job, 1),
1629 F77_XFCN (zgetrs, ZGETRS, (F77_CONST_CHAR_ARG2 (&job, 1), 1900 nc, tmp_data, nr, anorm,
1630 nr, b_nc, tmp_data, nr, 1901 rcond, pz, prz, info
1631 pipvt, result, b.rows(), info 1902 F77_CHAR_ARG_LEN (1)));
1632 F77_CHAR_ARG_LEN (1))); 1903
1633 1904 if (f77_exception_encountered)
1634 if (f77_exception_encountered) 1905 (*current_liboctave_error_handler)
1635 (*current_liboctave_error_handler) 1906 ("unrecoverable error in zgecon");
1636 ("unrecoverable error in zgetrs"); 1907
1908 if (info != 0)
1909 info = -2;
1910
1911 volatile double rcond_plus_one = rcond + 1.0;
1912
1913 if (rcond_plus_one == 1.0 || xisnan (rcond))
1914 {
1915 info = -2;
1916
1917 if (sing_handler)
1918 sing_handler (rcond);
1919 else
1920 (*current_liboctave_error_handler)
1921 ("matrix singular to machine precision, rcond = %g",
1922 rcond);
1923 }
1924 }
1925
1926 if (info == 0)
1927 {
1928 retval = b;
1929 Complex *result = retval.fortran_vec ();
1930
1931 octave_idx_type b_nc = b.cols ();
1932
1933 char job = 'N';
1934 F77_XFCN (zgetrs, ZGETRS, (F77_CONST_CHAR_ARG2 (&job, 1),
1935 nr, b_nc, tmp_data, nr,
1936 pipvt, result, b.rows(), info
1937 F77_CHAR_ARG_LEN (1)));
1938
1939 if (f77_exception_encountered)
1940 (*current_liboctave_error_handler)
1941 ("unrecoverable error in zgetrs");
1942 }
1943 else
1944 mattype.mark_as_rectangular ();
1637 } 1945 }
1638 } 1946 }
1639 } 1947 }
1640 } 1948 }
1641 1949
1642 return retval; 1950 return retval;
1643 } 1951 }
1644 1952
1953 ComplexMatrix
1954 ComplexMatrix::solve (MatrixType &typ, const Matrix& b) const
1955 {
1956 octave_idx_type info;
1957 double rcond;
1958 return solve (typ, b, info, rcond, 0);
1959 }
1960
1961 ComplexMatrix
1962 ComplexMatrix::solve (MatrixType &typ, const Matrix& b,
1963 octave_idx_type& info) const
1964 {
1965 double rcond;
1966 return solve (typ, b, info, rcond, 0);
1967 }
1968
1969 ComplexMatrix
1970 ComplexMatrix::solve (MatrixType &typ, const Matrix& b, octave_idx_type& info,
1971 double& rcond) const
1972 {
1973 return solve (typ, b, info, rcond, 0);
1974 }
1975
1976 ComplexMatrix
1977 ComplexMatrix::solve (MatrixType &typ, const Matrix& b, octave_idx_type& info,
1978 double& rcond, solve_singularity_handler sing_handler,
1979 bool singular_fallback) const
1980 {
1981 ComplexMatrix tmp (b);
1982 return solve (typ, tmp, info, rcond, sing_handler, singular_fallback);
1983 }
1984
1985 ComplexMatrix
1986 ComplexMatrix::solve (MatrixType &typ, const ComplexMatrix& b) const
1987 {
1988 octave_idx_type info;
1989 double rcond;
1990 return solve (typ, b, info, rcond, 0);
1991 }
1992
1993 ComplexMatrix
1994 ComplexMatrix::solve (MatrixType &typ, const ComplexMatrix& b,
1995 octave_idx_type& info) const
1996 {
1997 double rcond;
1998 return solve (typ, b, info, rcond, 0);
1999 }
2000
2001 ComplexMatrix
2002 ComplexMatrix::solve (MatrixType &typ, const ComplexMatrix& b,
2003 octave_idx_type& info, double& rcond) const
2004 {
2005 return solve (typ, b, info, rcond, 0);
2006 }
2007
2008 ComplexMatrix
2009 ComplexMatrix::solve (MatrixType &mattype, const ComplexMatrix& b,
2010 octave_idx_type& info, double& rcond,
2011 solve_singularity_handler sing_handler,
2012 bool singular_fallback) const
2013 {
2014 ComplexMatrix retval;
2015 int typ = mattype.type ();
2016
2017 if (typ == MatrixType::Unknown)
2018 typ = mattype.type (*this);
2019
2020 // Only calculate the condition number for LU/Cholesky
2021 if (typ == MatrixType::Upper || typ == MatrixType::Permuted_Upper)
2022 retval = utsolve (mattype, b, info, rcond, sing_handler, false);
2023 else if (typ == MatrixType::Lower || typ == MatrixType::Permuted_Lower)
2024 retval = ltsolve (mattype, b, info, rcond, sing_handler, false);
2025 else if (typ == MatrixType::Full || typ == MatrixType::Hermitian)
2026 retval = fsolve (mattype, b, info, rcond, sing_handler, true);
2027 else if (typ != MatrixType::Rectangular)
2028 {
2029 (*current_liboctave_error_handler) ("unknown matrix type");
2030 return ComplexMatrix ();
2031 }
2032
2033 // Rectangular or one of the above solvers flags a singular matrix
2034 if (singular_fallback && mattype.type () == MatrixType::Rectangular)
2035 {
2036 octave_idx_type rank;
2037 retval = lssolve (b, info, rank);
2038 }
2039
2040 return retval;
2041 }
2042
2043 ComplexColumnVector
2044 ComplexMatrix::solve (MatrixType &typ, const ColumnVector& b) const
2045 {
2046 octave_idx_type info;
2047 double rcond;
2048 return solve (typ, ComplexColumnVector (b), info, rcond, 0);
2049 }
2050
2051 ComplexColumnVector
2052 ComplexMatrix::solve (MatrixType &typ, const ColumnVector& b,
2053 octave_idx_type& info) const
2054 {
2055 double rcond;
2056 return solve (typ, ComplexColumnVector (b), info, rcond, 0);
2057 }
2058
2059 ComplexColumnVector
2060 ComplexMatrix::solve (MatrixType &typ, const ColumnVector& b,
2061 octave_idx_type& info, double& rcond) const
2062 {
2063 return solve (typ, ComplexColumnVector (b), info, rcond, 0);
2064 }
2065
2066 ComplexColumnVector
2067 ComplexMatrix::solve (MatrixType &typ, const ColumnVector& b,
2068 octave_idx_type& info, double& rcond,
2069 solve_singularity_handler sing_handler) const
2070 {
2071 return solve (typ, ComplexColumnVector (b), info, rcond, sing_handler);
2072 }
2073
2074 ComplexColumnVector
2075 ComplexMatrix::solve (MatrixType &typ, const ComplexColumnVector& b) const
2076 {
2077 octave_idx_type info;
2078 double rcond;
2079 return solve (typ, b, info, rcond, 0);
2080 }
2081
2082 ComplexColumnVector
2083 ComplexMatrix::solve (MatrixType &typ, const ComplexColumnVector& b,
2084 octave_idx_type& info) const
2085 {
2086 double rcond;
2087 return solve (typ, b, info, rcond, 0);
2088 }
2089
2090 ComplexColumnVector
2091 ComplexMatrix::solve (MatrixType &typ, const ComplexColumnVector& b,
2092 octave_idx_type& info, double& rcond) const
2093 {
2094 return solve (typ, b, info, rcond, 0);
2095 }
2096
2097 ComplexColumnVector
2098 ComplexMatrix::solve (MatrixType &typ, const ComplexColumnVector& b,
2099 octave_idx_type& info, double& rcond,
2100 solve_singularity_handler sing_handler) const
2101 {
2102
2103 ComplexMatrix tmp (b);
2104 return solve (typ, tmp, info, rcond, sing_handler).column(static_cast<octave_idx_type> (0));
2105 }
2106
2107 ComplexMatrix
2108 ComplexMatrix::solve (const Matrix& b) const
2109 {
2110 octave_idx_type info;
2111 double rcond;
2112 return solve (b, info, rcond, 0);
2113 }
2114
2115 ComplexMatrix
2116 ComplexMatrix::solve (const Matrix& b, octave_idx_type& info) const
2117 {
2118 double rcond;
2119 return solve (b, info, rcond, 0);
2120 }
2121
2122 ComplexMatrix
2123 ComplexMatrix::solve (const Matrix& b, octave_idx_type& info, double& rcond) const
2124 {
2125 return solve (b, info, rcond, 0);
2126 }
2127
2128 ComplexMatrix
2129 ComplexMatrix::solve (const Matrix& b, octave_idx_type& info, double& rcond,
2130 solve_singularity_handler sing_handler) const
2131 {
2132 ComplexMatrix tmp (b);
2133 return solve (tmp, info, rcond, sing_handler);
2134 }
2135
2136 ComplexMatrix
2137 ComplexMatrix::solve (const ComplexMatrix& b) const
2138 {
2139 octave_idx_type info;
2140 double rcond;
2141 return solve (b, info, rcond, 0);
2142 }
2143
2144 ComplexMatrix
2145 ComplexMatrix::solve (const ComplexMatrix& b, octave_idx_type& info) const
2146 {
2147 double rcond;
2148 return solve (b, info, rcond, 0);
2149 }
2150
2151 ComplexMatrix
2152 ComplexMatrix::solve (const ComplexMatrix& b, octave_idx_type& info, double& rcond) const
2153 {
2154 return solve (b, info, rcond, 0);
2155 }
2156
2157 ComplexMatrix
2158 ComplexMatrix::solve (const ComplexMatrix& b, octave_idx_type& info, double& rcond,
2159 solve_singularity_handler sing_handler) const
2160 {
2161 MatrixType mattype (*this);
2162 return solve (b, info, rcond, sing_handler);
2163 }
2164
1645 ComplexColumnVector 2165 ComplexColumnVector
1646 ComplexMatrix::solve (const ColumnVector& b) const 2166 ComplexMatrix::solve (const ColumnVector& b) const
1647 { 2167 {
1648 octave_idx_type info; 2168 octave_idx_type info;
1649 double rcond; 2169 double rcond;
1656 double rcond; 2176 double rcond;
1657 return solve (ComplexColumnVector (b), info, rcond, 0); 2177 return solve (ComplexColumnVector (b), info, rcond, 0);
1658 } 2178 }
1659 2179
1660 ComplexColumnVector 2180 ComplexColumnVector
1661 ComplexMatrix::solve (const ColumnVector& b, octave_idx_type& info, double& rcond) const 2181 ComplexMatrix::solve (const ColumnVector& b, octave_idx_type& info,
2182 double& rcond) const
1662 { 2183 {
1663 return solve (ComplexColumnVector (b), info, rcond, 0); 2184 return solve (ComplexColumnVector (b), info, rcond, 0);
1664 } 2185 }
1665 2186
1666 ComplexColumnVector 2187 ComplexColumnVector
1667 ComplexMatrix::solve (const ColumnVector& b, octave_idx_type& info, double& rcond, 2188 ComplexMatrix::solve (const ColumnVector& b, octave_idx_type& info,
2189 double& rcond,
1668 solve_singularity_handler sing_handler) const 2190 solve_singularity_handler sing_handler) const
1669 { 2191 {
1670 return solve (ComplexColumnVector (b), info, rcond, sing_handler); 2192 return solve (ComplexColumnVector (b), info, rcond, sing_handler);
1671 } 2193 }
1672 2194
1695 ComplexColumnVector 2217 ComplexColumnVector
1696 ComplexMatrix::solve (const ComplexColumnVector& b, octave_idx_type& info, 2218 ComplexMatrix::solve (const ComplexColumnVector& b, octave_idx_type& info,
1697 double& rcond, 2219 double& rcond,
1698 solve_singularity_handler sing_handler) const 2220 solve_singularity_handler sing_handler) const
1699 { 2221 {
1700 ComplexColumnVector retval; 2222 MatrixType mattype (*this);
1701 2223 return solve (mattype, b, info, rcond, sing_handler);
1702 octave_idx_type nr = rows ();
1703 octave_idx_type nc = cols ();
1704
1705 if (nr == 0 || nc == 0 || nr != nc || nr != b.length ())
1706 (*current_liboctave_error_handler)
1707 ("matrix dimension mismatch in solution of linear equations");
1708 else
1709 {
1710 info = 0;
1711
1712 Array<octave_idx_type> ipvt (nr);
1713 octave_idx_type *pipvt = ipvt.fortran_vec ();
1714
1715 ComplexMatrix atmp = *this;
1716 Complex *tmp_data = atmp.fortran_vec ();
1717
1718 Array<Complex> z (2 * nc);
1719 Complex *pz = z.fortran_vec ();
1720 Array<double> rz (2 * nc);
1721 double *prz = rz.fortran_vec ();
1722
1723 // Calculate the norm of the matrix, for later use.
1724 double anorm = atmp.abs().sum().row(static_cast<octave_idx_type>(0)).max();
1725
1726 F77_XFCN (zgetrf, ZGETRF, (nr, nr, tmp_data, nr, pipvt, info));
1727
1728 if (f77_exception_encountered)
1729 (*current_liboctave_error_handler) ("unrecoverable error in zgetrf");
1730 else
1731 {
1732 // Throw-away extra info LAPACK gives so as to not change output.
1733 rcond = 0.0;
1734 if (info != 0)
1735 {
1736 info = -2;
1737
1738 if (sing_handler)
1739 sing_handler (rcond);
1740 else
1741 (*current_liboctave_error_handler)
1742 ("matrix singular to machine precision, rcond = %g",
1743 rcond);
1744 }
1745 else
1746 {
1747 // Now calculate the condition number for non-singular matrix.
1748 char job = '1';
1749 F77_XFCN (zgecon, ZGECON, (F77_CONST_CHAR_ARG2 (&job, 1),
1750 nc, tmp_data, nr, anorm,
1751 rcond, pz, prz, info
1752 F77_CHAR_ARG_LEN (1)));
1753
1754 if (f77_exception_encountered)
1755 (*current_liboctave_error_handler)
1756 ("unrecoverable error in zgecon");
1757
1758 if (info != 0)
1759 info = -2;
1760
1761 volatile double rcond_plus_one = rcond + 1.0;
1762
1763 if (rcond_plus_one == 1.0 || xisnan (rcond))
1764 {
1765 info = -2;
1766
1767 if (sing_handler)
1768 sing_handler (rcond);
1769 else
1770 (*current_liboctave_error_handler)
1771 ("matrix singular to machine precision, rcond = %g",
1772 rcond);
1773 }
1774 else
1775 {
1776 retval = b;
1777 Complex *result = retval.fortran_vec ();
1778
1779 job = 'N';
1780 F77_XFCN (zgetrs, ZGETRS, (F77_CONST_CHAR_ARG2 (&job, 1),
1781 nr, 1, tmp_data, nr, pipvt,
1782 result, b.length(), info
1783 F77_CHAR_ARG_LEN (1)));
1784
1785 if (f77_exception_encountered)
1786 (*current_liboctave_error_handler)
1787 ("unrecoverable error in zgetrs");
1788
1789 }
1790 }
1791 }
1792 }
1793 return retval;
1794 } 2224 }
1795 2225
1796 ComplexMatrix 2226 ComplexMatrix
1797 ComplexMatrix::lssolve (const Matrix& b) const 2227 ComplexMatrix::lssolve (const Matrix& b) const
1798 { 2228 {