Mercurial > hg > octave-jordi

comparison liboctave/CMatrix.cc @ 4329:d53c33d93440

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
[project @ 2003-02-18 20:00:48 by jwe]
author jwe
date Tue, 18 Feb 2003 20:08:20 +0000
parents 236c10efcde2
children 9f86c2055b58
comparison
equal deleted inserted replaced
4328:f7b63f362168 4329:d53c33d93440
77 const Complex*, const int&, 77 const Complex*, const int&,
78 const Complex*, const int&, 78 const Complex*, const int&,
79 const Complex&, Complex*, const int&, 79 const Complex&, Complex*, const int&,
80 long, long); 80 long, long);
81 81
82 int F77_FUNC (zgeco, ZGECO) (Complex*, const int&, const int&, int*, 82 int F77_FUNC (zgetrf, ZGETRF) (const int&, const int&, Complex*, const int&,
83 double&, Complex*); 83 int*, int&);
84 84
85 int F77_FUNC (zgedi, ZGEDI) (Complex*, const int&, const int&, int*, 85 int F77_FUNC (zgetrs, ZGETRS) (const char*, const int&, const int&,
86 Complex*, Complex*, const int&); 86 Complex*, const int&,
87 87 const int*, Complex*, const int&, int&);
88 int F77_FUNC (zgesl, ZGESL) (Complex*, const int&, const int&, int*, 88
89 Complex*, const int&); 89 int F77_FUNC (zgetri, ZGETRI) (const int&, Complex*, const int&, const int*,
90 Complex*, const int&, int&);
91
92 int F77_FUNC (zgecon, ZGECON) (const char*, const int&, Complex*,
93 const int&, const double&, double&,
94 Complex*, double*, int&);
90 95
91 int F77_FUNC (zgelss, ZGELSS) (const int&, const int&, const int&, 96 int F77_FUNC (zgelss, ZGELSS) (const int&, const int&, const int&,
92 Complex*, const int&, Complex*, 97 Complex*, const int&, Complex*,
93 const int&, double*, double&, int&, 98 const int&, double*, double&, int&,
94 Complex*, const int&, double*, int&); 99 Complex*, const int&, double*, int&);
923 ComplexMatrix 928 ComplexMatrix
924 ComplexMatrix::inverse (void) const 929 ComplexMatrix::inverse (void) const
925 { 930 {
926 int info; 931 int info;
927 double rcond; 932 double rcond;
928 return inverse (info, rcond); 933 return inverse (info, rcond, 0, 0);
929 } 934 }
930 935
931 ComplexMatrix 936 ComplexMatrix
932 ComplexMatrix::inverse (int& info) const 937 ComplexMatrix::inverse (int& info) const
933 { 938 {
934 double rcond; 939 double rcond;
935 return inverse (info, rcond); 940 return inverse (info, rcond, 0, 0);
936 } 941 }
937 942
938 ComplexMatrix 943 ComplexMatrix
939 ComplexMatrix::inverse (int& info, double& rcond, int force) const 944 ComplexMatrix::inverse (int& info, double& rcond, int force,
945 int calc_cond) const
940 { 946 {
941 ComplexMatrix retval; 947 ComplexMatrix retval;
942 948
943 int nr = rows (); 949 int nr = rows ();
944 int nc = cols (); 950 int nc = cols ();
945 951
946 if (nr != nc) 952 if (nr != nc)
947 (*current_liboctave_error_handler) ("inverse requires square matrix"); 953 (*current_liboctave_error_handler) ("inverse requires square matrix");
948 else 954 else
949 { 955 {
950 info = 0;
951
952 Array<int> ipvt (nr); 956 Array<int> ipvt (nr);
953 int *pipvt = ipvt.fortran_vec (); 957 int *pipvt = ipvt.fortran_vec ();
954 958
955 Array<Complex> z (nr);
956 Complex *pz = z.fortran_vec ();
957
958 retval = *this; 959 retval = *this;
959 Complex *tmp_data = retval.fortran_vec (); 960 Complex *tmp_data = retval.fortran_vec ();
960 961
961 F77_XFCN (zgeco, ZGECO, (tmp_data, nr, nc, pipvt, rcond, pz)); 962 Array<Complex> z(1);
963 int lwork = 1;
964
965 // Query the optimum work array size
966
967 F77_XFCN (zgetri, ZGETRI, (nc, tmp_data, nr, pipvt,
968 z.fortran_vec (), lwork, info));
969
970 if (f77_exception_encountered)
971 {
972 (*current_liboctave_error_handler)
973 ("unrecoverable error in zgetri");
974 return retval;
975 }
976
977 lwork = static_cast<int> (real(z(0)));
978 lwork = (lwork < 2 *nc ? 2*nc : lwork);
979 z.resize (lwork);
980 Complex *pz = z.fortran_vec ();
981
982 info = 0;
983
984 /* Calculate the norm of the matrix, for later use */
985 double anorm;
986 if (calc_cond)
987 anorm = retval.abs().sum().row(0).max();
988
989 F77_XFCN (zgetrf, ZGETRF, (nc, nc, tmp_data, nr, pipvt, info));
962 990
963 if (f77_exception_encountered) 991 if (f77_exception_encountered)
964 (*current_liboctave_error_handler) ("unrecoverable error in zgeco"); 992 (*current_liboctave_error_handler) ("unrecoverable error in zgetrf");
965 else 993 else
966 { 994 {
967 volatile double rcond_plus_one = rcond + 1.0; 995 /* Throw-away extra info LAPACK gives so as to not change output */
968 996 rcond = 0.;
969 if (rcond_plus_one == 1.0 || xisnan (rcond)) 997 if ( info != 0)
970 info = -1; 998 info = -1;
999 else if (calc_cond)
1000 {
1001 /* Now calculate the condition number for non-singular matrix */
1002 char job = '1';
1003 Array<double> rz (2 * nc);
1004 double *prz = rz.fortran_vec ();
1005 F77_XFCN (zgecon, ZGECON, ( &job, nc, tmp_data, nr, anorm,
1006 rcond, pz, prz, info));
1007
1008 if (f77_exception_encountered)
1009 (*current_liboctave_error_handler)
1010 ("unrecoverable error in zgecon");
1011
1012 if ( info != 0)
1013 info = -1;
1014 }
971 1015
972 if (info == -1 && ! force) 1016 if (info == -1 && ! force)
973 retval = *this; // Restore contents. 1017 retval = *this; // Restore contents.
974 else 1018 else
975 { 1019 {
976 Complex *dummy = 0; 1020 F77_XFCN (zgetri, ZGETRI, (nc, tmp_data, nr, pipvt,
977 1021 pz, lwork, info));
978 F77_XFCN (zgedi, ZGEDI, (tmp_data, nr, nc, pipvt, dummy,
979 pz, 1));
980 1022
981 if (f77_exception_encountered) 1023 if (f77_exception_encountered)
982 (*current_liboctave_error_handler) 1024 (*current_liboctave_error_handler)
983 ("unrecoverable error in zgedi"); 1025 ("unrecoverable error in zgetri");
1026
1027 if ( info != 0)
1028 info = -1;
984 } 1029 }
985 } 1030 }
986 } 1031 }
987 1032
988 return retval; 1033 return retval;
989 } 1034 }
990 1035
991 ComplexMatrix 1036 ComplexMatrix
992 ComplexMatrix::pseudo_inverse (double tol) 1037 ComplexMatrix::pseudo_inverse (double tol)
1354 ComplexDET 1399 ComplexDET
1355 ComplexMatrix::determinant (void) const 1400 ComplexMatrix::determinant (void) const
1356 { 1401 {
1357 int info; 1402 int info;
1358 double rcond; 1403 double rcond;
1359 return determinant (info, rcond); 1404 return determinant (info, rcond, 0);
1360 } 1405 }
1361 1406
1362 ComplexDET 1407 ComplexDET
1363 ComplexMatrix::determinant (int& info) const 1408 ComplexMatrix::determinant (int& info) const
1364 { 1409 {
1365 double rcond; 1410 double rcond;
1366 return determinant (info, rcond); 1411 return determinant (info, rcond, 0);
1367 } 1412 }
1368 1413
1369 ComplexDET 1414 ComplexDET
1370 ComplexMatrix::determinant (int& info, double& rcond) const 1415 ComplexMatrix::determinant (int& info, double& rcond, int calc_cond) const
1371 { 1416 {
1372 ComplexDET retval; 1417 ComplexDET retval;
1373 1418
1374 int nr = rows (); 1419 int nr = rows ();
1375 int nc = cols (); 1420 int nc = cols ();
1381 d[1] = 0.0; 1426 d[1] = 0.0;
1382 retval = ComplexDET (d); 1427 retval = ComplexDET (d);
1383 } 1428 }
1384 else 1429 else
1385 { 1430 {
1386 info = 0;
1387
1388 Array<int> ipvt (nr); 1431 Array<int> ipvt (nr);
1389 int *pipvt = ipvt.fortran_vec (); 1432 int *pipvt = ipvt.fortran_vec ();
1390 1433
1391 Array<Complex> z (nr);
1392 Complex *pz = z.fortran_vec ();
1393
1394 ComplexMatrix atmp = *this; 1434 ComplexMatrix atmp = *this;
1395 Complex *tmp_data = atmp.fortran_vec (); 1435 Complex *tmp_data = atmp.fortran_vec ();
1396 1436
1397 F77_XFCN (zgeco, ZGECO, (tmp_data, nr, nr, pipvt, rcond, pz)); 1437 info = 0;
1438
1439 /* Calculate the norm of the matrix, for later use */
1440 double anorm = 0;
1441 if (calc_cond)
1442 anorm = atmp.abs().sum().row(0).max();
1443
1444 F77_XFCN (zgetrf, ZGETRF, (nr, nc, tmp_data, nr, pipvt, info));
1398 1445
1399 if (f77_exception_encountered) 1446 if (f77_exception_encountered)
1400 (*current_liboctave_error_handler) ("unrecoverable error in zgeco"); 1447 (*current_liboctave_error_handler) ("unrecoverable error in zgetrf");
1401 else 1448 else
1402 { 1449 {
1403 volatile double rcond_plus_one = rcond + 1.0; 1450 /* Throw-away extra info LAPACK gives so as to not change output */
1404 1451 rcond = 0.;
1405 if (rcond_plus_one == 1.0 || xisnan (rcond)) 1452 if ( info != 0)
1406 { 1453 {
1407 info = -1; 1454 info = -1;
1408 retval = ComplexDET (); 1455 retval = ComplexDET ();
1409 } 1456 }
1410 else 1457 else
1411 { 1458 {
1412 Complex d[2]; 1459 if (calc_cond)
1413 1460 {
1414 F77_XFCN (zgedi, ZGEDI, (tmp_data, nr, nr, pipvt, d, pz, 10)); 1461 /* Now calc the condition number for non-singular matrix */
1415 1462 char job = '1';
1416 if (f77_exception_encountered) 1463 Array<Complex> z (2*nr);
1417 (*current_liboctave_error_handler) 1464 Complex *pz = z.fortran_vec ();
1418 ("unrecoverable error in dgedi"); 1465 Array<double> rz (2*nr);
1419 else 1466 double *prz = rz.fortran_vec ();
1420 retval = ComplexDET (d); 1467
1468 F77_XFCN (zgecon, ZGECON, ( &job, nc, tmp_data, nr, anorm,
1469 rcond, pz, prz, info));
1470
1471 if (f77_exception_encountered)
1472 (*current_liboctave_error_handler)
1473 ("unrecoverable error in zgecon");
1474 }
1475
1476 if ( info != 0)
1477 {
1478 info = -1;
1479 retval = ComplexDET ();
1480 }
1481 else
1482 {
1483 Complex d[2] = { 1., 0.};
1484 for (int i=0; i<nc; i++)
1485 {
1486 if (ipvt(i) != (i+1)) d[0] = -d[0];
1487 d[0] = d[0] * atmp(i,i);
1488 if (d[0] == 0.) break;
1489 while (::abs(d[0]) < 1.)
1490 {
1491 d[0] = 10. * d[0];
1492 d[1] = d[1] - 1.;
1493 }
1494 while (::abs(d[0]) >= 10.)
1495 {
1496 d[0] = 0.1 * d[0];
1497 d[1] = d[1] + 1.;
1498 }
1499 }
1500 retval = ComplexDET (d);
1501 }
1421 } 1502 }
1422 } 1503 }
1423 } 1504 }
1424 1505
1425 return retval; 1506 return retval;
1426 } 1507 }
1427 1508
1428 ComplexMatrix 1509 ComplexMatrix
1429 ComplexMatrix::solve (const Matrix& b) const 1510 ComplexMatrix::solve (const Matrix& b) const
1492 info = 0; 1573 info = 0;
1493 1574
1494 Array<int> ipvt (nr); 1575 Array<int> ipvt (nr);
1495 int *pipvt = ipvt.fortran_vec (); 1576 int *pipvt = ipvt.fortran_vec ();
1496 1577
1497 Array<Complex> z (nr);
1498 Complex *pz = z.fortran_vec ();
1499
1500 ComplexMatrix atmp = *this; 1578 ComplexMatrix atmp = *this;
1501 Complex *tmp_data = atmp.fortran_vec (); 1579 Complex *tmp_data = atmp.fortran_vec ();
1502 1580
1503 F77_XFCN (zgeco, ZGECO, (tmp_data, nr, nr, pipvt, rcond, pz)); 1581 Array<Complex> z (2 * nc);
1582 Complex *pz = z.fortran_vec ();
1583 Array<double> rz (2 * nc);
1584 double *prz = rz.fortran_vec ();
1585
1586 /* Calculate the norm of the matrix, for later use */
1587 double anorm = atmp.abs().sum().row(0).max();
1588
1589 F77_XFCN (zgetrf, ZGETRF, (nr, nr, tmp_data, nr, pipvt, info));
1504 1590
1505 if (f77_exception_encountered) 1591 if (f77_exception_encountered)
1506 (*current_liboctave_error_handler) ("unrecoverable error in zgeco"); 1592 (*current_liboctave_error_handler) ("unrecoverable error in zgetrf");
1507 else 1593 else
1508 { 1594 {
1509 volatile double rcond_plus_one = rcond + 1.0; 1595 /* Throw-away extra info LAPACK gives so as to not change output */
1510 1596 rcond = 0.;
1511 if (rcond_plus_one == 1.0 || xisnan (rcond)) 1597 if ( info != 0)
1598 {
1599 info = -2;
1600
1601 if (sing_handler)
1602 sing_handler (rcond);
1603 else
1604 (*current_liboctave_error_handler)
1605 ("matrix singular to machine precision");
1606
1607 }
1608 else
1512 { 1609 {
1610 /* Now calculate the condition number for non-singular matrix */
1611 char job = '1';
1612 F77_XFCN (zgecon, ZGECON, ( &job, nc, tmp_data, nr, anorm,
1613 rcond, pz, prz, info));
1614
1615 if (f77_exception_encountered)
1616 (*current_liboctave_error_handler)
1617 ("unrecoverable error in zgecon");
1618
1619 if ( info != 0)
1620 info = -2;
1621
1622 volatile double rcond_plus_one = rcond + 1.0;
1623
1624 if (rcond_plus_one == 1.0 || xisnan (rcond))
1625 {
1626 info = -2;
1627
1628 if (sing_handler)
1629 sing_handler (rcond);
1630 else
1631 (*current_liboctave_error_handler)
1632 ("matrix singular to machine precision, rcond = %g",
1633 rcond);
1634 }
1635 else
1636 {
1637 retval = b;
1638 Complex *result = retval.fortran_vec ();
1639
1640 int b_nc = b.cols ();
1641
1642 char job = 'N';
1643 F77_XFCN (zgetrs, ZGETRS, (&job, nr, b_nc, tmp_data, nr,
1644 pipvt, result, b.rows(), info));
1645
1646 if (f77_exception_encountered)
1647 (*current_liboctave_error_handler)
1648 ("unrecoverable error in zgetrs");
1649 }
1650 }
1651 }
1652 }
1653
1654 return retval;
1655 }
1656
1657 ComplexColumnVector
1658 ComplexMatrix::solve (const ColumnVector& b) const
1659 {
1660 int info;
1661 double rcond;
1662 return solve (ComplexColumnVector (b), info, rcond, 0);
1663 }
1664
1665 ComplexColumnVector
1666 ComplexMatrix::solve (const ColumnVector& b, int& info) const
1667 {
1668 double rcond;
1669 return solve (ComplexColumnVector (b), info, rcond, 0);
1670 }
1671
1672 ComplexColumnVector
1673 ComplexMatrix::solve (const ColumnVector& b, int& info, double& rcond) const
1674 {
1675 return solve (ComplexColumnVector (b), info, rcond, 0);
1676 }
1677
1678 ComplexColumnVector
1679 ComplexMatrix::solve (const ColumnVector& b, int& info, double& rcond,
1680 solve_singularity_handler sing_handler) const
1681 {
1682 return solve (ComplexColumnVector (b), info, rcond, sing_handler);
1683 }
1684
1685 ComplexColumnVector
1686 ComplexMatrix::solve (const ComplexColumnVector& b) const
1687 {
1688 int info;
1689 double rcond;
1690 return solve (b, info, rcond, 0);
1691 }
1692
1693 ComplexColumnVector
1694 ComplexMatrix::solve (const ComplexColumnVector& b, int& info) const
1695 {
1696 double rcond;
1697 return solve (b, info, rcond, 0);
1698 }
1699
1700 ComplexColumnVector
1701 ComplexMatrix::solve (const ComplexColumnVector& b, int& info,
1702 double& rcond) const
1703 {
1704 return solve (b, info, rcond, 0);
1705 }
1706
1707 ComplexColumnVector
1708 ComplexMatrix::solve (const ComplexColumnVector& b, int& info,
1709 double& rcond,
1710 solve_singularity_handler sing_handler) const
1711 {
1712 ComplexColumnVector retval;
1713
1714 int nr = rows ();
1715 int nc = cols ();
1716
1717 if (nr == 0 || nc == 0 || nr != nc || nr != b.length ())
1718 (*current_liboctave_error_handler)
1719 ("matrix dimension mismatch in solution of linear equations");
1720 else
1721 {
1722 info = 0;
1723
1724 Array<int> ipvt (nr);
1725 int *pipvt = ipvt.fortran_vec ();
1726
1727 ComplexMatrix atmp = *this;
1728 Complex *tmp_data = atmp.fortran_vec ();
1729
1730 Array<Complex> z (2 * nc);
1731 Complex *pz = z.fortran_vec ();
1732 Array<double> rz (2 * nc);
1733 double *prz = rz.fortran_vec ();
1734
1735 /* Calculate the norm of the matrix, for later use */
1736 double anorm = atmp.abs().sum().row(0).max();
1737
1738 F77_XFCN (zgetrf, ZGETRF, (nr, nr, tmp_data, nr, pipvt, info));
1739
1740 if (f77_exception_encountered)
1741 (*current_liboctave_error_handler) ("unrecoverable error in zgetrf");
1742 else
1743 {
1744 /* Throw-away extra info LAPACK gives so as to not change output */
1745 rcond = 0.;
1746 if ( info != 0)
1747 {
1513 info = -2; 1748 info = -2;
1514 1749
1515 if (sing_handler) 1750 if (sing_handler)
1516 sing_handler (rcond); 1751 sing_handler (rcond);
1517 else 1752 else
1518 (*current_liboctave_error_handler) 1753 (*current_liboctave_error_handler)
1519 ("matrix singular to machine precision, rcond = %g", 1754 ("matrix singular to machine precision, rcond = %g",
1520 rcond); 1755 rcond);
1521 } 1756 }
1522 else 1757 else
1523 { 1758 {
1524 retval = b; 1759 /* Now calculate the condition number for non-singular matrix */
1525 Complex *result = retval.fortran_vec (); 1760 char job = '1';
1526 1761 F77_XFCN (zgecon, ZGECON, ( &job, nc, tmp_data, nr, anorm,
1527 int b_nc = b.cols (); 1762 rcond, pz, prz, info));
1528 1763
1529 for (volatile int j = 0; j < b_nc; j++) 1764 if (f77_exception_encountered)
1765 (*current_liboctave_error_handler)
1766 ("unrecoverable error in zgecon");
1767
1768 if ( info != 0)
1769 info = -2;
1770
1771 volatile double rcond_plus_one = rcond + 1.0;
1772
1773 if (rcond_plus_one == 1.0 || xisnan (rcond))
1530 { 1774 {
1531 F77_XFCN (zgesl, ZGESL, (tmp_data, nr, nr, pipvt, 1775 info = -2;
1532 &result[nr*j], 0)); 1776
1777 if (sing_handler)
1778 sing_handler (rcond);
1779 else
1780 (*current_liboctave_error_handler)
1781 ("matrix singular to machine precision, rcond = %g",
1782 rcond);
1783 }
1784 else
1785 {
1786 retval = b;
1787 Complex *result = retval.fortran_vec ();
1788
1789 char job = 'N';
1790 F77_XFCN (zgetrs, ZGETRS, (&job, nr, 1, tmp_data, nr, pipvt,
1791 result, b.length(), info));
1533 1792
1534 if (f77_exception_encountered) 1793 if (f77_exception_encountered)
1535 { 1794 (*current_liboctave_error_handler)
1536 (*current_liboctave_error_handler) 1795 ("unrecoverable error in zgetrs");
1537 ("unrecoverable error in dgesl"); 1796
1538
1539 break;
1540 }
1541 } 1797 }
1542 } 1798 }
1543 } 1799 }
1544 } 1800 }
1545
1546 return retval;
1547 }
1548
1549 ComplexColumnVector
1550 ComplexMatrix::solve (const ColumnVector& b) const
1551 {
1552 int info;
1553 double rcond;
1554 return solve (ComplexColumnVector (b), info, rcond, 0);
1555 }
1556
1557 ComplexColumnVector
1558 ComplexMatrix::solve (const ColumnVector& b, int& info) const
1559 {
1560 double rcond;
1561 return solve (ComplexColumnVector (b), info, rcond, 0);
1562 }
1563
1564 ComplexColumnVector
1565 ComplexMatrix::solve (const ColumnVector& b, int& info, double& rcond) const
1566 {
1567 return solve (ComplexColumnVector (b), info, rcond, 0);
1568 }
1569
1570 ComplexColumnVector
1571 ComplexMatrix::solve (const ColumnVector& b, int& info, double& rcond,
1572 solve_singularity_handler sing_handler) const
1573 {
1574 return solve (ComplexColumnVector (b), info, rcond, sing_handler);
1575 }
1576
1577 ComplexColumnVector
1578 ComplexMatrix::solve (const ComplexColumnVector& b) const
1579 {
1580 int info;
1581 double rcond;
1582 return solve (b, info, rcond, 0);
1583 }
1584
1585 ComplexColumnVector
1586 ComplexMatrix::solve (const ComplexColumnVector& b, int& info) const
1587 {
1588 double rcond;
1589 return solve (b, info, rcond, 0);
1590 }
1591
1592 ComplexColumnVector
1593 ComplexMatrix::solve (const ComplexColumnVector& b, int& info,
1594 double& rcond) const
1595 {
1596 return solve (b, info, rcond, 0);
1597 }
1598
1599 ComplexColumnVector
1600 ComplexMatrix::solve (const ComplexColumnVector& b, int& info,
1601 double& rcond,
1602 solve_singularity_handler sing_handler) const
1603 {
1604 ComplexColumnVector retval;
1605
1606 int nr = rows ();
1607 int nc = cols ();
1608
1609 if (nr == 0 || nc == 0 || nr != nc || nr != b.length ())
1610 (*current_liboctave_error_handler)
1611 ("matrix dimension mismatch in solution of linear equations");
1612 else
1613 {
1614 info = 0;
1615
1616 Array<int> ipvt (nr);
1617 int *pipvt = ipvt.fortran_vec ();
1618
1619 Array<Complex> z (nr);
1620 Complex *pz = z.fortran_vec ();
1621
1622 ComplexMatrix atmp = *this;
1623 Complex *tmp_data = atmp.fortran_vec ();
1624
1625 F77_XFCN (zgeco, ZGECO, (tmp_data, nr, nr, pipvt, rcond, pz));
1626
1627 if (f77_exception_encountered)
1628 (*current_liboctave_error_handler)
1629 ("unrecoverable error in zgeco");
1630 else
1631 {
1632 volatile double rcond_plus_one = rcond + 1.0;
1633
1634 if (rcond_plus_one == 1.0 || xisnan (rcond))
1635 {
1636 info = -2;
1637
1638 if (sing_handler)
1639 sing_handler (rcond);
1640 else
1641 (*current_liboctave_error_handler)
1642 ("matrix singular to machine precision, rcond = %g",
1643 rcond);
1644 }
1645 else
1646 {
1647 retval = b;
1648 Complex *result = retval.fortran_vec ();
1649
1650 F77_XFCN (zgesl, ZGESL, (tmp_data, nr, nr, pipvt, result, 0));
1651
1652 if (f77_exception_encountered)
1653 (*current_liboctave_error_handler)
1654 ("unrecoverable error in dgesl");
1655 }
1656 }
1657 }
1658
1659 return retval; 1801 return retval;
1660 } 1802 }
1661 1803
1662 ComplexMatrix 1804 ComplexMatrix
1663 ComplexMatrix::lssolve (const Matrix& b) const 1805 ComplexMatrix::lssolve (const Matrix& b) const
2484 2626
2485 MX_BASE_REDUCTION_OP (ComplexMatrix, ROW_EXPR, COL_EXPR, 0.0, 0.0); 2627 MX_BASE_REDUCTION_OP (ComplexMatrix, ROW_EXPR, COL_EXPR, 0.0, 0.0);
2486 2628
2487 #undef ROW_EXPR 2629 #undef ROW_EXPR
2488 #undef COL_EXPR 2630 #undef COL_EXPR
2631 }
2632
2633 Matrix ComplexMatrix::abs (void) const
2634 {
2635 int nr = rows ();
2636 int nc = cols ();
2637
2638 Matrix retval (nr, nc);
2639
2640 for (int j = 0; j < nc; j++)
2641 for (int i = 0; i < nr; i++)
2642 retval (i, j) = ::abs (elem (i, j));
2643
2644 return retval;
2489 } 2645 }
2490 2646
2491 ComplexColumnVector 2647 ComplexColumnVector
2492 ComplexMatrix::diag (void) const 2648 ComplexMatrix::diag (void) const
2493 { 2649 {
2598 2754
2599 Complex tmp_min = elem (i, idx_j); 2755 Complex tmp_min = elem (i, idx_j);
2600 2756
2601 bool real_only = row_is_real_only (i); 2757 bool real_only = row_is_real_only (i);
2602 2758
2603 double abs_min = real_only ? real (tmp_min) : abs (tmp_min); 2759 double abs_min = real_only ? real (tmp_min) : ::abs (tmp_min);
2604 2760
2605 if (xisnan (tmp_min)) 2761 if (xisnan (tmp_min))
2606 idx_j = -1; 2762 idx_j = -1;
2607 else 2763 else
2608 { 2764 {
2609 for (int j = 1; j < nc; j++) 2765 for (int j = 1; j < nc; j++)
2610 { 2766 {
2611 Complex tmp = elem (i, j); 2767 Complex tmp = elem (i, j);
2612 2768
2613 double abs_tmp = real_only ? real (tmp) : abs (tmp); 2769 double abs_tmp = real_only ? real (tmp) : ::abs (tmp);
2614 2770
2615 if (xisnan (tmp)) 2771 if (xisnan (tmp))
2616 { 2772 {
2617 idx_j = -1; 2773 idx_j = -1;
2618 break; 2774 break;
2660 2816
2661 Complex tmp_max = elem (i, idx_j); 2817 Complex tmp_max = elem (i, idx_j);
2662 2818
2663 bool real_only = row_is_real_only (i); 2819 bool real_only = row_is_real_only (i);
2664 2820
2665 double abs_max = real_only ? real (tmp_max) : abs (tmp_max); 2821 double abs_max = real_only ? real (tmp_max) : ::abs (tmp_max);
2666 2822
2667 if (xisnan (tmp_max)) 2823 if (xisnan (tmp_max))
2668 idx_j = -1; 2824 idx_j = -1;
2669 else 2825 else
2670 { 2826 {
2671 for (int j = 1; j < nc; j++) 2827 for (int j = 1; j < nc; j++)
2672 { 2828 {
2673 Complex tmp = elem (i, j); 2829 Complex tmp = elem (i, j);
2674 2830
2675 double abs_tmp = real_only ? real (tmp) : abs (tmp); 2831 double abs_tmp = real_only ? real (tmp) : ::abs (tmp);
2676 2832
2677 if (xisnan (tmp)) 2833 if (xisnan (tmp))
2678 { 2834 {
2679 idx_j = -1; 2835 idx_j = -1;
2680 break; 2836 break;
2722 2878
2723 Complex tmp_min = elem (idx_i, j); 2879 Complex tmp_min = elem (idx_i, j);
2724 2880
2725 bool real_only = column_is_real_only (j); 2881 bool real_only = column_is_real_only (j);
2726 2882
2727 double abs_min = real_only ? real (tmp_min) : abs (tmp_min); 2883 double abs_min = real_only ? real (tmp_min) : ::abs (tmp_min);
2728 2884
2729 if (xisnan (tmp_min)) 2885 if (xisnan (tmp_min))
2730 idx_i = -1; 2886 idx_i = -1;
2731 else 2887 else
2732 { 2888 {
2733 for (int i = 1; i < nr; i++) 2889 for (int i = 1; i < nr; i++)
2734 { 2890 {
2735 Complex tmp = elem (i, j); 2891 Complex tmp = elem (i, j);
2736 2892
2737 double abs_tmp = real_only ? real (tmp) : abs (tmp); 2893 double abs_tmp = real_only ? real (tmp) : ::abs (tmp);
2738 2894
2739 if (xisnan (tmp)) 2895 if (xisnan (tmp))
2740 { 2896 {
2741 idx_i = -1; 2897 idx_i = -1;
2742 break; 2898 break;
2784 2940
2785 Complex tmp_max = elem (idx_i, j); 2941 Complex tmp_max = elem (idx_i, j);
2786 2942
2787 bool real_only = column_is_real_only (j); 2943 bool real_only = column_is_real_only (j);
2788 2944
2789 double abs_max = real_only ? real (tmp_max) : abs (tmp_max); 2945 double abs_max = real_only ? real (tmp_max) : ::abs (tmp_max);
2790 2946
2791 if (xisnan (tmp_max)) 2947 if (xisnan (tmp_max))
2792 idx_i = -1; 2948 idx_i = -1;
2793 else 2949 else
2794 { 2950 {
2795 for (int i = 1; i < nr; i++) 2951 for (int i = 1; i < nr; i++)
2796 { 2952 {
2797 Complex tmp = elem (i, j); 2953 Complex tmp = elem (i, j);
2798 2954
2799 double abs_tmp = real_only ? real (tmp) : abs (tmp); 2955 double abs_tmp = real_only ? real (tmp) : ::abs (tmp);
2800 2956
2801 if (xisnan (tmp)) 2957 if (xisnan (tmp))
2802 { 2958 {
2803 idx_i = -1; 2959 idx_i = -1;
2804 break; 2960 break;