Mercurial > hg > octave-lyh

comparison liboctave/dMatrix.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 c038c2947ee1
comparison
equal deleted inserted replaced
5784:70f7659d0fb9 5785:6b9cec830d72
105 F77_RET_T 105 F77_RET_T
106 F77_FUNC (dgelss, DGELSS) (const octave_idx_type&, const octave_idx_type&, const octave_idx_type&, 106 F77_FUNC (dgelss, DGELSS) (const octave_idx_type&, const octave_idx_type&, const octave_idx_type&,
107 double*, const octave_idx_type&, double*, 107 double*, const octave_idx_type&, double*,
108 const octave_idx_type&, double*, double&, octave_idx_type&, 108 const octave_idx_type&, double*, double&, octave_idx_type&,
109 double*, const octave_idx_type&, octave_idx_type&); 109 double*, const octave_idx_type&, octave_idx_type&);
110
111 F77_RET_T
112 F77_FUNC (dpotrf, DPOTRF) (F77_CONST_CHAR_ARG_DECL, const octave_idx_type&,
113 double *, const octave_idx_type&,
114 octave_idx_type& F77_CHAR_ARG_LEN_DECL);
115
116 F77_RET_T
117 F77_FUNC (dpocon, DPOCON) (F77_CONST_CHAR_ARG_DECL, const octave_idx_type&,
118 double*, const octave_idx_type&, const double&,
119 double&, double*, octave_idx_type*,
120 octave_idx_type& F77_CHAR_ARG_LEN_DECL);
121 F77_RET_T
122 F77_FUNC (dpotrs, DPOTRS) (F77_CONST_CHAR_ARG_DECL, const octave_idx_type&,
123 const octave_idx_type&, const double*,
124 const octave_idx_type&, double*,
125 const octave_idx_type&, octave_idx_type&
126 F77_CHAR_ARG_LEN_DECL);
127
128 F77_RET_T
129 F77_FUNC (dtrcon, DTRCON) (F77_CONST_CHAR_ARG_DECL, F77_CONST_CHAR_ARG_DECL,
130 F77_CONST_CHAR_ARG_DECL, const octave_idx_type&,
131 const double*, const octave_idx_type&, double&,
132 double*, octave_idx_type*, octave_idx_type&
133 F77_CHAR_ARG_LEN_DECL
134 F77_CHAR_ARG_LEN_DECL
135 F77_CHAR_ARG_LEN_DECL);
136 F77_RET_T
137 F77_FUNC (dtrtrs, DTRTRS) (F77_CONST_CHAR_ARG_DECL, F77_CONST_CHAR_ARG_DECL,
138 F77_CONST_CHAR_ARG_DECL, const octave_idx_type&,
139 const octave_idx_type&, const double*,
140 const octave_idx_type&, double*,
141 const octave_idx_type&, octave_idx_type&
142 F77_CHAR_ARG_LEN_DECL
143 F77_CHAR_ARG_LEN_DECL
144 F77_CHAR_ARG_LEN_DECL);
110 145
111 // Note that the original complex fft routines were not written for 146 // Note that the original complex fft routines were not written for
112 // double complex arguments. They have been modified by adding an 147 // double complex arguments. They have been modified by adding an
113 // implicit double precision (a-h,o-z) statement at the beginning of 148 // implicit double precision (a-h,o-z) statement at the beginning of
114 // each subroutine. 149 // each subroutine.
1152 1187
1153 return retval; 1188 return retval;
1154 } 1189 }
1155 1190
1156 Matrix 1191 Matrix
1157 Matrix::solve (const Matrix& b) const 1192 Matrix::utsolve (MatrixType &mattype, const Matrix& b, octave_idx_type& info,
1158 { 1193 double& rcond, solve_singularity_handler sing_handler,
1159 octave_idx_type info; 1194 bool calc_cond) const
1160 double rcond; 1195 {
1161 return solve (b, info, rcond, 0); 1196 Matrix retval;
1162 } 1197
1163 1198 octave_idx_type nr = rows ();
1164 Matrix 1199 octave_idx_type nc = cols ();
1165 Matrix::solve (const Matrix& b, octave_idx_type& info) const 1200
1166 { 1201 if (nr == 0 || nc == 0 || nr != b.rows ())
1167 double rcond; 1202 (*current_liboctave_error_handler)
1168 return solve (b, info, rcond, 0); 1203 ("matrix dimension mismatch solution of linear equations");
1169 } 1204 else
1170 1205 {
1171 Matrix 1206 volatile int typ = mattype.type ();
1172 Matrix::solve (const Matrix& b, octave_idx_type& info, double& rcond) const 1207
1173 { 1208 if (typ == MatrixType::Permuted_Upper ||
1174 return solve (b, info, rcond, 0); 1209 typ == MatrixType::Upper)
1175 } 1210 {
1176 1211 octave_idx_type b_nc = b.cols ();
1177 Matrix 1212 rcond = 1.;
1178 Matrix::solve (const Matrix& b, octave_idx_type& info, double& rcond, 1213 info = 0;
1179 solve_singularity_handler sing_handler) const 1214
1215 if (typ == MatrixType::Permuted_Upper)
1216 {
1217 (*current_liboctave_error_handler)
1218 ("Permuted triangular matrix not implemented");
1219 }
1220 else
1221 {
1222 const double *tmp_data = fortran_vec ();
1223
1224 if (calc_cond)
1225 {
1226 char norm = '1';
1227 char uplo = 'U';
1228 char dia = 'N';
1229
1230 Array<double> z (3 * nc);
1231 double *pz = z.fortran_vec ();
1232 Array<octave_idx_type> iz (nc);
1233 octave_idx_type *piz = iz.fortran_vec ();
1234
1235 F77_XFCN (dtrcon, DTRCON, (F77_CONST_CHAR_ARG2 (&norm, 1),
1236 F77_CONST_CHAR_ARG2 (&uplo, 1),
1237 F77_CONST_CHAR_ARG2 (&dia, 1),
1238 nr, tmp_data, nr, rcond,
1239 pz, piz, info
1240 F77_CHAR_ARG_LEN (1)
1241 F77_CHAR_ARG_LEN (1)
1242 F77_CHAR_ARG_LEN (1)));
1243
1244 if (f77_exception_encountered)
1245 (*current_liboctave_error_handler)
1246 ("unrecoverable error in dtrcon");
1247
1248 if (info != 0)
1249 info = -2;
1250
1251 volatile double rcond_plus_one = rcond + 1.0;
1252
1253 if (rcond_plus_one == 1.0 || xisnan (rcond))
1254 {
1255 info = -2;
1256
1257 if (sing_handler)
1258 sing_handler (rcond);
1259 else
1260 (*current_liboctave_error_handler)
1261 ("matrix singular to machine precision, rcond = %g",
1262 rcond);
1263 }
1264 }
1265
1266 if (info == 0)
1267 {
1268 retval = b;
1269 double *result = retval.fortran_vec ();
1270
1271 char uplo = 'U';
1272 char trans = 'N';
1273 char dia = 'N';
1274
1275 F77_XFCN (dtrtrs, DTRTRS, (F77_CONST_CHAR_ARG2 (&uplo, 1),
1276 F77_CONST_CHAR_ARG2 (&trans, 1),
1277 F77_CONST_CHAR_ARG2 (&dia, 1),
1278 nr, b_nc, tmp_data, nr,
1279 result, nr, info
1280 F77_CHAR_ARG_LEN (1)
1281 F77_CHAR_ARG_LEN (1)
1282 F77_CHAR_ARG_LEN (1)));
1283
1284 if (f77_exception_encountered)
1285 (*current_liboctave_error_handler)
1286 ("unrecoverable error in dtrtrs");
1287 }
1288 }
1289 }
1290 else
1291 (*current_liboctave_error_handler) ("incorrect matrix type");
1292 }
1293
1294 return retval;
1295 }
1296
1297 Matrix
1298 Matrix::ltsolve (MatrixType &mattype, const Matrix& b, octave_idx_type& info,
1299 double& rcond, solve_singularity_handler sing_handler,
1300 bool calc_cond) const
1301 {
1302 Matrix retval;
1303
1304 octave_idx_type nr = rows ();
1305 octave_idx_type nc = cols ();
1306
1307 if (nr == 0 || nc == 0 || nr != b.rows ())
1308 (*current_liboctave_error_handler)
1309 ("matrix dimension mismatch solution of linear equations");
1310 else
1311 {
1312 volatile int typ = mattype.type ();
1313
1314 if (typ == MatrixType::Permuted_Lower ||
1315 typ == MatrixType::Lower)
1316 {
1317 octave_idx_type b_nc = b.cols ();
1318 rcond = 1.;
1319 info = 0;
1320
1321 if (typ == MatrixType::Permuted_Lower)
1322 {
1323 (*current_liboctave_error_handler)
1324 ("Permuted triangular matrix not implemented");
1325 }
1326 else
1327 {
1328 const double *tmp_data = fortran_vec ();
1329
1330 if (calc_cond)
1331 {
1332 char norm = '1';
1333 char uplo = 'L';
1334 char dia = 'N';
1335
1336 Array<double> z (3 * nc);
1337 double *pz = z.fortran_vec ();
1338 Array<octave_idx_type> iz (nc);
1339 octave_idx_type *piz = iz.fortran_vec ();
1340
1341 F77_XFCN (dtrcon, DTRCON, (F77_CONST_CHAR_ARG2 (&norm, 1),
1342 F77_CONST_CHAR_ARG2 (&uplo, 1),
1343 F77_CONST_CHAR_ARG2 (&dia, 1),
1344 nr, tmp_data, nr, rcond,
1345 pz, piz, info
1346 F77_CHAR_ARG_LEN (1)
1347 F77_CHAR_ARG_LEN (1)
1348 F77_CHAR_ARG_LEN (1)));
1349
1350 if (f77_exception_encountered)
1351 (*current_liboctave_error_handler)
1352 ("unrecoverable error in dtrcon");
1353
1354 if (info != 0)
1355 info = -2;
1356
1357 volatile double rcond_plus_one = rcond + 1.0;
1358
1359 if (rcond_plus_one == 1.0 || xisnan (rcond))
1360 {
1361 info = -2;
1362
1363 if (sing_handler)
1364 sing_handler (rcond);
1365 else
1366 (*current_liboctave_error_handler)
1367 ("matrix singular to machine precision, rcond = %g",
1368 rcond);
1369 }
1370 }
1371
1372 if (info == 0)
1373 {
1374 retval = b;
1375 double *result = retval.fortran_vec ();
1376
1377 char uplo = 'L';
1378 char trans = 'N';
1379 char dia = 'N';
1380
1381 F77_XFCN (dtrtrs, DTRTRS, (F77_CONST_CHAR_ARG2 (&uplo, 1),
1382 F77_CONST_CHAR_ARG2 (&trans, 1),
1383 F77_CONST_CHAR_ARG2 (&dia, 1),
1384 nr, b_nc, tmp_data, nr,
1385 result, nr, info
1386 F77_CHAR_ARG_LEN (1)
1387 F77_CHAR_ARG_LEN (1)
1388 F77_CHAR_ARG_LEN (1)));
1389
1390 if (f77_exception_encountered)
1391 (*current_liboctave_error_handler)
1392 ("unrecoverable error in dtrtrs");
1393 }
1394 }
1395 }
1396 else
1397 (*current_liboctave_error_handler) ("incorrect matrix type");
1398 }
1399
1400 return retval;
1401 }
1402
1403 Matrix
1404 Matrix::fsolve (MatrixType &mattype, const Matrix& b, octave_idx_type& info,
1405 double& rcond, solve_singularity_handler sing_handler,
1406 bool calc_cond) const
1180 { 1407 {
1181 Matrix retval; 1408 Matrix retval;
1182 1409
1183 octave_idx_type nr = rows (); 1410 octave_idx_type nr = rows ();
1184 octave_idx_type nc = cols (); 1411 octave_idx_type nc = cols ();
1186 if (nr == 0 || nc == 0 || nr != nc || nr != b.rows ()) 1413 if (nr == 0 || nc == 0 || nr != nc || nr != b.rows ())
1187 (*current_liboctave_error_handler) 1414 (*current_liboctave_error_handler)
1188 ("matrix dimension mismatch solution of linear equations"); 1415 ("matrix dimension mismatch solution of linear equations");
1189 else 1416 else
1190 { 1417 {
1191 info = 0; 1418 volatile int typ = mattype.type ();
1192 1419
1193 Array<octave_idx_type> ipvt (nr); 1420 // Calculate the norm of the matrix, for later use.
1194 octave_idx_type *pipvt = ipvt.fortran_vec (); 1421 double anorm = -1.;
1195 1422
1196 Matrix atmp = *this; 1423 if (typ == MatrixType::Hermitian)
1197 double *tmp_data = atmp.fortran_vec ();
1198
1199 Array<double> z (4 * nc);
1200 double *pz = z.fortran_vec ();
1201 Array<octave_idx_type> iz (nc);
1202 octave_idx_type *piz = iz.fortran_vec ();
1203
1204 // Calculate the norm of the matrix, for later use.
1205 double anorm = atmp.abs().sum().row(static_cast<octave_idx_type>(0)).max();
1206
1207 F77_XFCN (dgetrf, DGETRF, (nr, nr, tmp_data, nr, pipvt, info));
1208
1209 if (f77_exception_encountered)
1210 (*current_liboctave_error_handler) ("unrecoverable error in dgetrf");
1211 else
1212 { 1424 {
1213 // Throw-away extra info LAPACK gives so as to not change output. 1425 info = 0;
1214 rcond = 0.0; 1426 char job = 'L';
1215 if (info != 0) 1427 Matrix atmp = *this;
1428 double *tmp_data = atmp.fortran_vec ();
1429 anorm = atmp.abs().sum().row(static_cast<octave_idx_type>(0)).max();
1430
1431 F77_XFCN (dpotrf, DPOTRF, (F77_CONST_CHAR_ARG2 (&job, 1), nr,
1432 tmp_data, nr, info
1433 F77_CHAR_ARG_LEN (1)));
1434
1435 if (f77_exception_encountered)
1436 (*current_liboctave_error_handler)
1437 ("unrecoverable error in dpotrf");
1438 else
1216 { 1439 {
1217 info = -2; 1440 // Throw-away extra info LAPACK gives so as to not change output.
1218 1441 rcond = 0.0;
1219 if (sing_handler) 1442 if (info != 0)
1220 sing_handler (rcond); 1443 {
1221 else 1444 info = -2;
1222 (*current_liboctave_error_handler) 1445
1223 ("matrix singular to machine precision"); 1446 mattype.mark_as_unsymmetric ();
1224 1447 typ = MatrixType::Full;
1225 } 1448 }
1226 else 1449 else
1450 {
1451 if (calc_cond)
1452 {
1453 Array<double> z (3 * nc);
1454 double *pz = z.fortran_vec ();
1455 Array<octave_idx_type> iz (nc);
1456 octave_idx_type *piz = iz.fortran_vec ();
1457
1458 F77_XFCN (dpocon, DPOCON, (F77_CONST_CHAR_ARG2 (&job, 1),
1459 nr, tmp_data, nr, anorm,
1460 rcond, pz, piz, info
1461 F77_CHAR_ARG_LEN (1)));
1462
1463 if (f77_exception_encountered)
1464 (*current_liboctave_error_handler)
1465 ("unrecoverable error in dpocon");
1466
1467 if (info != 0)
1468 info = -2;
1469
1470 volatile double rcond_plus_one = rcond + 1.0;
1471
1472 if (rcond_plus_one == 1.0 || xisnan (rcond))
1473 {
1474 info = -2;
1475
1476 if (sing_handler)
1477 sing_handler (rcond);
1478 else
1479 (*current_liboctave_error_handler)
1480 ("matrix singular to machine precision, rcond = %g",
1481 rcond);
1482 }
1483 }
1484
1485 if (info == 0)
1486 {
1487 retval = b;
1488 double *result = retval.fortran_vec ();
1489
1490 octave_idx_type b_nc = b.cols ();
1491
1492 F77_XFCN (dpotrs, DPOTRS, (F77_CONST_CHAR_ARG2 (&job, 1),
1493 nr, b_nc, tmp_data, nr,
1494 result, b.rows(), info
1495 F77_CHAR_ARG_LEN (1)));
1496
1497 if (f77_exception_encountered)
1498 (*current_liboctave_error_handler)
1499 ("unrecoverable error in dpotrs");
1500 }
1501 else
1502 {
1503 mattype.mark_as_unsymmetric ();
1504 typ = MatrixType::Full;
1505 }
1506 }
1507 }
1508 }
1509
1510 if (typ == MatrixType::Full)
1511 {
1512 info = 0;
1513
1514 Array<octave_idx_type> ipvt (nr);
1515 octave_idx_type *pipvt = ipvt.fortran_vec ();
1516
1517 Matrix atmp = *this;
1518 double *tmp_data = atmp.fortran_vec ();
1519 if(anorm < 0.)
1520 anorm = atmp.abs().sum().row(static_cast<octave_idx_type>(0)).max();
1521
1522 Array<double> z (4 * nc);
1523 double *pz = z.fortran_vec ();
1524 Array<octave_idx_type> iz (nc);
1525 octave_idx_type *piz = iz.fortran_vec ();
1526
1527 F77_XFCN (dgetrf, DGETRF, (nr, nr, tmp_data, nr, pipvt, info));
1528
1529 if (f77_exception_encountered)
1530 (*current_liboctave_error_handler)
1531 ("unrecoverable error in dgetrf");
1532 else
1227 { 1533 {
1228 // Now calculate the condition number for non-singular matrix. 1534 // Throw-away extra info LAPACK gives so as to not change output.
1229 char job = '1'; 1535 rcond = 0.0;
1230 F77_XFCN (dgecon, DGECON, (F77_CONST_CHAR_ARG2 (&job, 1),
1231 nc, tmp_data, nr, anorm,
1232 rcond, pz, piz, info
1233 F77_CHAR_ARG_LEN (1)));
1234
1235 if (f77_exception_encountered)
1236 (*current_liboctave_error_handler)
1237 ("unrecoverable error in dgecon");
1238
1239 if (info != 0) 1536 if (info != 0)
1240 info = -2;
1241
1242 volatile double rcond_plus_one = rcond + 1.0;
1243
1244 if (rcond_plus_one == 1.0 || xisnan (rcond))
1245 { 1537 {
1246 info = -2; 1538 info = -2;
1247 1539
1248 if (sing_handler) 1540 if (sing_handler)
1249 sing_handler (rcond); 1541 sing_handler (rcond);
1250 else 1542 else
1251 (*current_liboctave_error_handler) 1543 (*current_liboctave_error_handler)
1252 ("matrix singular to machine precision, rcond = %g", 1544 ("matrix singular to machine precision");
1253 rcond); 1545
1546 mattype.mark_as_rectangular ();
1254 } 1547 }
1255 else 1548 else
1256 { 1549 {
1257 retval = b; 1550 if (calc_cond)
1258 double *result = retval.fortran_vec (); 1551 {
1259 1552 // Now calculate the condition number for
1260 octave_idx_type b_nc = b.cols (); 1553 // non-singular matrix.
1261 1554 char job = '1';
1262 job = 'N'; 1555 F77_XFCN (dgecon, DGECON, (F77_CONST_CHAR_ARG2 (&job, 1),
1263 F77_XFCN (dgetrs, DGETRS, (F77_CONST_CHAR_ARG2 (&job, 1), 1556 nc, tmp_data, nr, anorm,
1264 nr, b_nc, tmp_data, nr, 1557 rcond, pz, piz, info
1265 pipvt, result, b.rows(), info 1558 F77_CHAR_ARG_LEN (1)));
1266 F77_CHAR_ARG_LEN (1))); 1559
1560 if (f77_exception_encountered)
1561 (*current_liboctave_error_handler)
1562 ("unrecoverable error in dgecon");
1563
1564 if (info != 0)
1565 info = -2;
1566
1567 volatile double rcond_plus_one = rcond + 1.0;
1568
1569 if (rcond_plus_one == 1.0 || xisnan (rcond))
1570 {
1571 info = -2;
1572
1573 if (sing_handler)
1574 sing_handler (rcond);
1575 else
1576 (*current_liboctave_error_handler)
1577 ("matrix singular to machine precision, rcond = %g",
1578 rcond);
1579 }
1580 }
1581
1582 if (info == 0)
1583 {
1584 retval = b;
1585 double *result = retval.fortran_vec ();
1586
1587 octave_idx_type b_nc = b.cols ();
1588
1589 char job = 'N';
1590 F77_XFCN (dgetrs, DGETRS, (F77_CONST_CHAR_ARG2 (&job, 1),
1591 nr, b_nc, tmp_data, nr,
1592 pipvt, result, b.rows(), info
1593 F77_CHAR_ARG_LEN (1)));
1267 1594
1268 if (f77_exception_encountered) 1595 if (f77_exception_encountered)
1269 (*current_liboctave_error_handler) 1596 (*current_liboctave_error_handler)
1270 ("unrecoverable error in dgetrs"); 1597 ("unrecoverable error in dgetrs");
1598 }
1599 else
1600 mattype.mark_as_rectangular ();
1271 } 1601 }
1272 } 1602 }
1273 } 1603 }
1274 } 1604 else if (typ != MatrixType::Hermitian)
1275 1605 (*current_liboctave_error_handler) ("incorrect matrix type");
1276 return retval; 1606 }
1607
1608 return retval;
1609 }
1610
1611 Matrix
1612 Matrix::solve (MatrixType &typ, const Matrix& b) const
1613 {
1614 octave_idx_type info;
1615 double rcond;
1616 return solve (typ, b, info, rcond, 0);
1617 }
1618
1619 Matrix
1620 Matrix::solve (MatrixType &typ, const Matrix& b, octave_idx_type& info,
1621 double& rcond) const
1622 {
1623 return solve (typ, b, info, rcond, 0);
1624 }
1625
1626 Matrix
1627 Matrix::solve (MatrixType &mattype, const Matrix& b, octave_idx_type& info,
1628 double& rcond, solve_singularity_handler sing_handler,
1629 bool singular_fallback) const
1630 {
1631 Matrix retval;
1632 int typ = mattype.type ();
1633
1634 if (typ == MatrixType::Unknown)
1635 typ = mattype.type (*this);
1636
1637 // Only calculate the condition number for LU/Cholesky
1638 if (typ == MatrixType::Upper || typ == MatrixType::Permuted_Upper)
1639 retval = utsolve (mattype, b, info, rcond, sing_handler, false);
1640 else if (typ == MatrixType::Lower || typ == MatrixType::Permuted_Lower)
1641 retval = ltsolve (mattype, b, info, rcond, sing_handler, false);
1642 else if (typ == MatrixType::Full || typ == MatrixType::Hermitian)
1643 retval = fsolve (mattype, b, info, rcond, sing_handler, true);
1644 else if (typ != MatrixType::Rectangular)
1645 {
1646 (*current_liboctave_error_handler) ("unknown matrix type");
1647 return Matrix ();
1648 }
1649
1650 // Rectangular or one of the above solvers flags a singular matrix
1651 if (singular_fallback && mattype.type () == MatrixType::Rectangular)
1652 {
1653 octave_idx_type rank;
1654 retval = lssolve (b, info, rank);
1655 }
1656
1657 return retval;
1658 }
1659
1660 ComplexMatrix
1661 Matrix::solve (MatrixType &typ, const ComplexMatrix& b) const
1662 {
1663 ComplexMatrix tmp (*this);
1664 return tmp.solve (typ, b);
1665 }
1666
1667 ComplexMatrix
1668 Matrix::solve (MatrixType &typ, const ComplexMatrix& b,
1669 octave_idx_type& info) const
1670 {
1671 ComplexMatrix tmp (*this);
1672 return tmp.solve (typ, b, info);
1673 }
1674
1675 ComplexMatrix
1676 Matrix::solve (MatrixType &typ, const ComplexMatrix& b, octave_idx_type& info,
1677 double& rcond) const
1678 {
1679 ComplexMatrix tmp (*this);
1680 return tmp.solve (typ, b, info, rcond);
1681 }
1682
1683 ComplexMatrix
1684 Matrix::solve (MatrixType &typ, const ComplexMatrix& b, octave_idx_type& info,
1685 double& rcond, solve_singularity_handler sing_handler,
1686 bool singular_fallback) const
1687 {
1688 ComplexMatrix tmp (*this);
1689 return tmp.solve (typ, b, info, rcond, sing_handler, singular_fallback);
1690 }
1691
1692 ColumnVector
1693 Matrix::solve (MatrixType &typ, const ColumnVector& b) const
1694 {
1695 octave_idx_type info; double rcond;
1696 return solve (typ, b, info, rcond);
1697 }
1698
1699 ColumnVector
1700 Matrix::solve (MatrixType &typ, const ColumnVector& b,
1701 octave_idx_type& info) const
1702 {
1703 double rcond;
1704 return solve (typ, b, info, rcond);
1705 }
1706
1707 ColumnVector
1708 Matrix::solve (MatrixType &typ, const ColumnVector& b, octave_idx_type& info,
1709 double& rcond) const
1710 {
1711 return solve (typ, b, info, rcond, 0);
1712 }
1713
1714 ColumnVector
1715 Matrix::solve (MatrixType &typ, const ColumnVector& b, octave_idx_type& info,
1716 double& rcond, solve_singularity_handler sing_handler) const
1717 {
1718 Matrix tmp (b);
1719 return solve (typ, tmp, info, rcond, sing_handler).column(static_cast<octave_idx_type> (0));
1720 }
1721
1722 ComplexColumnVector
1723 Matrix::solve (MatrixType &typ, const ComplexColumnVector& b) const
1724 {
1725 ComplexMatrix tmp (*this);
1726 return tmp.solve (typ, b);
1727 }
1728
1729 ComplexColumnVector
1730 Matrix::solve (MatrixType &typ, const ComplexColumnVector& b,
1731 octave_idx_type& info) const
1732 {
1733 ComplexMatrix tmp (*this);
1734 return tmp.solve (typ, b, info);
1735 }
1736
1737 ComplexColumnVector
1738 Matrix::solve (MatrixType &typ, const ComplexColumnVector& b,
1739 octave_idx_type& info, double& rcond) const
1740 {
1741 ComplexMatrix tmp (*this);
1742 return tmp.solve (typ, b, info, rcond);
1743 }
1744
1745 ComplexColumnVector
1746 Matrix::solve (MatrixType &typ, const ComplexColumnVector& b,
1747 octave_idx_type& info, double& rcond,
1748 solve_singularity_handler sing_handler) const
1749 {
1750 ComplexMatrix tmp (*this);
1751 return tmp.solve(typ, b, info, rcond, sing_handler);
1752 }
1753
1754 Matrix
1755 Matrix::solve (const Matrix& b) const
1756 {
1757 octave_idx_type info;
1758 double rcond;
1759 return solve (b, info, rcond, 0);
1760 }
1761
1762 Matrix
1763 Matrix::solve (const Matrix& b, octave_idx_type& info) const
1764 {
1765 double rcond;
1766 return solve (b, info, rcond, 0);
1767 }
1768
1769 Matrix
1770 Matrix::solve (const Matrix& b, octave_idx_type& info, double& rcond) const
1771 {
1772 return solve (b, info, rcond, 0);
1773 }
1774
1775 Matrix
1776 Matrix::solve (const Matrix& b, octave_idx_type& info,
1777 double& rcond, solve_singularity_handler sing_handler) const
1778 {
1779 MatrixType mattype (*this);
1780 return solve (mattype, b, info, rcond, sing_handler);
1277 } 1781 }
1278 1782
1279 ComplexMatrix 1783 ComplexMatrix
1280 Matrix::solve (const ComplexMatrix& b) const 1784 Matrix::solve (const ComplexMatrix& b) const
1281 { 1785 {
1327 1831
1328 ColumnVector 1832 ColumnVector
1329 Matrix::solve (const ColumnVector& b, octave_idx_type& info, double& rcond, 1833 Matrix::solve (const ColumnVector& b, octave_idx_type& info, double& rcond,
1330 solve_singularity_handler sing_handler) const 1834 solve_singularity_handler sing_handler) const
1331 { 1835 {
1332 ColumnVector retval; 1836 MatrixType mattype (*this);
1333 1837 return solve (mattype, b, info, rcond, sing_handler);
1334 octave_idx_type nr = rows ();
1335 octave_idx_type nc = cols ();
1336
1337 if (nr == 0 || nc == 0 || nr != nc || nr != b.length ())
1338 (*current_liboctave_error_handler)
1339 ("matrix dimension mismatch solution of linear equations");
1340 else
1341 {
1342 info = 0;
1343
1344 Array<octave_idx_type> ipvt (nr);
1345 octave_idx_type *pipvt = ipvt.fortran_vec ();
1346
1347 Matrix atmp = *this;
1348 double *tmp_data = atmp.fortran_vec ();
1349
1350 Array<double> z (4 * nc);
1351 double *pz = z.fortran_vec ();
1352 Array<octave_idx_type> iz (nc);
1353 octave_idx_type *piz = iz.fortran_vec ();
1354
1355 // Calculate the norm of the matrix, for later use.
1356 double anorm = atmp.abs().sum().row(static_cast<octave_idx_type>(0)).max();
1357
1358 F77_XFCN (dgetrf, DGETRF, (nr, nr, tmp_data, nr, pipvt, info));
1359
1360 if (f77_exception_encountered)
1361 (*current_liboctave_error_handler) ("unrecoverable error in dgetrf");
1362 else
1363 {
1364 // Throw-away extra info LAPACK gives so as to not change output.
1365 rcond = 0.0;
1366 if (info > 0)
1367 {
1368 info = -2;
1369
1370 if (sing_handler)
1371 sing_handler (rcond);
1372 else
1373 (*current_liboctave_error_handler)
1374 ("matrix singular to machine precision");
1375
1376 }
1377 else
1378 {
1379 // Now calculate the condition number for non-singular matrix.
1380 char job = '1';
1381 F77_XFCN (dgecon, DGECON, (F77_CONST_CHAR_ARG2 (&job, 1),
1382 nc, tmp_data, nr, anorm,
1383 rcond, pz, piz, info
1384 F77_CHAR_ARG_LEN (1)));
1385
1386 if (f77_exception_encountered)
1387 (*current_liboctave_error_handler)
1388 ("unrecoverable error in dgecon");
1389
1390 if (info != 0)
1391 info = -2;
1392
1393 volatile double rcond_plus_one = rcond + 1.0;
1394
1395 if (rcond_plus_one == 1.0 || xisnan (rcond))
1396 {
1397 info = -2;
1398
1399 if (sing_handler)
1400 sing_handler (rcond);
1401 else
1402 (*current_liboctave_error_handler)
1403 ("matrix singular to machine precision, rcond = %g",
1404 rcond);
1405 }
1406 else
1407 {
1408 retval = b;
1409 double *result = retval.fortran_vec ();
1410
1411 job = 'N';
1412 F77_XFCN (dgetrs, DGETRS, (F77_CONST_CHAR_ARG2 (&job, 1),
1413 nr, 1, tmp_data, nr, pipvt,
1414 result, b.length(), info
1415 F77_CHAR_ARG_LEN (1)));
1416
1417 if (f77_exception_encountered)
1418 (*current_liboctave_error_handler)
1419 ("unrecoverable error in dgetrs");
1420 }
1421 }
1422 }
1423 }
1424
1425 return retval;
1426 } 1838 }
1427 1839
1428 ComplexColumnVector 1840 ComplexColumnVector
1429 Matrix::solve (const ComplexColumnVector& b) const 1841 Matrix::solve (const ComplexColumnVector& b) const
1430 { 1842 {