Mercurial > octave-nkf

comparison liboctave/CSparse.cc @ 5322:22994a5730f9

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
[project @ 2005-04-29 13:04:24 by dbateman]
author dbateman
date Fri, 29 Apr 2005 13:04:25 +0000
parents 4c8a2e4e0717
children a103c41e68b2
comparison
equal deleted inserted replaced
5321:84b72a402b86 5322:22994a5730f9
96 const Complex*, const Complex*, const octave_idx_type*, 96 const Complex*, const Complex*, const octave_idx_type*,
97 Complex *, const octave_idx_type&, octave_idx_type& 97 Complex *, const octave_idx_type&, octave_idx_type&
98 F77_CHAR_ARG_LEN_DECL); 98 F77_CHAR_ARG_LEN_DECL);
99 99
100 F77_RET_T 100 F77_RET_T
101 F77_FUNC (zptsv, ZPTSV) (const octave_idx_type&, const octave_idx_type&, Complex*, Complex*, 101 F77_FUNC (zptsv, ZPTSV) (const octave_idx_type&, const octave_idx_type&, double*, Complex*,
102 Complex*, const octave_idx_type&, octave_idx_type&); 102 Complex*, const octave_idx_type&, octave_idx_type&);
103 103
104 F77_RET_T 104 F77_RET_T
105 F77_FUNC (zgtsv, ZGTSV) (const octave_idx_type&, const octave_idx_type&, Complex*, Complex*, 105 F77_FUNC (zgtsv, ZGTSV) (const octave_idx_type&, const octave_idx_type&, Complex*, Complex*,
106 Complex*, Complex*, const octave_idx_type&, octave_idx_type&); 106 Complex*, Complex*, const octave_idx_type&, octave_idx_type&);
647 err = 0; 647 err = 0;
648 648
649 // Setup the control parameters 649 // Setup the control parameters
650 Matrix Control (UMFPACK_CONTROL, 1); 650 Matrix Control (UMFPACK_CONTROL, 1);
651 double *control = Control.fortran_vec (); 651 double *control = Control.fortran_vec ();
652 umfpack_zi_defaults (control); 652 UMFPACK_ZNAME (defaults) (control);
653 653
654 double tmp = Voctave_sparse_controls.get_key ("spumoni"); 654 double tmp = Voctave_sparse_controls.get_key ("spumoni");
655 if (!xisnan (tmp)) 655 if (!xisnan (tmp))
656 Control (UMFPACK_PRL) = tmp; 656 Control (UMFPACK_PRL) = tmp;
657 657
668 Control (UMFPACK_FIXQ) = tmp; 668 Control (UMFPACK_FIXQ) = tmp;
669 669
670 // Turn-off UMFPACK scaling for LU 670 // Turn-off UMFPACK scaling for LU
671 Control (UMFPACK_SCALE) = UMFPACK_SCALE_NONE; 671 Control (UMFPACK_SCALE) = UMFPACK_SCALE_NONE;
672 672
673 umfpack_zi_report_control (control); 673 UMFPACK_ZNAME (report_control) (control);
674 674
675 const octave_idx_type *Ap = cidx (); 675 const octave_idx_type *Ap = cidx ();
676 const octave_idx_type *Ai = ridx (); 676 const octave_idx_type *Ai = ridx ();
677 const Complex *Ax = data (); 677 const Complex *Ax = data ();
678 678
679 umfpack_zi_report_matrix (nr, nc, Ap, Ai, 679 UMFPACK_ZNAME (report_matrix) (nr, nc, Ap, Ai,
680 X_CAST (const double *, Ax), 680 X_CAST (const double *, Ax),
681 NULL, 1, control); 681 NULL, 1, control);
682 682
683 void *Symbolic; 683 void *Symbolic;
684 Matrix Info (1, UMFPACK_INFO); 684 Matrix Info (1, UMFPACK_INFO);
685 double *info = Info.fortran_vec (); 685 double *info = Info.fortran_vec ();
686 int status = umfpack_zi_qsymbolic 686 int status = UMFPACK_ZNAME (qsymbolic)
687 (nr, nc, Ap, Ai, X_CAST (const double *, Ax), NULL, 687 (nr, nc, Ap, Ai, X_CAST (const double *, Ax), NULL,
688 NULL, &Symbolic, control, info); 688 NULL, &Symbolic, control, info);
689 689
690 if (status < 0) 690 if (status < 0)
691 { 691 {
692 (*current_liboctave_error_handler) 692 (*current_liboctave_error_handler)
693 ("SparseComplexMatrix::determinant symbolic factorization failed"); 693 ("SparseComplexMatrix::determinant symbolic factorization failed");
694 694
695 umfpack_zi_report_status (control, status); 695 UMFPACK_ZNAME (report_status) (control, status);
696 umfpack_zi_report_info (control, info); 696 UMFPACK_ZNAME (report_info) (control, info);
697 697
698 umfpack_zi_free_symbolic (&Symbolic) ; 698 UMFPACK_ZNAME (free_symbolic) (&Symbolic) ;
699 } 699 }
700 else 700 else
701 { 701 {
702 umfpack_zi_report_symbolic (Symbolic, control); 702 UMFPACK_ZNAME (report_symbolic) (Symbolic, control);
703 703
704 void *Numeric; 704 void *Numeric;
705 status = umfpack_zi_numeric (Ap, Ai, X_CAST (const double *, Ax), 705 status = UMFPACK_ZNAME (numeric) (Ap, Ai,
706 NULL, Symbolic, &Numeric, 706 X_CAST (const double *, Ax), NULL,
707 control, info) ; 707 Symbolic, &Numeric, control, info) ;
708 umfpack_zi_free_symbolic (&Symbolic) ; 708 UMFPACK_ZNAME (free_symbolic) (&Symbolic) ;
709 709
710 rcond = Info (UMFPACK_RCOND); 710 rcond = Info (UMFPACK_RCOND);
711 711
712 if (status < 0) 712 if (status < 0)
713 { 713 {
714 (*current_liboctave_error_handler) 714 (*current_liboctave_error_handler)
715 ("SparseComplexMatrix::determinant numeric factorization failed"); 715 ("SparseComplexMatrix::determinant numeric factorization failed");
716 716
717 umfpack_zi_report_status (control, status); 717 UMFPACK_ZNAME (report_status) (control, status);
718 umfpack_zi_report_info (control, info); 718 UMFPACK_ZNAME (report_info) (control, info);
719 719
720 umfpack_zi_free_numeric (&Numeric); 720 UMFPACK_ZNAME (free_numeric) (&Numeric);
721 } 721 }
722 else 722 else
723 { 723 {
724 umfpack_zi_report_numeric (Numeric, control); 724 UMFPACK_ZNAME (report_numeric) (Numeric, control);
725 725
726 Complex d[2]; 726 Complex d[2];
727 double d_exponent; 727 double d_exponent;
728 728
729 status = umfpack_zi_get_determinant 729 status = UMFPACK_ZNAME (get_determinant)
730 (X_CAST (double *, &d[0]), NULL, &d_exponent, 730 (X_CAST (double *, &d[0]), NULL, &d_exponent,
731 Numeric, info); 731 Numeric, info);
732 d[1] = d_exponent; 732 d[1] = d_exponent;
733 733
734 if (status < 0) 734 if (status < 0)
735 { 735 {
736 (*current_liboctave_error_handler) 736 (*current_liboctave_error_handler)
737 ("SparseComplexMatrix::determinant error calculating determinant"); 737 ("SparseComplexMatrix::determinant error calculating determinant");
738 738
739 umfpack_zi_report_status (control, status); 739 UMFPACK_ZNAME (report_status) (control, status);
740 umfpack_zi_report_info (control, info); 740 UMFPACK_ZNAME (report_info) (control, info);
741 741
742 umfpack_zi_free_numeric (&Numeric); 742 UMFPACK_ZNAME (free_numeric) (&Numeric);
743 } 743 }
744 else 744 else
745 retval = ComplexDET (d); 745 retval = ComplexDET (d);
746 } 746 }
747 } 747 }
1050 anorm = atmp; 1050 anorm = atmp;
1051 } 1051 }
1052 1052
1053 if (typ == SparseType::Permuted_Upper) 1053 if (typ == SparseType::Permuted_Upper)
1054 { 1054 {
1055 retval.resize (b.rows (), b.cols ()); 1055 retval.resize (nr, b_cols);
1056 octave_idx_type *perm = mattype.triangular_perm ();
1056 OCTAVE_LOCAL_BUFFER (Complex, work, nr); 1057 OCTAVE_LOCAL_BUFFER (Complex, work, nr);
1057 octave_idx_type *p_perm = mattype.triangular_row_perm ();
1058 octave_idx_type *q_perm = mattype.triangular_col_perm ();
1059
1060 (*current_liboctave_warning_handler)
1061 ("SparseComplexMatrix::solve XXX FIXME XXX permuted triangular code not tested");
1062 1058
1063 for (octave_idx_type j = 0; j < b_cols; j++) 1059 for (octave_idx_type j = 0; j < b_cols; j++)
1064 { 1060 {
1065 for (octave_idx_type i = 0; i < nr; i++) 1061 for (octave_idx_type i = 0; i < nr; i++)
1066 work[i] = b(i,j); 1062 work[i] = b(i,j);
1067 1063
1068 for (octave_idx_type k = nr-1; k >= 0; k--) 1064 for (octave_idx_type k = nr-1; k >= 0; k--)
1069 { 1065 {
1070 octave_idx_type iidx = q_perm[k]; 1066 octave_idx_type kidx = perm[k];
1071 if (work[iidx] != 0.) 1067
1068 if (work[k] != 0.)
1072 { 1069 {
1073 if (ridx(cidx(iidx+1)-1) != iidx) 1070 if (ridx(cidx(kidx+1)-1) != k)
1074 { 1071 {
1075 err = -2; 1072 err = -2;
1076 goto triangular_error; 1073 goto triangular_error;
1077 } 1074 }
1078 1075
1079 Complex tmp = work[iidx] / data(cidx(iidx+1)-1); 1076 Complex tmp = work[k] / data(cidx(kidx+1)-1);
1080 work[iidx] = tmp; 1077 work[k] = tmp;
1081 for (octave_idx_type i = cidx(iidx); i < cidx(iidx+1)-1; i++) 1078 for (octave_idx_type i = cidx(kidx);
1079 i < cidx(kidx+1)-1; i++)
1082 { 1080 {
1083 octave_idx_type idx2 = q_perm[ridx(i)]; 1081 octave_idx_type iidx = ridx(i);
1084 work[idx2] = 1082 work[iidx] = work[iidx] - tmp * data(i);
1085 work[idx2] - tmp * data(i);
1086 } 1083 }
1087 } 1084 }
1088 } 1085 }
1089 1086
1090 for (octave_idx_type i = 0; i < nr; i++) 1087 for (octave_idx_type i = 0; i < nr; i++)
1091 retval (i, j) = work[p_perm[i]]; 1088 retval (perm[i], j) = work[i];
1092 } 1089 }
1093 1090
1094 // Calculation of 1-norm of inv(*this) 1091 // Calculation of 1-norm of inv(*this)
1095 for (octave_idx_type i = 0; i < nr; i++) 1092 for (octave_idx_type i = 0; i < nr; i++)
1096 work[i] = 0.; 1093 work[i] = 0.;
1097 1094
1098 for (octave_idx_type j = 0; j < nr; j++) 1095 for (octave_idx_type j = 0; j < nr; j++)
1099 { 1096 {
1100 work[q_perm[j]] = 1.; 1097 work[j] = 1.;
1101 1098
1102 for (octave_idx_type k = j; k >= 0; k--) 1099 for (octave_idx_type k = j; k >= 0; k--)
1103 { 1100 {
1104 octave_idx_type iidx = q_perm[k]; 1101 octave_idx_type iidx = perm[k];
1105 1102
1106 if (work[iidx] != 0.) 1103 if (work[k] != 0.)
1107 { 1104 {
1108 Complex tmp = work[iidx] / data(cidx(iidx+1)-1); 1105 Complex tmp = work[k] / data(cidx(iidx+1)-1);
1109 work[iidx] = tmp; 1106 work[k] = tmp;
1110 for (octave_idx_type i = cidx(iidx); i < cidx(iidx+1)-1; i++) 1107 for (octave_idx_type i = cidx(iidx);
1108 i < cidx(iidx+1)-1; i++)
1111 { 1109 {
1112 octave_idx_type idx2 = q_perm[ridx(i)]; 1110 octave_idx_type idx2 = ridx(i);
1113 work[idx2] = work[idx2] - tmp * data(i); 1111 work[idx2] = work[idx2] - tmp * data(i);
1114 } 1112 }
1115 } 1113 }
1116 } 1114 }
1117 double atmp = 0; 1115 double atmp = 0;
1267 octave_idx_type ii = 0; 1265 octave_idx_type ii = 0;
1268 octave_idx_type x_nz = b_nz; 1266 octave_idx_type x_nz = b_nz;
1269 1267
1270 if (typ == SparseType::Permuted_Upper) 1268 if (typ == SparseType::Permuted_Upper)
1271 { 1269 {
1270 octave_idx_type *perm = mattype.triangular_perm ();
1272 OCTAVE_LOCAL_BUFFER (Complex, work, nr); 1271 OCTAVE_LOCAL_BUFFER (Complex, work, nr);
1273 octave_idx_type *p_perm = mattype.triangular_row_perm (); 1272
1274 octave_idx_type *q_perm = mattype.triangular_col_perm (); 1273 OCTAVE_LOCAL_BUFFER (octave_idx_type, rperm, nr);
1275 1274 for (octave_idx_type i = 0; i < nr; i++)
1276 (*current_liboctave_warning_handler) 1275 rperm[perm[i]] = i;
1277 ("SparseComplexMatrix::solve XXX FIXME XXX permuted triangular code not tested");
1278 1276
1279 for (octave_idx_type j = 0; j < b_nc; j++) 1277 for (octave_idx_type j = 0; j < b_nc; j++)
1280 { 1278 {
1281 for (octave_idx_type i = 0; i < nr; i++) 1279 for (octave_idx_type i = 0; i < nr; i++)
1282 work[i] = 0.; 1280 work[i] = 0.;
1283 for (octave_idx_type i = b.cidx(j); i < b.cidx(j+1); i++) 1281 for (octave_idx_type i = b.cidx(j); i < b.cidx(j+1); i++)
1284 work[b.ridx(i)] = b.data(i); 1282 work[b.ridx(i)] = b.data(i);
1285 1283
1286 for (octave_idx_type k = nr-1; k >= 0; k--) 1284 for (octave_idx_type k = nr-1; k >= 0; k--)
1287 { 1285 {
1288 octave_idx_type iidx = q_perm[k]; 1286 octave_idx_type kidx = perm[k];
1289 if (work[iidx] != 0.) 1287
1288 if (work[k] != 0.)
1290 { 1289 {
1291 if (ridx(cidx(iidx+1)-1) != iidx) 1290 if (ridx(cidx(kidx+1)-1) != k)
1292 { 1291 {
1293 err = -2; 1292 err = -2;
1294 goto triangular_error; 1293 goto triangular_error;
1295 } 1294 }
1296 1295
1297 Complex tmp = work[iidx] / data(cidx(iidx+1)-1); 1296 Complex tmp = work[k] / data(cidx(kidx+1)-1);
1298 work[iidx] = tmp; 1297 work[k] = tmp;
1299 for (octave_idx_type i = cidx(iidx); i < cidx(iidx+1)-1; i++) 1298 for (octave_idx_type i = cidx(kidx);
1299 i < cidx(kidx+1)-1; i++)
1300 { 1300 {
1301 octave_idx_type idx2 = q_perm[ridx(i)]; 1301 octave_idx_type iidx = ridx(i);
1302 work[idx2] = 1302 work[iidx] = work[iidx] - tmp * data(i);
1303 work[idx2] - tmp * data(i);
1304 } 1303 }
1305 } 1304 }
1306 } 1305 }
1307 1306
1308 // Count non-zeros in work vector and adjust space in 1307 // Count non-zeros in work vector and adjust space in
1319 retval.change_capacity (sz); 1318 retval.change_capacity (sz);
1320 x_nz = sz; 1319 x_nz = sz;
1321 } 1320 }
1322 1321
1323 for (octave_idx_type i = 0; i < nr; i++) 1322 for (octave_idx_type i = 0; i < nr; i++)
1324 if (work[p_perm[i]] != 0.) 1323 if (work[rperm[i]] != 0.)
1325 { 1324 {
1326 retval.xridx(ii) = i; 1325 retval.xridx(ii) = i;
1327 retval.xdata(ii++) = work[p_perm[i]]; 1326 retval.xdata(ii++) = work[rperm[i]];
1328 } 1327 }
1329 retval.xcidx(j+1) = ii; 1328 retval.xcidx(j+1) = ii;
1330 } 1329 }
1331 1330
1332 retval.maybe_compress (); 1331 retval.maybe_compress ();
1335 for (octave_idx_type i = 0; i < nr; i++) 1334 for (octave_idx_type i = 0; i < nr; i++)
1336 work[i] = 0.; 1335 work[i] = 0.;
1337 1336
1338 for (octave_idx_type j = 0; j < nr; j++) 1337 for (octave_idx_type j = 0; j < nr; j++)
1339 { 1338 {
1340 work[q_perm[j]] = 1.; 1339 work[j] = 1.;
1341 1340
1342 for (octave_idx_type k = j; k >= 0; k--) 1341 for (octave_idx_type k = j; k >= 0; k--)
1343 { 1342 {
1344 octave_idx_type iidx = q_perm[k]; 1343 octave_idx_type iidx = perm[k];
1345 1344
1346 if (work[iidx] != 0.) 1345 if (work[k] != 0.)
1347 { 1346 {
1348 Complex tmp = work[iidx] / data(cidx(iidx+1)-1); 1347 Complex tmp = work[k] / data(cidx(iidx+1)-1);
1349 work[iidx] = tmp; 1348 work[k] = tmp;
1350 for (octave_idx_type i = cidx(iidx); i < cidx(iidx+1)-1; i++) 1349 for (octave_idx_type i = cidx(iidx);
1350 i < cidx(iidx+1)-1; i++)
1351 { 1351 {
1352 octave_idx_type idx2 = q_perm[ridx(i)]; 1352 octave_idx_type idx2 = ridx(i);
1353 work[idx2] = work[idx2] - tmp * data(i); 1353 work[idx2] = work[idx2] - tmp * data(i);
1354 } 1354 }
1355 } 1355 }
1356 } 1356 }
1357 double atmp = 0; 1357 double atmp = 0;
1524 anorm = atmp; 1524 anorm = atmp;
1525 } 1525 }
1526 1526
1527 if (typ == SparseType::Permuted_Upper) 1527 if (typ == SparseType::Permuted_Upper)
1528 { 1528 {
1529 retval.resize (b.rows (), b.cols ()); 1529 retval.resize (nr, b_nc);
1530 octave_idx_type *perm = mattype.triangular_perm ();
1530 OCTAVE_LOCAL_BUFFER (Complex, work, nr); 1531 OCTAVE_LOCAL_BUFFER (Complex, work, nr);
1531 octave_idx_type *p_perm = mattype.triangular_row_perm ();
1532 octave_idx_type *q_perm = mattype.triangular_col_perm ();
1533
1534 (*current_liboctave_warning_handler)
1535 ("SparseComplexMatrix::solve XXX FIXME XXX permuted triangular code not tested");
1536 1532
1537 for (octave_idx_type j = 0; j < b_nc; j++) 1533 for (octave_idx_type j = 0; j < b_nc; j++)
1538 { 1534 {
1539 for (octave_idx_type i = 0; i < nr; i++) 1535 for (octave_idx_type i = 0; i < nr; i++)
1540 work[i] = b(i,j); 1536 work[i] = b(i,j);
1541 1537
1542 for (octave_idx_type k = nr-1; k >= 0; k--) 1538 for (octave_idx_type k = nr-1; k >= 0; k--)
1543 { 1539 {
1544 octave_idx_type iidx = q_perm[k]; 1540 octave_idx_type kidx = perm[k];
1545 if (work[iidx] != 0.) 1541
1542 if (work[k] != 0.)
1546 { 1543 {
1547 if (ridx(cidx(iidx+1)-1) != iidx) 1544 if (ridx(cidx(kidx+1)-1) != k)
1548 { 1545 {
1549 err = -2; 1546 err = -2;
1550 goto triangular_error; 1547 goto triangular_error;
1551 } 1548 }
1552 1549
1553 Complex tmp = work[iidx] / data(cidx(iidx+1)-1); 1550 Complex tmp = work[k] / data(cidx(kidx+1)-1);
1554 work[iidx] = tmp; 1551 work[k] = tmp;
1555 for (octave_idx_type i = cidx(iidx); i < cidx(iidx+1)-1; i++) 1552 for (octave_idx_type i = cidx(kidx);
1553 i < cidx(kidx+1)-1; i++)
1556 { 1554 {
1557 octave_idx_type idx2 = q_perm[ridx(i)]; 1555 octave_idx_type iidx = ridx(i);
1558 work[idx2] = 1556 work[iidx] = work[iidx] - tmp * data(i);
1559 work[idx2] - tmp * data(i);
1560 } 1557 }
1561 } 1558 }
1562 } 1559 }
1563 1560
1564 for (octave_idx_type i = 0; i < nr; i++) 1561 for (octave_idx_type i = 0; i < nr; i++)
1565 retval (i, j) = work[p_perm[i]]; 1562 retval (perm[i], j) = work[i];
1566
1567 } 1563 }
1568 1564
1569 // Calculation of 1-norm of inv(*this) 1565 // Calculation of 1-norm of inv(*this)
1570 for (octave_idx_type i = 0; i < nr; i++) 1566 for (octave_idx_type i = 0; i < nr; i++)
1571 work[i] = 0.; 1567 work[i] = 0.;
1572 1568
1573 for (octave_idx_type j = 0; j < nr; j++) 1569 for (octave_idx_type j = 0; j < nr; j++)
1574 { 1570 {
1575 work[q_perm[j]] = 1.; 1571 work[j] = 1.;
1576 1572
1577 for (octave_idx_type k = j; k >= 0; k--) 1573 for (octave_idx_type k = j; k >= 0; k--)
1578 { 1574 {
1579 octave_idx_type iidx = q_perm[k]; 1575 octave_idx_type iidx = perm[k];
1580 1576
1581 if (work[iidx] != 0.) 1577 if (work[k] != 0.)
1582 { 1578 {
1583 Complex tmp = work[iidx] / data(cidx(iidx+1)-1); 1579 Complex tmp = work[k] / data(cidx(iidx+1)-1);
1584 work[iidx] = tmp; 1580 work[k] = tmp;
1585 for (octave_idx_type i = cidx(iidx); i < cidx(iidx+1)-1; i++) 1581 for (octave_idx_type i = cidx(iidx);
1582 i < cidx(iidx+1)-1; i++)
1586 { 1583 {
1587 octave_idx_type idx2 = q_perm[ridx(i)]; 1584 octave_idx_type idx2 = ridx(i);
1588 work[idx2] = work[idx2] - tmp * data(i); 1585 work[idx2] = work[idx2] - tmp * data(i);
1589 } 1586 }
1590 } 1587 }
1591 } 1588 }
1592 double atmp = 0; 1589 double atmp = 0;
1742 octave_idx_type ii = 0; 1739 octave_idx_type ii = 0;
1743 octave_idx_type x_nz = b_nz; 1740 octave_idx_type x_nz = b_nz;
1744 1741
1745 if (typ == SparseType::Permuted_Upper) 1742 if (typ == SparseType::Permuted_Upper)
1746 { 1743 {
1744 octave_idx_type *perm = mattype.triangular_perm ();
1747 OCTAVE_LOCAL_BUFFER (Complex, work, nr); 1745 OCTAVE_LOCAL_BUFFER (Complex, work, nr);
1748 octave_idx_type *p_perm = mattype.triangular_row_perm (); 1746
1749 octave_idx_type *q_perm = mattype.triangular_col_perm (); 1747 OCTAVE_LOCAL_BUFFER (octave_idx_type, rperm, nr);
1750 1748 for (octave_idx_type i = 0; i < nr; i++)
1751 (*current_liboctave_warning_handler) 1749 rperm[perm[i]] = i;
1752 ("SparseComplexMatrix::solve XXX FIXME XXX permuted triangular code not tested");
1753 1750
1754 for (octave_idx_type j = 0; j < b_nc; j++) 1751 for (octave_idx_type j = 0; j < b_nc; j++)
1755 { 1752 {
1756 for (octave_idx_type i = 0; i < nr; i++) 1753 for (octave_idx_type i = 0; i < nr; i++)
1757 work[i] = 0.; 1754 work[i] = 0.;
1758 for (octave_idx_type i = b.cidx(j); i < b.cidx(j+1); i++) 1755 for (octave_idx_type i = b.cidx(j); i < b.cidx(j+1); i++)
1759 work[b.ridx(i)] = b.data(i); 1756 work[b.ridx(i)] = b.data(i);
1760 1757
1761 for (octave_idx_type k = nr-1; k >= 0; k--) 1758 for (octave_idx_type k = nr-1; k >= 0; k--)
1762 { 1759 {
1763 octave_idx_type iidx = q_perm[k]; 1760 octave_idx_type kidx = perm[k];
1764 if (work[iidx] != 0.) 1761
1762 if (work[k] != 0.)
1765 { 1763 {
1766 if (ridx(cidx(iidx+1)-1) != iidx) 1764 if (ridx(cidx(kidx+1)-1) != k)
1767 { 1765 {
1768 err = -2; 1766 err = -2;
1769 goto triangular_error; 1767 goto triangular_error;
1770 } 1768 }
1771 1769
1772 Complex tmp = work[iidx] / data(cidx(iidx+1)-1); 1770 Complex tmp = work[k] / data(cidx(kidx+1)-1);
1773 work[iidx] = tmp; 1771 work[k] = tmp;
1774 for (octave_idx_type i = cidx(iidx); i < cidx(iidx+1)-1; i++) 1772 for (octave_idx_type i = cidx(kidx);
1773 i < cidx(kidx+1)-1; i++)
1775 { 1774 {
1776 octave_idx_type idx2 = q_perm[ridx(i)]; 1775 octave_idx_type iidx = ridx(i);
1777 work[idx2] = 1776 work[iidx] = work[iidx] - tmp * data(i);
1778 work[idx2] - tmp * data(i);
1779 } 1777 }
1780 } 1778 }
1781 } 1779 }
1782 1780
1783 // Count non-zeros in work vector and adjust space in 1781 // Count non-zeros in work vector and adjust space in
1794 retval.change_capacity (sz); 1792 retval.change_capacity (sz);
1795 x_nz = sz; 1793 x_nz = sz;
1796 } 1794 }
1797 1795
1798 for (octave_idx_type i = 0; i < nr; i++) 1796 for (octave_idx_type i = 0; i < nr; i++)
1799 if (work[p_perm[i]] != 0.) 1797 if (work[rperm[i]] != 0.)
1800 { 1798 {
1801 retval.xridx(ii) = i; 1799 retval.xridx(ii) = i;
1802 retval.xdata(ii++) = work[p_perm[i]]; 1800 retval.xdata(ii++) = work[rperm[i]];
1803 } 1801 }
1804 retval.xcidx(j+1) = ii; 1802 retval.xcidx(j+1) = ii;
1805 } 1803 }
1806 1804
1807 retval.maybe_compress (); 1805 retval.maybe_compress ();
1810 for (octave_idx_type i = 0; i < nr; i++) 1808 for (octave_idx_type i = 0; i < nr; i++)
1811 work[i] = 0.; 1809 work[i] = 0.;
1812 1810
1813 for (octave_idx_type j = 0; j < nr; j++) 1811 for (octave_idx_type j = 0; j < nr; j++)
1814 { 1812 {
1815 work[q_perm[j]] = 1.; 1813 work[j] = 1.;
1816 1814
1817 for (octave_idx_type k = j; k >= 0; k--) 1815 for (octave_idx_type k = j; k >= 0; k--)
1818 { 1816 {
1819 octave_idx_type iidx = q_perm[k]; 1817 octave_idx_type iidx = perm[k];
1820 1818
1821 if (work[iidx] != 0.) 1819 if (work[k] != 0.)
1822 { 1820 {
1823 Complex tmp = work[iidx] / data(cidx(iidx+1)-1); 1821 Complex tmp = work[k] / data(cidx(iidx+1)-1);
1824 work[iidx] = tmp; 1822 work[k] = tmp;
1825 for (octave_idx_type i = cidx(iidx); i < cidx(iidx+1)-1; i++) 1823 for (octave_idx_type i = cidx(iidx);
1824 i < cidx(iidx+1)-1; i++)
1826 { 1825 {
1827 octave_idx_type idx2 = q_perm[ridx(i)]; 1826 octave_idx_type idx2 = ridx(i);
1828 work[idx2] = work[idx2] - tmp * data(i); 1827 work[idx2] = work[idx2] - tmp * data(i);
1829 } 1828 }
1830 } 1829 }
1831 } 1830 }
1832 double atmp = 0; 1831 double atmp = 0;
2001 2000
2002 if (typ == SparseType::Permuted_Lower) 2001 if (typ == SparseType::Permuted_Lower)
2003 { 2002 {
2004 retval.resize (b.rows (), b.cols ()); 2003 retval.resize (b.rows (), b.cols ());
2005 OCTAVE_LOCAL_BUFFER (Complex, work, nr); 2004 OCTAVE_LOCAL_BUFFER (Complex, work, nr);
2006 octave_idx_type *p_perm = mattype.triangular_row_perm (); 2005 octave_idx_type *perm = mattype.triangular_perm ();
2007 octave_idx_type *q_perm = mattype.triangular_col_perm ();
2008
2009 (*current_liboctave_warning_handler)
2010 ("SparseComplexMatrix::solve XXX FIXME XXX permuted triangular code not tested");
2011 2006
2012 for (octave_idx_type j = 0; j < b_cols; j++) 2007 for (octave_idx_type j = 0; j < b_cols; j++)
2013 { 2008 {
2014 for (octave_idx_type i = 0; i < nr; i++) 2009 for (octave_idx_type i = 0; i < nr; i++)
2015 work[i] = b(i,j); 2010 work[perm[i]] = b(i,j);
2016 2011
2017 for (octave_idx_type k = 0; k < nr; k++) 2012 for (octave_idx_type k = 0; k < nr; k++)
2018 { 2013 {
2019 octave_idx_type iidx = q_perm[k]; 2014 if (work[k] != 0.)
2020 if (work[iidx] != 0.)
2021 { 2015 {
2022 if (ridx(cidx(iidx)) != iidx) 2016 octave_idx_type minr = nr;
2017 octave_idx_type mini = 0;
2018
2019 for (octave_idx_type i = cidx(k); i < cidx(k+1); i++)
2020 if (perm[ridx(i)] < minr)
2021 {
2022 minr = perm[ridx(i)];
2023 mini = i;
2024 }
2025
2026 if (minr != k)
2023 { 2027 {
2024 err = -2; 2028 err = -2;
2025 goto triangular_error; 2029 goto triangular_error;
2026 } 2030 }
2027 2031
2028 Complex tmp = work[iidx] / data(cidx(iidx+1)-1); 2032 Complex tmp = work[k] / data(mini);
2029 work[iidx] = tmp; 2033 work[k] = tmp;
2030 for (octave_idx_type i = cidx(iidx)+1; i < cidx(iidx+1); i++) 2034 for (octave_idx_type i = cidx(k); i < cidx(k+1); i++)
2031 { 2035 {
2032 octave_idx_type idx2 = q_perm[ridx(i)]; 2036 if (i == mini)
2033 work[idx2] = 2037 continue;
2034 work[idx2] - tmp * data(i); 2038
2039 octave_idx_type iidx = perm[ridx(i)];
2040 work[iidx] = work[iidx] - tmp * data(i);
2035 } 2041 }
2036 } 2042 }
2037 } 2043 }
2038 2044
2039 for (octave_idx_type i = 0; i < nr; i++) 2045 for (octave_idx_type i = 0; i < nr; i++)
2040 retval (i, j) = work[p_perm[i]]; 2046 retval (i, j) = work[i];
2041
2042 } 2047 }
2043 2048
2044 // Calculation of 1-norm of inv(*this) 2049 // Calculation of 1-norm of inv(*this)
2045 for (octave_idx_type i = 0; i < nr; i++) 2050 for (octave_idx_type i = 0; i < nr; i++)
2046 work[i] = 0.; 2051 work[i] = 0.;
2047 2052
2048 for (octave_idx_type j = 0; j < nr; j++) 2053 for (octave_idx_type j = 0; j < nr; j++)
2049 { 2054 {
2050 work[q_perm[j]] = 1.; 2055 work[j] = 1.;
2051 2056
2052 for (octave_idx_type k = 0; k < nr; k++) 2057 for (octave_idx_type k = 0; k < nr; k++)
2053 { 2058 {
2054 octave_idx_type iidx = q_perm[k]; 2059 if (work[k] != 0.)
2055
2056 if (work[iidx] != 0.)
2057 { 2060 {
2058 Complex tmp = work[iidx] / data(cidx(iidx+1)-1); 2061 octave_idx_type minr = nr;
2059 work[iidx] = tmp; 2062 octave_idx_type mini = 0;
2060 for (octave_idx_type i = cidx(iidx)+1; i < cidx(iidx+1); i++) 2063
2064 for (octave_idx_type i = cidx(k); i < cidx(k+1); i++)
2065 if (perm[ridx(i)] < minr)
2066 {
2067 minr = perm[ridx(i)];
2068 mini = i;
2069 }
2070
2071 Complex tmp = work[k] / data(mini);
2072 work[k] = tmp;
2073 for (octave_idx_type i = cidx(k); i < cidx(k+1); i++)
2061 { 2074 {
2062 octave_idx_type idx2 = q_perm[ridx(i)]; 2075 if (i == mini)
2063 work[idx2] = work[idx2] - tmp * data(i); 2076 continue;
2077
2078 octave_idx_type iidx = perm[ridx(i)];
2079 work[iidx] = work[iidx] - tmp * data(i);
2064 } 2080 }
2065 } 2081 }
2066 } 2082 }
2083
2067 double atmp = 0; 2084 double atmp = 0;
2068 for (octave_idx_type i = 0; i < j+1; i++) 2085 for (octave_idx_type i = j; i < nr; i++)
2069 { 2086 {
2070 atmp += std::abs(work[i]); 2087 atmp += std::abs(work[i]);
2071 work[i] = 0.; 2088 work[i] = 0.;
2072 } 2089 }
2073 if (atmp > ainvnorm) 2090 if (atmp > ainvnorm)
2219 octave_idx_type x_nz = b_nz; 2236 octave_idx_type x_nz = b_nz;
2220 2237
2221 if (typ == SparseType::Permuted_Lower) 2238 if (typ == SparseType::Permuted_Lower)
2222 { 2239 {
2223 OCTAVE_LOCAL_BUFFER (Complex, work, nr); 2240 OCTAVE_LOCAL_BUFFER (Complex, work, nr);
2224 octave_idx_type *p_perm = mattype.triangular_row_perm (); 2241 octave_idx_type *perm = mattype.triangular_perm ();
2225 octave_idx_type *q_perm = mattype.triangular_col_perm ();
2226
2227 (*current_liboctave_warning_handler)
2228 ("SparseComplexMatrix::solve XXX FIXME XXX permuted triangular code not tested");
2229 2242
2230 for (octave_idx_type j = 0; j < b_nc; j++) 2243 for (octave_idx_type j = 0; j < b_nc; j++)
2231 { 2244 {
2232 for (octave_idx_type i = 0; i < nr; i++) 2245 for (octave_idx_type i = 0; i < nr; i++)
2233 work[i] = 0.; 2246 work[i] = 0.;
2234 for (octave_idx_type i = b.cidx(j); i < b.cidx(j+1); i++) 2247 for (octave_idx_type i = b.cidx(j); i < b.cidx(j+1); i++)
2235 work[b.ridx(i)] = b.data(i); 2248 work[perm[b.ridx(i)]] = b.data(i);
2236 2249
2237 for (octave_idx_type k = 0; k < nr; k++) 2250 for (octave_idx_type k = 0; k < nr; k++)
2238 { 2251 {
2239 octave_idx_type iidx = q_perm[k]; 2252 if (work[k] != 0.)
2240 if (work[iidx] != 0.)
2241 { 2253 {
2242 if (ridx(cidx(iidx)) != iidx) 2254 octave_idx_type minr = nr;
2255 octave_idx_type mini = 0;
2256
2257 for (octave_idx_type i = cidx(k); i < cidx(k+1); i++)
2258 if (perm[ridx(i)] < minr)
2259 {
2260 minr = perm[ridx(i)];
2261 mini = i;
2262 }
2263
2264 if (minr != k)
2243 { 2265 {
2244 err = -2; 2266 err = -2;
2245 goto triangular_error; 2267 goto triangular_error;
2246 } 2268 }
2247 2269
2248 Complex tmp = work[iidx] / data(cidx(iidx+1)-1); 2270 Complex tmp = work[k] / data(mini);
2249 work[iidx] = tmp; 2271 work[k] = tmp;
2250 for (octave_idx_type i = cidx(iidx)+1; i < cidx(iidx+1); i++) 2272 for (octave_idx_type i = cidx(k); i < cidx(k+1); i++)
2251 { 2273 {
2252 octave_idx_type idx2 = q_perm[ridx(i)]; 2274 if (i == mini)
2253 work[idx2] = 2275 continue;
2254 work[idx2] - tmp * data(i); 2276
2277 octave_idx_type iidx = perm[ridx(i)];
2278 work[iidx] = work[iidx] - tmp * data(i);
2255 } 2279 }
2256 } 2280 }
2257 } 2281 }
2258 2282
2259 // Count non-zeros in work vector and adjust space in 2283 // Count non-zeros in work vector and adjust space in
2270 retval.change_capacity (sz); 2294 retval.change_capacity (sz);
2271 x_nz = sz; 2295 x_nz = sz;
2272 } 2296 }
2273 2297
2274 for (octave_idx_type i = 0; i < nr; i++) 2298 for (octave_idx_type i = 0; i < nr; i++)
2275 if (work[p_perm[i]] != 0.) 2299 if (work[i] != 0.)
2276 { 2300 {
2277 retval.xridx(ii) = i; 2301 retval.xridx(ii) = i;
2278 retval.xdata(ii++) = work[p_perm[i]]; 2302 retval.xdata(ii++) = work[i];
2279 } 2303 }
2280 retval.xcidx(j+1) = ii; 2304 retval.xcidx(j+1) = ii;
2281 } 2305 }
2282 2306
2283 retval.maybe_compress (); 2307 retval.maybe_compress ();
2286 for (octave_idx_type i = 0; i < nr; i++) 2310 for (octave_idx_type i = 0; i < nr; i++)
2287 work[i] = 0.; 2311 work[i] = 0.;
2288 2312
2289 for (octave_idx_type j = 0; j < nr; j++) 2313 for (octave_idx_type j = 0; j < nr; j++)
2290 { 2314 {
2291 work[q_perm[j]] = 1.; 2315 work[j] = 1.;
2292 2316
2293 for (octave_idx_type k = 0; k < nr; k++) 2317 for (octave_idx_type k = 0; k < nr; k++)
2294 { 2318 {
2295 octave_idx_type iidx = q_perm[k]; 2319 if (work[k] != 0.)
2296
2297 if (work[iidx] != 0.)
2298 { 2320 {
2299 Complex tmp = work[iidx] / data(cidx(iidx+1)-1); 2321 octave_idx_type minr = nr;
2300 work[iidx] = tmp; 2322 octave_idx_type mini = 0;
2301 for (octave_idx_type i = cidx(iidx)+1; i < cidx(iidx+1); i++) 2323
2324 for (octave_idx_type i = cidx(k); i < cidx(k+1); i++)
2325 if (perm[ridx(i)] < minr)
2326 {
2327 minr = perm[ridx(i)];
2328 mini = i;
2329 }
2330
2331 Complex tmp = work[k] / data(mini);
2332 work[k] = tmp;
2333 for (octave_idx_type i = cidx(k); i < cidx(k+1); i++)
2302 { 2334 {
2303 octave_idx_type idx2 = q_perm[ridx(i)]; 2335 if (i == mini)
2304 work[idx2] = work[idx2] - tmp * data(i); 2336 continue;
2337
2338 octave_idx_type iidx = perm[ridx(i)];
2339 work[iidx] = work[iidx] - tmp * data(i);
2305 } 2340 }
2306 } 2341 }
2307 } 2342 }
2343
2308 double atmp = 0; 2344 double atmp = 0;
2309 for (octave_idx_type i = 0; i < j+1; i++) 2345 for (octave_idx_type i = j; i < nr; i++)
2310 { 2346 {
2311 atmp += std::abs(work[i]); 2347 atmp += std::abs(work[i]);
2312 work[i] = 0.; 2348 work[i] = 0.;
2313 } 2349 }
2314 if (atmp > ainvnorm) 2350 if (atmp > ainvnorm)
2480 2516
2481 if (typ == SparseType::Permuted_Lower) 2517 if (typ == SparseType::Permuted_Lower)
2482 { 2518 {
2483 retval.resize (b.rows (), b.cols ()); 2519 retval.resize (b.rows (), b.cols ());
2484 OCTAVE_LOCAL_BUFFER (Complex, work, nr); 2520 OCTAVE_LOCAL_BUFFER (Complex, work, nr);
2485 octave_idx_type *p_perm = mattype.triangular_row_perm (); 2521 octave_idx_type *perm = mattype.triangular_perm ();
2486 octave_idx_type *q_perm = mattype.triangular_col_perm ();
2487
2488 (*current_liboctave_warning_handler)
2489 ("SparseComplexMatrix::solve XXX FIXME XXX permuted triangular code not tested");
2490 2522
2491 for (octave_idx_type j = 0; j < b_nc; j++) 2523 for (octave_idx_type j = 0; j < b_nc; j++)
2492 { 2524 {
2493 for (octave_idx_type i = 0; i < nr; i++) 2525 for (octave_idx_type i = 0; i < nr; i++)
2494 work[i] = b(i,j); 2526 work[perm[i]] = b(i,j);
2495 2527
2496 for (octave_idx_type k = 0; k < nr; k++) 2528 for (octave_idx_type k = 0; k < nr; k++)
2497 { 2529 {
2498 octave_idx_type iidx = q_perm[k]; 2530 if (work[k] != 0.)
2499 if (work[iidx] != 0.)
2500 { 2531 {
2501 if (ridx(cidx(iidx)) != iidx) 2532 octave_idx_type minr = nr;
2533 octave_idx_type mini = 0;
2534
2535 for (octave_idx_type i = cidx(k); i < cidx(k+1); i++)
2536 if (perm[ridx(i)] < minr)
2537 {
2538 minr = perm[ridx(i)];
2539 mini = i;
2540 }
2541
2542 if (minr != k)
2502 { 2543 {
2503 err = -2; 2544 err = -2;
2504 goto triangular_error; 2545 goto triangular_error;
2505 } 2546 }
2506 2547
2507 Complex tmp = work[iidx] / data(cidx(iidx+1)-1); 2548 Complex tmp = work[k] / data(mini);
2508 work[iidx] = tmp; 2549 work[k] = tmp;
2509 for (octave_idx_type i = cidx(iidx)+1; i < cidx(iidx+1); i++) 2550 for (octave_idx_type i = cidx(k); i < cidx(k+1); i++)
2510 { 2551 {
2511 octave_idx_type idx2 = q_perm[ridx(i)]; 2552 if (i == mini)
2512 work[idx2] = 2553 continue;
2513 work[idx2] - tmp * data(i); 2554
2555 octave_idx_type iidx = perm[ridx(i)];
2556 work[iidx] = work[iidx] - tmp * data(i);
2514 } 2557 }
2515 } 2558 }
2516 } 2559 }
2517 2560
2518 for (octave_idx_type i = 0; i < nr; i++) 2561 for (octave_idx_type i = 0; i < nr; i++)
2519 retval (i, j) = work[p_perm[i]]; 2562 retval (i, j) = work[i];
2520
2521 } 2563 }
2522 2564
2523 // Calculation of 1-norm of inv(*this) 2565 // Calculation of 1-norm of inv(*this)
2524 for (octave_idx_type i = 0; i < nr; i++) 2566 for (octave_idx_type i = 0; i < nr; i++)
2525 work[i] = 0.; 2567 work[i] = 0.;
2526 2568
2527 for (octave_idx_type j = 0; j < nr; j++) 2569 for (octave_idx_type j = 0; j < nr; j++)
2528 { 2570 {
2529 work[q_perm[j]] = 1.; 2571 work[j] = 1.;
2530 2572
2531 for (octave_idx_type k = 0; k < nr; k++) 2573 for (octave_idx_type k = 0; k < nr; k++)
2532 { 2574 {
2533 octave_idx_type iidx = q_perm[k]; 2575 if (work[k] != 0.)
2534
2535 if (work[iidx] != 0.)
2536 { 2576 {
2537 Complex tmp = work[iidx] / data(cidx(iidx+1)-1); 2577 octave_idx_type minr = nr;
2538 work[iidx] = tmp; 2578 octave_idx_type mini = 0;
2539 for (octave_idx_type i = cidx(iidx)+1; i < cidx(iidx+1); i++) 2579
2580 for (octave_idx_type i = cidx(k); i < cidx(k+1); i++)
2581 if (perm[ridx(i)] < minr)
2582 {
2583 minr = perm[ridx(i)];
2584 mini = i;
2585 }
2586
2587 Complex tmp = work[k] / data(mini);
2588 work[k] = tmp;
2589 for (octave_idx_type i = cidx(k); i < cidx(k+1); i++)
2540 { 2590 {
2541 octave_idx_type idx2 = q_perm[ridx(i)]; 2591 if (i == mini)
2542 work[idx2] = work[idx2] - tmp * data(i); 2592 continue;
2593
2594 octave_idx_type iidx = perm[ridx(i)];
2595 work[iidx] = work[iidx] - tmp * data(i);
2543 } 2596 }
2544 } 2597 }
2545 } 2598 }
2599
2546 double atmp = 0; 2600 double atmp = 0;
2547 for (octave_idx_type i = 0; i < j+1; i++) 2601 for (octave_idx_type i = j; i < nr; i++)
2548 { 2602 {
2549 atmp += std::abs(work[i]); 2603 atmp += std::abs(work[i]);
2550 work[i] = 0.; 2604 work[i] = 0.;
2551 } 2605 }
2552 if (atmp > ainvnorm) 2606 if (atmp > ainvnorm)
2699 octave_idx_type x_nz = b_nz; 2753 octave_idx_type x_nz = b_nz;
2700 2754
2701 if (typ == SparseType::Permuted_Lower) 2755 if (typ == SparseType::Permuted_Lower)
2702 { 2756 {
2703 OCTAVE_LOCAL_BUFFER (Complex, work, nr); 2757 OCTAVE_LOCAL_BUFFER (Complex, work, nr);
2704 octave_idx_type *p_perm = mattype.triangular_row_perm (); 2758 octave_idx_type *perm = mattype.triangular_perm ();
2705 octave_idx_type *q_perm = mattype.triangular_col_perm ();
2706
2707 (*current_liboctave_warning_handler)
2708 ("SparseComplexMatrix::solve XXX FIXME XXX permuted triangular code not tested");
2709 2759
2710 for (octave_idx_type j = 0; j < b_nc; j++) 2760 for (octave_idx_type j = 0; j < b_nc; j++)
2711 { 2761 {
2712 for (octave_idx_type i = 0; i < nr; i++) 2762 for (octave_idx_type i = 0; i < nr; i++)
2713 work[i] = 0.; 2763 work[i] = 0.;
2714 for (octave_idx_type i = b.cidx(j); i < b.cidx(j+1); i++) 2764 for (octave_idx_type i = b.cidx(j); i < b.cidx(j+1); i++)
2715 work[b.ridx(i)] = b.data(i); 2765 work[perm[b.ridx(i)]] = b.data(i);
2716 2766
2717 for (octave_idx_type k = 0; k < nr; k++) 2767 for (octave_idx_type k = 0; k < nr; k++)
2718 { 2768 {
2719 octave_idx_type iidx = q_perm[k]; 2769 if (work[k] != 0.)
2720 if (work[iidx] != 0.)
2721 { 2770 {
2722 if (ridx(cidx(iidx)) != iidx) 2771 octave_idx_type minr = nr;
2772 octave_idx_type mini = 0;
2773
2774 for (octave_idx_type i = cidx(k); i < cidx(k+1); i++)
2775 if (perm[ridx(i)] < minr)
2776 {
2777 minr = perm[ridx(i)];
2778 mini = i;
2779 }
2780
2781 if (minr != k)
2723 { 2782 {
2724 err = -2; 2783 err = -2;
2725 goto triangular_error; 2784 goto triangular_error;
2726 } 2785 }
2727 2786
2728 Complex tmp = work[iidx] / data(cidx(iidx+1)-1); 2787 Complex tmp = work[k] / data(mini);
2729 work[iidx] = tmp; 2788 work[k] = tmp;
2730 for (octave_idx_type i = cidx(iidx)+1; i < cidx(iidx+1); i++) 2789 for (octave_idx_type i = cidx(k); i < cidx(k+1); i++)
2731 { 2790 {
2732 octave_idx_type idx2 = q_perm[ridx(i)]; 2791 if (i == mini)
2733 work[idx2] = 2792 continue;
2734 work[idx2] - tmp * data(i); 2793
2794 octave_idx_type iidx = perm[ridx(i)];
2795 work[iidx] = work[iidx] - tmp * data(i);
2735 } 2796 }
2736 } 2797 }
2737 } 2798 }
2738 2799
2739 // Count non-zeros in work vector and adjust space in 2800 // Count non-zeros in work vector and adjust space in
2750 retval.change_capacity (sz); 2811 retval.change_capacity (sz);
2751 x_nz = sz; 2812 x_nz = sz;
2752 } 2813 }
2753 2814
2754 for (octave_idx_type i = 0; i < nr; i++) 2815 for (octave_idx_type i = 0; i < nr; i++)
2755 if (work[p_perm[i]] != 0.) 2816 if (work[i] != 0.)
2756 { 2817 {
2757 retval.xridx(ii) = i; 2818 retval.xridx(ii) = i;
2758 retval.xdata(ii++) = work[p_perm[i]]; 2819 retval.xdata(ii++) = work[i];
2759 } 2820 }
2760 retval.xcidx(j+1) = ii; 2821 retval.xcidx(j+1) = ii;
2761 } 2822 }
2762 2823
2763 retval.maybe_compress (); 2824 retval.maybe_compress ();
2766 for (octave_idx_type i = 0; i < nr; i++) 2827 for (octave_idx_type i = 0; i < nr; i++)
2767 work[i] = 0.; 2828 work[i] = 0.;
2768 2829
2769 for (octave_idx_type j = 0; j < nr; j++) 2830 for (octave_idx_type j = 0; j < nr; j++)
2770 { 2831 {
2771 work[q_perm[j]] = 1.; 2832 work[j] = 1.;
2772 2833
2773 for (octave_idx_type k = 0; k < nr; k++) 2834 for (octave_idx_type k = 0; k < nr; k++)
2774 { 2835 {
2775 octave_idx_type iidx = q_perm[k]; 2836 if (work[k] != 0.)
2776
2777 if (work[iidx] != 0.)
2778 { 2837 {
2779 Complex tmp = work[iidx] / data(cidx(iidx+1)-1); 2838 octave_idx_type minr = nr;
2780 work[iidx] = tmp; 2839 octave_idx_type mini = 0;
2781 for (octave_idx_type i = cidx(iidx)+1; i < cidx(iidx+1); i++) 2840
2841 for (octave_idx_type i = cidx(k); i < cidx(k+1); i++)
2842 if (perm[ridx(i)] < minr)
2843 {
2844 minr = perm[ridx(i)];
2845 mini = i;
2846 }
2847
2848 Complex tmp = work[k] / data(mini);
2849 work[k] = tmp;
2850 for (octave_idx_type i = cidx(k); i < cidx(k+1); i++)
2782 { 2851 {
2783 octave_idx_type idx2 = q_perm[ridx(i)]; 2852 if (i == mini)
2784 work[idx2] = work[idx2] - tmp * data(i); 2853 continue;
2854
2855 octave_idx_type iidx = perm[ridx(i)];
2856 work[iidx] = work[iidx] - tmp * data(i);
2785 } 2857 }
2786 } 2858 }
2787 } 2859 }
2860
2788 double atmp = 0; 2861 double atmp = 0;
2789 for (octave_idx_type i = 0; i < j+1; i++) 2862 for (octave_idx_type i = j; i < nr; i++)
2790 { 2863 {
2791 atmp += std::abs(work[i]); 2864 atmp += std::abs(work[i]);
2792 work[i] = 0.; 2865 work[i] = 0.;
2793 } 2866 }
2794 if (atmp > ainvnorm) 2867 if (atmp > ainvnorm)
2940 volatile int typ = mattype.type (); 3013 volatile int typ = mattype.type ();
2941 mattype.info (); 3014 mattype.info ();
2942 3015
2943 if (typ == SparseType::Tridiagonal_Hermitian) 3016 if (typ == SparseType::Tridiagonal_Hermitian)
2944 { 3017 {
2945 OCTAVE_LOCAL_BUFFER (Complex, D, nr); 3018 OCTAVE_LOCAL_BUFFER (double, D, nr);
2946 OCTAVE_LOCAL_BUFFER (Complex, DL, nr - 1); 3019 OCTAVE_LOCAL_BUFFER (Complex, DL, nr - 1);
2947 3020
2948 if (mattype.is_dense ()) 3021 if (mattype.is_dense ())
2949 { 3022 {
2950 octave_idx_type ii = 0; 3023 octave_idx_type ii = 0;
2951 3024
2952 for (octave_idx_type j = 0; j < nc-1; j++) 3025 for (octave_idx_type j = 0; j < nc-1; j++)
2953 { 3026 {
2954 D[j] = data(ii++); 3027 D[j] = std::real(data(ii++));
2955 DL[j] = data(ii); 3028 DL[j] = data(ii);
2956 ii += 2; 3029 ii += 2;
2957 } 3030 }
2958 D[nc-1] = data(ii); 3031 D[nc-1] = std::real(data(ii));
2959 } 3032 }
2960 else 3033 else
2961 { 3034 {
2962 D[0] = 0.; 3035 D[0] = 0.;
2963 for (octave_idx_type i = 0; i < nr - 1; i++) 3036 for (octave_idx_type i = 0; i < nr - 1; i++)
2968 3041
2969 for (octave_idx_type j = 0; j < nc; j++) 3042 for (octave_idx_type j = 0; j < nc; j++)
2970 for (octave_idx_type i = cidx(j); i < cidx(j+1); i++) 3043 for (octave_idx_type i = cidx(j); i < cidx(j+1); i++)
2971 { 3044 {
2972 if (ridx(i) == j) 3045 if (ridx(i) == j)
2973 D[j] = data(i); 3046 D[j] = std::real(data(i));
2974 else if (ridx(i) == j + 1) 3047 else if (ridx(i) == j + 1)
2975 DL[j] = data(i); 3048 DL[j] = data(i);
2976 } 3049 }
2977 } 3050 }
2978 3051
3030 if (ridx(i) == j) 3103 if (ridx(i) == j)
3031 D[j] = data(i); 3104 D[j] = data(i);
3032 else if (ridx(i) == j + 1) 3105 else if (ridx(i) == j + 1)
3033 DL[j] = data(i); 3106 DL[j] = data(i);
3034 else if (ridx(i) == j - 1) 3107 else if (ridx(i) == j - 1)
3035 DU[j] = data(i); 3108 DU[j-1] = data(i);
3036 } 3109 }
3037 } 3110 }
3038 3111
3039 octave_idx_type b_nc = b.cols(); 3112 octave_idx_type b_nc = b.cols();
3040 retval = ComplexMatrix (b); 3113 retval = ComplexMatrix (b);
3127 if (ridx(i) == j) 3200 if (ridx(i) == j)
3128 D[j] = data(i); 3201 D[j] = data(i);
3129 else if (ridx(i) == j + 1) 3202 else if (ridx(i) == j + 1)
3130 DL[j] = data(i); 3203 DL[j] = data(i);
3131 else if (ridx(i) == j - 1) 3204 else if (ridx(i) == j - 1)
3132 DU[j] = data(i); 3205 DU[j-1] = data(i);
3133 } 3206 }
3134 } 3207 }
3135 3208
3136 F77_XFCN (zgttrf, ZGTTRF, (nr, DL, D, DU, DU2, pipvt, err)); 3209 F77_XFCN (zgttrf, ZGTTRF, (nr, DL, D, DU, DU2, pipvt, err));
3137 3210
3236 { 3309 {
3237 // Print spparms("spumoni") info if requested 3310 // Print spparms("spumoni") info if requested
3238 volatile int typ = mattype.type (); 3311 volatile int typ = mattype.type ();
3239 mattype.info (); 3312 mattype.info ();
3240 3313
3241 // Note can't treat symmetric case as there is no dpttrf function
3242 if (typ == SparseType::Tridiagonal_Hermitian) 3314 if (typ == SparseType::Tridiagonal_Hermitian)
3243 { 3315 {
3244 OCTAVE_LOCAL_BUFFER (Complex, D, nr); 3316 OCTAVE_LOCAL_BUFFER (double, D, nr);
3245 OCTAVE_LOCAL_BUFFER (Complex, DL, nr - 1); 3317 OCTAVE_LOCAL_BUFFER (Complex, DL, nr - 1);
3246 3318
3247 if (mattype.is_dense ()) 3319 if (mattype.is_dense ())
3248 { 3320 {
3249 octave_idx_type ii = 0; 3321 octave_idx_type ii = 0;
3250 3322
3251 for (octave_idx_type j = 0; j < nc-1; j++) 3323 for (octave_idx_type j = 0; j < nc-1; j++)
3252 { 3324 {
3253 D[j] = data(ii++); 3325 D[j] = std::real(data(ii++));
3254 DL[j] = data(ii); 3326 DL[j] = data(ii);
3255 ii += 2; 3327 ii += 2;
3256 } 3328 }
3257 D[nc-1] = data(ii); 3329 D[nc-1] = std::real(data(ii));
3258 } 3330 }
3259 else 3331 else
3260 { 3332 {
3261 D[0] = 0.; 3333 D[0] = 0.;
3262 for (octave_idx_type i = 0; i < nr - 1; i++) 3334 for (octave_idx_type i = 0; i < nr - 1; i++)
3267 3339
3268 for (octave_idx_type j = 0; j < nc; j++) 3340 for (octave_idx_type j = 0; j < nc; j++)
3269 for (octave_idx_type i = cidx(j); i < cidx(j+1); i++) 3341 for (octave_idx_type i = cidx(j); i < cidx(j+1); i++)
3270 { 3342 {
3271 if (ridx(i) == j) 3343 if (ridx(i) == j)
3272 D[j] = data(i); 3344 D[j] = std::real (data(i));
3273 else if (ridx(i) == j + 1) 3345 else if (ridx(i) == j + 1)
3274 DL[j] = data(i); 3346 DL[j] = data(i);
3275 } 3347 }
3276 } 3348 }
3277 3349
3333 if (ridx(i) == j) 3405 if (ridx(i) == j)
3334 D[j] = data(i); 3406 D[j] = data(i);
3335 else if (ridx(i) == j + 1) 3407 else if (ridx(i) == j + 1)
3336 DL[j] = data(i); 3408 DL[j] = data(i);
3337 else if (ridx(i) == j - 1) 3409 else if (ridx(i) == j - 1)
3338 DU[j] = data(i); 3410 DU[j-1] = data(i);
3339 } 3411 }
3340 } 3412 }
3341 3413
3342 octave_idx_type b_nr = b.rows(); 3414 octave_idx_type b_nr = b.rows();
3343 octave_idx_type b_nc = b.cols(); 3415 octave_idx_type b_nc = b.cols();
3433 if (ridx(i) == j) 3505 if (ridx(i) == j)
3434 D[j] = data(i); 3506 D[j] = data(i);
3435 else if (ridx(i) == j + 1) 3507 else if (ridx(i) == j + 1)
3436 DL[j] = data(i); 3508 DL[j] = data(i);
3437 else if (ridx(i) == j - 1) 3509 else if (ridx(i) == j - 1)
3438 DU[j] = data(i); 3510 DU[j-1] = data(i);
3439 } 3511 }
3440 } 3512 }
3441 3513
3442 F77_XFCN (zgttrf, ZGTTRF, (nr, DL, D, DU, DU2, pipvt, err)); 3514 F77_XFCN (zgttrf, ZGTTRF, (nr, DL, D, DU, DU2, pipvt, err));
3443 3515
4455 4527
4456 #ifdef HAVE_UMFPACK 4528 #ifdef HAVE_UMFPACK
4457 // Setup the control parameters 4529 // Setup the control parameters
4458 Control = Matrix (UMFPACK_CONTROL, 1); 4530 Control = Matrix (UMFPACK_CONTROL, 1);
4459 double *control = Control.fortran_vec (); 4531 double *control = Control.fortran_vec ();
4460 umfpack_zi_defaults (control); 4532 UMFPACK_ZNAME (defaults) (control);
4461 4533
4462 double tmp = Voctave_sparse_controls.get_key ("spumoni"); 4534 double tmp = Voctave_sparse_controls.get_key ("spumoni");
4463 if (!xisnan (tmp)) 4535 if (!xisnan (tmp))
4464 Control (UMFPACK_PRL) = tmp; 4536 Control (UMFPACK_PRL) = tmp;
4465 tmp = Voctave_sparse_controls.get_key ("piv_tol"); 4537 tmp = Voctave_sparse_controls.get_key ("piv_tol");
4472 // Set whether we are allowed to modify Q or not 4544 // Set whether we are allowed to modify Q or not
4473 tmp = Voctave_sparse_controls.get_key ("autoamd"); 4545 tmp = Voctave_sparse_controls.get_key ("autoamd");
4474 if (!xisnan (tmp)) 4546 if (!xisnan (tmp))
4475 Control (UMFPACK_FIXQ) = tmp; 4547 Control (UMFPACK_FIXQ) = tmp;
4476 4548
4477 umfpack_zi_report_control (control); 4549 UMFPACK_ZNAME (report_control) (control);
4478 4550
4479 const octave_idx_type *Ap = cidx (); 4551 const octave_idx_type *Ap = cidx ();
4480 const octave_idx_type *Ai = ridx (); 4552 const octave_idx_type *Ai = ridx ();
4481 const Complex *Ax = data (); 4553 const Complex *Ax = data ();
4482 octave_idx_type nr = rows (); 4554 octave_idx_type nr = rows ();
4483 octave_idx_type nc = cols (); 4555 octave_idx_type nc = cols ();
4484 4556
4485 umfpack_zi_report_matrix (nr, nc, Ap, Ai, X_CAST (const double *, Ax), 4557 UMFPACK_ZNAME (report_matrix) (nr, nc, Ap, Ai,
4486 NULL, 1, control); 4558 X_CAST (const double *, Ax), NULL, 1, control);
4487 4559
4488 void *Symbolic; 4560 void *Symbolic;
4489 Info = Matrix (1, UMFPACK_INFO); 4561 Info = Matrix (1, UMFPACK_INFO);
4490 double *info = Info.fortran_vec (); 4562 double *info = Info.fortran_vec ();
4491 int status = umfpack_zi_qsymbolic (nr, nc, Ap, Ai, 4563 int status = UMFPACK_ZNAME (qsymbolic) (nr, nc, Ap, Ai,
4492 X_CAST (const double *, Ax), 4564 X_CAST (const double *, Ax),
4493 NULL, NULL, &Symbolic, control, info); 4565 NULL, NULL, &Symbolic, control, info);
4494 4566
4495 if (status < 0) 4567 if (status < 0)
4496 { 4568 {
4497 (*current_liboctave_error_handler) 4569 (*current_liboctave_error_handler)
4498 ("SparseComplexMatrix::solve symbolic factorization failed"); 4570 ("SparseComplexMatrix::solve symbolic factorization failed");
4499 err = -1; 4571 err = -1;
4500 4572
4501 umfpack_zi_report_status (control, status); 4573 UMFPACK_ZNAME (report_status) (control, status);
4502 umfpack_zi_report_info (control, info); 4574 UMFPACK_ZNAME (report_info) (control, info);
4503 4575
4504 umfpack_zi_free_symbolic (&Symbolic) ; 4576 UMFPACK_ZNAME (free_symbolic) (&Symbolic) ;
4505 } 4577 }
4506 else 4578 else
4507 { 4579 {
4508 umfpack_zi_report_symbolic (Symbolic, control); 4580 UMFPACK_ZNAME (report_symbolic) (Symbolic, control);
4509 4581
4510 status = umfpack_zi_numeric (Ap, Ai, X_CAST (const double *, Ax), NULL, 4582 status = UMFPACK_ZNAME (numeric) (Ap, Ai,
4583 X_CAST (const double *, Ax), NULL,
4511 Symbolic, &Numeric, control, info) ; 4584 Symbolic, &Numeric, control, info) ;
4512 umfpack_zi_free_symbolic (&Symbolic) ; 4585 UMFPACK_ZNAME (free_symbolic) (&Symbolic) ;
4513 4586
4514 #ifdef HAVE_LSSOLVE 4587 #ifdef HAVE_LSSOLVE
4515 rcond = Info (UMFPACK_RCOND); 4588 rcond = Info (UMFPACK_RCOND);
4516 volatile double rcond_plus_one = rcond + 1.0; 4589 volatile double rcond_plus_one = rcond + 1.0;
4517 4590
4518 if (status == UMFPACK_WARNING_singular_matrix || 4591 if (status == UMFPACK_WARNING_singular_matrix ||
4519 rcond_plus_one == 1.0 || xisnan (rcond)) 4592 rcond_plus_one == 1.0 || xisnan (rcond))
4520 { 4593 {
4521 umfpack_zi_report_numeric (Numeric, control); 4594 UMFPACK_ZNAME (report_numeric) (Numeric, control);
4522 4595
4523 err = -2; 4596 err = -2;
4524 4597
4525 if (sing_handler) 4598 if (sing_handler)
4526 sing_handler (rcond); 4599 sing_handler (rcond);
4535 if (status < 0) 4608 if (status < 0)
4536 { 4609 {
4537 (*current_liboctave_error_handler) 4610 (*current_liboctave_error_handler)
4538 ("SparseComplexMatrix::solve numeric factorization failed"); 4611 ("SparseComplexMatrix::solve numeric factorization failed");
4539 4612
4540 umfpack_zi_report_status (control, status); 4613 UMFPACK_ZNAME (report_status) (control, status);
4541 umfpack_zi_report_info (control, info); 4614 UMFPACK_ZNAME (report_info) (control, info);
4542 4615
4543 err = -1; 4616 err = -1;
4544 } 4617 }
4545 else 4618 else
4546 { 4619 {
4547 umfpack_zi_report_numeric (Numeric, control); 4620 UMFPACK_ZNAME (report_numeric) (Numeric, control);
4548 } 4621 }
4549 } 4622 }
4550 4623
4551 if (err != 0) 4624 if (err != 0)
4552 umfpack_zi_free_numeric (&Numeric); 4625 UMFPACK_ZNAME (free_numeric) (&Numeric);
4553 #else 4626 #else
4554 (*current_liboctave_error_handler) ("UMFPACK not installed"); 4627 (*current_liboctave_error_handler) ("UMFPACK not installed");
4555 #endif 4628 #endif
4556 4629
4557 return Numeric; 4630 return Numeric;
4618 Complex *Xx = retval.fortran_vec (); 4691 Complex *Xx = retval.fortran_vec ();
4619 4692
4620 for (octave_idx_type j = 0, iidx = 0; j < b_nc; j++, iidx += b_nr) 4693 for (octave_idx_type j = 0, iidx = 0; j < b_nc; j++, iidx += b_nr)
4621 { 4694 {
4622 #ifdef UMFPACK_SEPARATE_SPLIT 4695 #ifdef UMFPACK_SEPARATE_SPLIT
4623 status = umfpack_zi_solve (UMFPACK_A, Ap, Ai, 4696 status = UMFPACK_ZNAME (solve) (UMFPACK_A, Ap,
4624 X_CAST (const double *, Ax), 4697 Ai, X_CAST (const double *, Ax),
4625 NULL, 4698 NULL,
4626 X_CAST (double *, &Xx[iidx]), 4699 X_CAST (double *, &Xx[iidx]),
4627 NULL, 4700 NULL,
4628 &Bx[iidx], Bz, Numeric, 4701 &Bx[iidx], Bz, Numeric,
4629 control, info); 4702 control, info);
4630 #else 4703 #else
4631 for (octave_idx_type i = 0; i < b_nr; i++) 4704 for (octave_idx_type i = 0; i < b_nr; i++)
4632 Bz[i] = b.elem (i, j); 4705 Bz[i] = b.elem (i, j);
4633 4706
4634 status = umfpack_zi_solve (UMFPACK_A, Ap, Ai, 4707 status = UMFPACK_ZNAME (solve) (UMFPACK_A, Ap,
4635 X_CAST (const double *, Ax), 4708 Ai, X_CAST (const double *, Ax),
4636 NULL, 4709 NULL,
4637 X_CAST (double *, &Xx[iidx]), 4710 X_CAST (double *, &Xx[iidx]),
4638 NULL, 4711 NULL,
4639 X_CAST (const double *, Bz), 4712 X_CAST (const double *, Bz),
4640 NULL, Numeric, 4713 NULL, Numeric,
4644 if (status < 0) 4717 if (status < 0)
4645 { 4718 {
4646 (*current_liboctave_error_handler) 4719 (*current_liboctave_error_handler)
4647 ("SparseComplexMatrix::solve solve failed"); 4720 ("SparseComplexMatrix::solve solve failed");
4648 4721
4649 umfpack_zi_report_status (control, status); 4722 UMFPACK_ZNAME (report_status) (control, status);
4650 4723
4651 err = -1; 4724 err = -1;
4652 4725
4653 break; 4726 break;
4654 } 4727 }
4671 rcond); 4744 rcond);
4672 4745
4673 } 4746 }
4674 #endif 4747 #endif
4675 4748
4676 umfpack_zi_report_info (control, info); 4749 UMFPACK_ZNAME (report_info) (control, info);
4677 4750
4678 umfpack_zi_free_numeric (&Numeric); 4751 UMFPACK_ZNAME (free_numeric) (&Numeric);
4679 } 4752 }
4680 #else 4753 #else
4681 (*current_liboctave_error_handler) ("UMFPACK not installed"); 4754 (*current_liboctave_error_handler) ("UMFPACK not installed");
4682 #endif 4755 #endif
4683 } 4756 }
4760 4833
4761 #ifdef UMFPACK_SEPARATE_SPLIT 4834 #ifdef UMFPACK_SEPARATE_SPLIT
4762 for (octave_idx_type i = 0; i < b_nr; i++) 4835 for (octave_idx_type i = 0; i < b_nr; i++)
4763 Bx[i] = b.elem (i, j); 4836 Bx[i] = b.elem (i, j);
4764 4837
4765 status = umfpack_zi_solve (UMFPACK_A, Ap, Ai, 4838 status = UMFPACK_ZNAME (solve) (UMFPACK_A, Ap,
4766 X_CAST (const double *, Ax), 4839 Ai, X_CAST (const double *, Ax),
4767 NULL, 4840 NULL,
4768 X_CAST (double *, Xx), NULL, 4841 X_CAST (double *, Xx), NULL,
4769 Bx, Bz, Numeric, control, 4842 Bx, Bz, Numeric, control,
4770 info); 4843 info);
4771 #else 4844 #else
4772 for (octave_idx_type i = 0; i < b_nr; i++) 4845 for (octave_idx_type i = 0; i < b_nr; i++)
4773 Bz[i] = b.elem (i, j); 4846 Bz[i] = b.elem (i, j);
4774 4847
4775 status = umfpack_zi_solve (UMFPACK_A, Ap, Ai, 4848 status = UMFPACK_ZNAME (solve) (UMFPACK_A, Ap, Ai,
4776 X_CAST (const double *, Ax), 4849 X_CAST (const double *, Ax),
4777 NULL, 4850 NULL,
4778 X_CAST (double *, Xx), NULL, 4851 X_CAST (double *, Xx), NULL,
4779 X_CAST (double *, Bz), NULL, 4852 X_CAST (double *, Bz), NULL,
4780 Numeric, control, 4853 Numeric, control,
4783 if (status < 0) 4856 if (status < 0)
4784 { 4857 {
4785 (*current_liboctave_error_handler) 4858 (*current_liboctave_error_handler)
4786 ("SparseComplexMatrix::solve solve failed"); 4859 ("SparseComplexMatrix::solve solve failed");
4787 4860
4788 umfpack_zi_report_status (control, status); 4861 UMFPACK_ZNAME (report_status) (control, status);
4789 4862
4790 err = -1; 4863 err = -1;
4791 4864
4792 break; 4865 break;
4793 } 4866 }
4831 rcond); 4904 rcond);
4832 4905
4833 } 4906 }
4834 #endif 4907 #endif
4835 4908
4836 umfpack_zi_report_info (control, info); 4909 UMFPACK_ZNAME (report_info) (control, info);
4837 4910
4838 umfpack_zi_free_numeric (&Numeric); 4911 UMFPACK_ZNAME (free_numeric) (&Numeric);
4839 } 4912 }
4840 #else 4913 #else
4841 (*current_liboctave_error_handler) ("UMFPACK not installed"); 4914 (*current_liboctave_error_handler) ("UMFPACK not installed");
4842 #endif 4915 #endif
4843 } 4916 }
4902 Complex *Xx = retval.fortran_vec (); 4975 Complex *Xx = retval.fortran_vec ();
4903 4976
4904 for (octave_idx_type j = 0, iidx = 0; j < b_nc; j++, iidx += b_nr) 4977 for (octave_idx_type j = 0, iidx = 0; j < b_nc; j++, iidx += b_nr)
4905 { 4978 {
4906 status = 4979 status =
4907 umfpack_zi_solve (UMFPACK_A, Ap, Ai, 4980 UMFPACK_ZNAME (solve) (UMFPACK_A, Ap, Ai,
4908 X_CAST (const double *, Ax), 4981 X_CAST (const double *, Ax),
4909 NULL, X_CAST (double *, &Xx[iidx]), 4982 NULL, X_CAST (double *, &Xx[iidx]),
4910 NULL, X_CAST (const double *, &Bx[iidx]), 4983 NULL, X_CAST (const double *, &Bx[iidx]),
4911 NULL, Numeric, control, info); 4984 NULL, Numeric, control, info);
4912 4985
4913 if (status < 0) 4986 if (status < 0)
4914 { 4987 {
4915 (*current_liboctave_error_handler) 4988 (*current_liboctave_error_handler)
4916 ("SparseComplexMatrix::solve solve failed"); 4989 ("SparseComplexMatrix::solve solve failed");
4917 4990
4918 umfpack_zi_report_status (control, status); 4991 UMFPACK_ZNAME (report_status) (control, status);
4919 4992
4920 err = -1; 4993 err = -1;
4921 4994
4922 break; 4995 break;
4923 } 4996 }
4940 rcond); 5013 rcond);
4941 5014
4942 } 5015 }
4943 #endif 5016 #endif
4944 5017
4945 umfpack_zi_report_info (control, info); 5018 UMFPACK_ZNAME (report_info) (control, info);
4946 5019
4947 umfpack_zi_free_numeric (&Numeric); 5020 UMFPACK_ZNAME (free_numeric) (&Numeric);
4948 } 5021 }
4949 #else 5022 #else
4950 (*current_liboctave_error_handler) ("UMFPACK not installed"); 5023 (*current_liboctave_error_handler) ("UMFPACK not installed");
4951 #endif 5024 #endif
4952 } 5025 }
5020 for (octave_idx_type j = 0; j < b_nc; j++) 5093 for (octave_idx_type j = 0; j < b_nc; j++)
5021 { 5094 {
5022 for (octave_idx_type i = 0; i < b_nr; i++) 5095 for (octave_idx_type i = 0; i < b_nr; i++)
5023 Bx[i] = b (i,j); 5096 Bx[i] = b (i,j);
5024 5097
5025 status = umfpack_zi_solve (UMFPACK_A, Ap, Ai, 5098 status = UMFPACK_ZNAME (solve) (UMFPACK_A, Ap,
5026 X_CAST (const double *, Ax), 5099 Ai, X_CAST (const double *, Ax),
5027 NULL, X_CAST (double *, Xx), 5100 NULL, X_CAST (double *, Xx),
5028 NULL, X_CAST (double *, Bx), 5101 NULL, X_CAST (double *, Bx),
5029 NULL, Numeric, control, info); 5102 NULL, Numeric, control, info);
5030 5103
5031 if (status < 0) 5104 if (status < 0)
5032 { 5105 {
5033 (*current_liboctave_error_handler) 5106 (*current_liboctave_error_handler)
5034 ("SparseComplexMatrix::solve solve failed"); 5107 ("SparseComplexMatrix::solve solve failed");
5035 5108
5036 umfpack_zi_report_status (control, status); 5109 UMFPACK_ZNAME (report_status) (control, status);
5037 5110
5038 err = -1; 5111 err = -1;
5039 5112
5040 break; 5113 break;
5041 } 5114 }
5079 rcond); 5152 rcond);
5080 5153
5081 } 5154 }
5082 #endif 5155 #endif
5083 5156
5084 umfpack_zi_report_info (control, info); 5157 UMFPACK_ZNAME (report_info) (control, info);
5085 5158
5086 umfpack_zi_free_numeric (&Numeric); 5159 UMFPACK_ZNAME (free_numeric) (&Numeric);
5087 } 5160 }
5088 #else 5161 #else
5089 (*current_liboctave_error_handler) ("UMFPACK not installed"); 5162 (*current_liboctave_error_handler) ("UMFPACK not installed");
5090 #endif 5163 #endif
5091 } 5164 }
5122 ComplexMatrix 5195 ComplexMatrix
5123 SparseComplexMatrix::solve (SparseType &mattype, const Matrix& b, octave_idx_type& err, 5196 SparseComplexMatrix::solve (SparseType &mattype, const Matrix& b, octave_idx_type& err,
5124 double& rcond, 5197 double& rcond,
5125 solve_singularity_handler sing_handler) const 5198 solve_singularity_handler sing_handler) const
5126 { 5199 {
5127 int typ = mattype.type (); 5200 int typ = mattype.type (false);
5128 5201
5129 if (typ == SparseType::Unknown) 5202 if (typ == SparseType::Unknown)
5130 typ = mattype.type (*this); 5203 typ = mattype.type (*this);
5131 5204
5132 if (typ == SparseType::Diagonal || typ == SparseType::Permuted_Diagonal) 5205 if (typ == SparseType::Diagonal || typ == SparseType::Permuted_Diagonal)
5176 SparseComplexMatrix 5249 SparseComplexMatrix
5177 SparseComplexMatrix::solve (SparseType &mattype, const SparseMatrix& b, 5250 SparseComplexMatrix::solve (SparseType &mattype, const SparseMatrix& b,
5178 octave_idx_type& err, double& rcond, 5251 octave_idx_type& err, double& rcond,
5179 solve_singularity_handler sing_handler) const 5252 solve_singularity_handler sing_handler) const
5180 { 5253 {
5181 int typ = mattype.type (); 5254 int typ = mattype.type (false);
5182 5255
5183 if (typ == SparseType::Unknown) 5256 if (typ == SparseType::Unknown)
5184 typ = mattype.type (*this); 5257 typ = mattype.type (*this);
5185 5258
5186 if (typ == SparseType::Diagonal || typ == SparseType::Permuted_Diagonal) 5259 if (typ == SparseType::Diagonal || typ == SparseType::Permuted_Diagonal)
5230 ComplexMatrix 5303 ComplexMatrix
5231 SparseComplexMatrix::solve (SparseType &mattype, const ComplexMatrix& b, 5304 SparseComplexMatrix::solve (SparseType &mattype, const ComplexMatrix& b,
5232 octave_idx_type& err, double& rcond, 5305 octave_idx_type& err, double& rcond,
5233 solve_singularity_handler sing_handler) const 5306 solve_singularity_handler sing_handler) const
5234 { 5307 {
5235 int typ = mattype.type (); 5308 int typ = mattype.type (false);
5236 5309
5237 if (typ == SparseType::Unknown) 5310 if (typ == SparseType::Unknown)
5238 typ = mattype.type (*this); 5311 typ = mattype.type (*this);
5239 5312
5240 if (typ == SparseType::Diagonal || typ == SparseType::Permuted_Diagonal) 5313 if (typ == SparseType::Diagonal || typ == SparseType::Permuted_Diagonal)
5285 SparseComplexMatrix 5358 SparseComplexMatrix
5286 SparseComplexMatrix::solve (SparseType &mattype, const SparseComplexMatrix& b, 5359 SparseComplexMatrix::solve (SparseType &mattype, const SparseComplexMatrix& b,
5287 octave_idx_type& err, double& rcond, 5360 octave_idx_type& err, double& rcond,
5288 solve_singularity_handler sing_handler) const 5361 solve_singularity_handler sing_handler) const
5289 { 5362 {
5290 int typ = mattype.type (); 5363 int typ = mattype.type (false);
5291 5364
5292 if (typ == SparseType::Unknown) 5365 if (typ == SparseType::Unknown)
5293 typ = mattype.type (*this); 5366 typ = mattype.type (*this);
5294 5367
5295 if (typ == SparseType::Diagonal || typ == SparseType::Permuted_Diagonal) 5368 if (typ == SparseType::Diagonal || typ == SparseType::Permuted_Diagonal)