Mercurial > hg > octave-nkf

comparison liboctave/dMatrix.cc @ 11570:57632dea2446

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
attempt better backward compatibility for Array constructors
author John W. Eaton <jwe@octave.org>
date Wed, 19 Jan 2011 17:55:56 -0500
parents fd0a3ac60b0e
children a83bad07f7e3
comparison
equal deleted inserted replaced
11569:7e9a111cae20 11570:57632dea2446
248 : MArray<double> (cv) 248 : MArray<double> (cv)
249 { 249 {
250 } 250 }
251 251
252 Matrix::Matrix (const DiagMatrix& a) 252 Matrix::Matrix (const DiagMatrix& a)
253 : MArray<double> (a.rows (), a.cols (), 0.0) 253 : MArray<double> (a.dims (), 0.0)
254 { 254 {
255 for (octave_idx_type i = 0; i < a.length (); i++) 255 for (octave_idx_type i = 0; i < a.length (); i++)
256 elem (i, i) = a.elem (i, i); 256 elem (i, i) = a.elem (i, i);
257 } 257 }
258 258
259 Matrix::Matrix (const PermMatrix& a) 259 Matrix::Matrix (const PermMatrix& a)
260 : MArray<double> (a.rows (), a.cols (), 0.0) 260 : MArray<double> (a.dims (), 0.0)
261 { 261 {
262 const Array<octave_idx_type> ia (a.pvec ()); 262 const Array<octave_idx_type> ia (a.pvec ());
263 octave_idx_type len = a.rows (); 263 octave_idx_type len = a.rows ();
264 if (a.is_col_perm ()) 264 if (a.is_col_perm ())
265 for (octave_idx_type i = 0; i < len; i++) 265 for (octave_idx_type i = 0; i < len; i++)
276 : MArray<double> (a) 276 : MArray<double> (a)
277 { 277 {
278 } 278 }
279 279
280 Matrix::Matrix (const charMatrix& a) 280 Matrix::Matrix (const charMatrix& a)
281 : MArray<double> (a.rows (), a.cols ()) 281 : MArray<double> (a.dims ())
282 { 282 {
283 for (octave_idx_type i = 0; i < a.rows (); i++) 283 for (octave_idx_type i = 0; i < a.rows (); i++)
284 for (octave_idx_type j = 0; j < a.cols (); j++) 284 for (octave_idx_type j = 0; j < a.cols (); j++)
285 elem (i, j) = static_cast<unsigned char> (a.elem (i, j)); 285 elem (i, j) = static_cast<unsigned char> (a.elem (i, j));
286 } 286 }
738 738
739 if (nr != nc || nr == 0 || nc == 0) 739 if (nr != nc || nr == 0 || nc == 0)
740 (*current_liboctave_error_handler) ("inverse requires square matrix"); 740 (*current_liboctave_error_handler) ("inverse requires square matrix");
741 else 741 else
742 { 742 {
743 Array<octave_idx_type> ipvt (nr, 1); 743 Array<octave_idx_type> ipvt (dim_vector (nr, 1));
744 octave_idx_type *pipvt = ipvt.fortran_vec (); 744 octave_idx_type *pipvt = ipvt.fortran_vec ();
745 745
746 retval = *this; 746 retval = *this;
747 double *tmp_data = retval.fortran_vec (); 747 double *tmp_data = retval.fortran_vec ();
748 748
749 Array<double> z(1, 1); 749 Array<double> z (dim_vector (1, 1));
750 octave_idx_type lwork = -1; 750 octave_idx_type lwork = -1;
751 751
752 // Query the optimum work array size. 752 // Query the optimum work array size.
753 F77_XFCN (dgetri, DGETRI, (nc, tmp_data, nr, pipvt, 753 F77_XFCN (dgetri, DGETRI, (nc, tmp_data, nr, pipvt,
754 z.fortran_vec (), lwork, info)); 754 z.fortran_vec (), lwork, info));
775 { 775 {
776 octave_idx_type dgecon_info = 0; 776 octave_idx_type dgecon_info = 0;
777 777
778 // Now calculate the condition number for non-singular matrix. 778 // Now calculate the condition number for non-singular matrix.
779 char job = '1'; 779 char job = '1';
780 Array<octave_idx_type> iz (nc, 1); 780 Array<octave_idx_type> iz (dim_vector (nc, 1));
781 octave_idx_type *piz = iz.fortran_vec (); 781 octave_idx_type *piz = iz.fortran_vec ();
782 F77_XFCN (dgecon, DGECON, (F77_CONST_CHAR_ARG2 (&job, 1), 782 F77_XFCN (dgecon, DGECON, (F77_CONST_CHAR_ARG2 (&job, 1),
783 nc, tmp_data, nr, anorm, 783 nc, tmp_data, nr, anorm,
784 rcon, pz, piz, dgecon_info 784 rcon, pz, piz, dgecon_info
785 F77_CHAR_ARG_LEN (1))); 785 F77_CHAR_ARG_LEN (1)));
1010 nsamples = nc; 1010 nsamples = nc;
1011 } 1011 }
1012 1012
1013 octave_idx_type nn = 4*npts+15; 1013 octave_idx_type nn = 4*npts+15;
1014 1014
1015 Array<Complex> wsave (nn, 1); 1015 Array<Complex> wsave (dim_vector (nn, 1));
1016 Complex *pwsave = wsave.fortran_vec (); 1016 Complex *pwsave = wsave.fortran_vec ();
1017 1017
1018 retval = ComplexMatrix (*this); 1018 retval = ComplexMatrix (*this);
1019 Complex *tmp_data = retval.fortran_vec (); 1019 Complex *tmp_data = retval.fortran_vec ();
1020 1020
1051 nsamples = nc; 1051 nsamples = nc;
1052 } 1052 }
1053 1053
1054 octave_idx_type nn = 4*npts+15; 1054 octave_idx_type nn = 4*npts+15;
1055 1055
1056 Array<Complex> wsave (nn, 1); 1056 Array<Complex> wsave (dim_vector (nn, 1));
1057 Complex *pwsave = wsave.fortran_vec (); 1057 Complex *pwsave = wsave.fortran_vec ();
1058 1058
1059 retval = ComplexMatrix (*this); 1059 retval = ComplexMatrix (*this);
1060 Complex *tmp_data = retval.fortran_vec (); 1060 Complex *tmp_data = retval.fortran_vec ();
1061 1061
1095 nsamples = nc; 1095 nsamples = nc;
1096 } 1096 }
1097 1097
1098 octave_idx_type nn = 4*npts+15; 1098 octave_idx_type nn = 4*npts+15;
1099 1099
1100 Array<Complex> wsave (nn, 1); 1100 Array<Complex> wsave (dim_vector (nn, 1));
1101 Complex *pwsave = wsave.fortran_vec (); 1101 Complex *pwsave = wsave.fortran_vec ();
1102 1102
1103 retval = ComplexMatrix (*this); 1103 retval = ComplexMatrix (*this);
1104 Complex *tmp_data = retval.fortran_vec (); 1104 Complex *tmp_data = retval.fortran_vec ();
1105 1105
1117 nn = 4*npts+15; 1117 nn = 4*npts+15;
1118 1118
1119 wsave.resize (nn, 1); 1119 wsave.resize (nn, 1);
1120 pwsave = wsave.fortran_vec (); 1120 pwsave = wsave.fortran_vec ();
1121 1121
1122 Array<Complex> tmp (npts, 1); 1122 Array<Complex> tmp (dim_vector (npts, 1));
1123 Complex *prow = tmp.fortran_vec (); 1123 Complex *prow = tmp.fortran_vec ();
1124 1124
1125 F77_FUNC (zffti, ZFFTI) (npts, pwsave); 1125 F77_FUNC (zffti, ZFFTI) (npts, pwsave);
1126 1126
1127 for (octave_idx_type j = 0; j < nsamples; j++) 1127 for (octave_idx_type j = 0; j < nsamples; j++)
1161 nsamples = nc; 1161 nsamples = nc;
1162 } 1162 }
1163 1163
1164 octave_idx_type nn = 4*npts+15; 1164 octave_idx_type nn = 4*npts+15;
1165 1165
1166 Array<Complex> wsave (nn, 1); 1166 Array<Complex> wsave (dim_vector (nn, 1));
1167 Complex *pwsave = wsave.fortran_vec (); 1167 Complex *pwsave = wsave.fortran_vec ();
1168 1168
1169 retval = ComplexMatrix (*this); 1169 retval = ComplexMatrix (*this);
1170 Complex *tmp_data = retval.fortran_vec (); 1170 Complex *tmp_data = retval.fortran_vec ();
1171 1171
1186 nn = 4*npts+15; 1186 nn = 4*npts+15;
1187 1187
1188 wsave.resize (nn, 1); 1188 wsave.resize (nn, 1);
1189 pwsave = wsave.fortran_vec (); 1189 pwsave = wsave.fortran_vec ();
1190 1190
1191 Array<Complex> tmp (npts, 1); 1191 Array<Complex> tmp (dim_vector (npts, 1));
1192 Complex *prow = tmp.fortran_vec (); 1192 Complex *prow = tmp.fortran_vec ();
1193 1193
1194 F77_FUNC (zffti, ZFFTI) (npts, pwsave); 1194 F77_FUNC (zffti, ZFFTI) (npts, pwsave);
1195 1195
1196 for (octave_idx_type j = 0; j < nsamples; j++) 1196 for (octave_idx_type j = 0; j < nsamples; j++)
1283 mattype.mark_as_unsymmetric (); 1283 mattype.mark_as_unsymmetric ();
1284 typ = MatrixType::Full; 1284 typ = MatrixType::Full;
1285 } 1285 }
1286 else 1286 else
1287 { 1287 {
1288 Array<double> z (3 * nc, 1); 1288 Array<double> z (dim_vector (3 * nc, 1));
1289 double *pz = z.fortran_vec (); 1289 double *pz = z.fortran_vec ();
1290 Array<octave_idx_type> iz (nc, 1); 1290 Array<octave_idx_type> iz (dim_vector (nc, 1));
1291 octave_idx_type *piz = iz.fortran_vec (); 1291 octave_idx_type *piz = iz.fortran_vec ();
1292 1292
1293 F77_XFCN (dpocon, DPOCON, (F77_CONST_CHAR_ARG2 (&job, 1), 1293 F77_XFCN (dpocon, DPOCON, (F77_CONST_CHAR_ARG2 (&job, 1),
1294 nr, tmp_data, nr, anorm, 1294 nr, tmp_data, nr, anorm,
1295 rcon, pz, piz, info 1295 rcon, pz, piz, info
1307 else if (typ != MatrixType::Full) 1307 else if (typ != MatrixType::Full)
1308 (*current_liboctave_error_handler) ("det: invalid dense matrix type"); 1308 (*current_liboctave_error_handler) ("det: invalid dense matrix type");
1309 1309
1310 if (typ == MatrixType::Full) 1310 if (typ == MatrixType::Full)
1311 { 1311 {
1312 Array<octave_idx_type> ipvt (nr, 1); 1312 Array<octave_idx_type> ipvt (dim_vector (nr, 1));
1313 octave_idx_type *pipvt = ipvt.fortran_vec (); 1313 octave_idx_type *pipvt = ipvt.fortran_vec ();
1314 1314
1315 Matrix atmp = *this; 1315 Matrix atmp = *this;
1316 double *tmp_data = atmp.fortran_vec (); 1316 double *tmp_data = atmp.fortran_vec ();
1317 1317
1334 { 1334 {
1335 if (calc_cond) 1335 if (calc_cond)
1336 { 1336 {
1337 // Now calc the condition number for non-singular matrix. 1337 // Now calc the condition number for non-singular matrix.
1338 char job = '1'; 1338 char job = '1';
1339 Array<double> z (4 * nc, 1); 1339 Array<double> z (dim_vector (4 * nc, 1));
1340 double *pz = z.fortran_vec (); 1340 double *pz = z.fortran_vec ();
1341 Array<octave_idx_type> iz (nc, 1); 1341 Array<octave_idx_type> iz (dim_vector (nc, 1));
1342 octave_idx_type *piz = iz.fortran_vec (); 1342 octave_idx_type *piz = iz.fortran_vec ();
1343 1343
1344 F77_XFCN (dgecon, DGECON, (F77_CONST_CHAR_ARG2 (&job, 1), 1344 F77_XFCN (dgecon, DGECON, (F77_CONST_CHAR_ARG2 (&job, 1),
1345 nc, tmp_data, nr, anorm, 1345 nc, tmp_data, nr, anorm,
1346 rcon, pz, piz, info 1346 rcon, pz, piz, info
1399 octave_idx_type info = 0; 1399 octave_idx_type info = 0;
1400 char norm = '1'; 1400 char norm = '1';
1401 char uplo = 'U'; 1401 char uplo = 'U';
1402 char dia = 'N'; 1402 char dia = 'N';
1403 1403
1404 Array<double> z (3 * nc, 1); 1404 Array<double> z (dim_vector (3 * nc, 1));
1405 double *pz = z.fortran_vec (); 1405 double *pz = z.fortran_vec ();
1406 Array<octave_idx_type> iz (nc, 1); 1406 Array<octave_idx_type> iz (dim_vector (nc, 1));
1407 octave_idx_type *piz = iz.fortran_vec (); 1407 octave_idx_type *piz = iz.fortran_vec ();
1408 1408
1409 F77_XFCN (dtrcon, DTRCON, (F77_CONST_CHAR_ARG2 (&norm, 1), 1409 F77_XFCN (dtrcon, DTRCON, (F77_CONST_CHAR_ARG2 (&norm, 1),
1410 F77_CONST_CHAR_ARG2 (&uplo, 1), 1410 F77_CONST_CHAR_ARG2 (&uplo, 1),
1411 F77_CONST_CHAR_ARG2 (&dia, 1), 1411 F77_CONST_CHAR_ARG2 (&dia, 1),
1427 octave_idx_type info = 0; 1427 octave_idx_type info = 0;
1428 char norm = '1'; 1428 char norm = '1';
1429 char uplo = 'L'; 1429 char uplo = 'L';
1430 char dia = 'N'; 1430 char dia = 'N';
1431 1431
1432 Array<double> z (3 * nc, 1); 1432 Array<double> z (dim_vector (3 * nc, 1));
1433 double *pz = z.fortran_vec (); 1433 double *pz = z.fortran_vec ();
1434 Array<octave_idx_type> iz (nc, 1); 1434 Array<octave_idx_type> iz (dim_vector (nc, 1));
1435 octave_idx_type *piz = iz.fortran_vec (); 1435 octave_idx_type *piz = iz.fortran_vec ();
1436 1436
1437 F77_XFCN (dtrcon, DTRCON, (F77_CONST_CHAR_ARG2 (&norm, 1), 1437 F77_XFCN (dtrcon, DTRCON, (F77_CONST_CHAR_ARG2 (&norm, 1),
1438 F77_CONST_CHAR_ARG2 (&uplo, 1), 1438 F77_CONST_CHAR_ARG2 (&uplo, 1),
1439 F77_CONST_CHAR_ARG2 (&dia, 1), 1439 F77_CONST_CHAR_ARG2 (&dia, 1),
1472 mattype.mark_as_unsymmetric (); 1472 mattype.mark_as_unsymmetric ();
1473 typ = MatrixType::Full; 1473 typ = MatrixType::Full;
1474 } 1474 }
1475 else 1475 else
1476 { 1476 {
1477 Array<double> z (3 * nc, 1); 1477 Array<double> z (dim_vector (3 * nc, 1));
1478 double *pz = z.fortran_vec (); 1478 double *pz = z.fortran_vec ();
1479 Array<octave_idx_type> iz (nc, 1); 1479 Array<octave_idx_type> iz (dim_vector (nc, 1));
1480 octave_idx_type *piz = iz.fortran_vec (); 1480 octave_idx_type *piz = iz.fortran_vec ();
1481 1481
1482 F77_XFCN (dpocon, DPOCON, (F77_CONST_CHAR_ARG2 (&job, 1), 1482 F77_XFCN (dpocon, DPOCON, (F77_CONST_CHAR_ARG2 (&job, 1),
1483 nr, tmp_data, nr, anorm, 1483 nr, tmp_data, nr, anorm,
1484 rcon, pz, piz, info 1484 rcon, pz, piz, info
1491 1491
1492 if (typ == MatrixType::Full) 1492 if (typ == MatrixType::Full)
1493 { 1493 {
1494 octave_idx_type info = 0; 1494 octave_idx_type info = 0;
1495 1495
1496 Array<octave_idx_type> ipvt (nr, 1); 1496 Array<octave_idx_type> ipvt (dim_vector (nr, 1));
1497 octave_idx_type *pipvt = ipvt.fortran_vec (); 1497 octave_idx_type *pipvt = ipvt.fortran_vec ();
1498 1498
1499 if(anorm < 0.) 1499 if(anorm < 0.)
1500 anorm = atmp.abs().sum(). 1500 anorm = atmp.abs().sum().
1501 row(static_cast<octave_idx_type>(0)).max(); 1501 row(static_cast<octave_idx_type>(0)).max();
1502 1502
1503 Array<double> z (4 * nc, 1); 1503 Array<double> z (dim_vector (4 * nc, 1));
1504 double *pz = z.fortran_vec (); 1504 double *pz = z.fortran_vec ();
1505 Array<octave_idx_type> iz (nc, 1); 1505 Array<octave_idx_type> iz (dim_vector (nc, 1));
1506 octave_idx_type *piz = iz.fortran_vec (); 1506 octave_idx_type *piz = iz.fortran_vec ();
1507 1507
1508 F77_XFCN (dgetrf, DGETRF, (nr, nr, tmp_data, nr, pipvt, info)); 1508 F77_XFCN (dgetrf, DGETRF, (nr, nr, tmp_data, nr, pipvt, info));
1509 1509
1510 if (info != 0) 1510 if (info != 0)
1571 { 1571 {
1572 char norm = '1'; 1572 char norm = '1';
1573 char uplo = 'U'; 1573 char uplo = 'U';
1574 char dia = 'N'; 1574 char dia = 'N';
1575 1575
1576 Array<double> z (3 * nc, 1); 1576 Array<double> z (dim_vector (3 * nc, 1));
1577 double *pz = z.fortran_vec (); 1577 double *pz = z.fortran_vec ();
1578 Array<octave_idx_type> iz (nc, 1); 1578 Array<octave_idx_type> iz (dim_vector (nc, 1));
1579 octave_idx_type *piz = iz.fortran_vec (); 1579 octave_idx_type *piz = iz.fortran_vec ();
1580 1580
1581 F77_XFCN (dtrcon, DTRCON, (F77_CONST_CHAR_ARG2 (&norm, 1), 1581 F77_XFCN (dtrcon, DTRCON, (F77_CONST_CHAR_ARG2 (&norm, 1),
1582 F77_CONST_CHAR_ARG2 (&uplo, 1), 1582 F77_CONST_CHAR_ARG2 (&uplo, 1),
1583 F77_CONST_CHAR_ARG2 (&dia, 1), 1583 F77_CONST_CHAR_ARG2 (&dia, 1),
1671 { 1671 {
1672 char norm = '1'; 1672 char norm = '1';
1673 char uplo = 'L'; 1673 char uplo = 'L';
1674 char dia = 'N'; 1674 char dia = 'N';
1675 1675
1676 Array<double> z (3 * nc, 1); 1676 Array<double> z (dim_vector (3 * nc, 1));
1677 double *pz = z.fortran_vec (); 1677 double *pz = z.fortran_vec ();
1678 Array<octave_idx_type> iz (nc, 1); 1678 Array<octave_idx_type> iz (dim_vector (nc, 1));
1679 octave_idx_type *piz = iz.fortran_vec (); 1679 octave_idx_type *piz = iz.fortran_vec ();
1680 1680
1681 F77_XFCN (dtrcon, DTRCON, (F77_CONST_CHAR_ARG2 (&norm, 1), 1681 F77_XFCN (dtrcon, DTRCON, (F77_CONST_CHAR_ARG2 (&norm, 1),
1682 F77_CONST_CHAR_ARG2 (&uplo, 1), 1682 F77_CONST_CHAR_ARG2 (&uplo, 1),
1683 F77_CONST_CHAR_ARG2 (&dia, 1), 1683 F77_CONST_CHAR_ARG2 (&dia, 1),
1777 } 1777 }
1778 else 1778 else
1779 { 1779 {
1780 if (calc_cond) 1780 if (calc_cond)
1781 { 1781 {
1782 Array<double> z (3 * nc, 1); 1782 Array<double> z (dim_vector (3 * nc, 1));
1783 double *pz = z.fortran_vec (); 1783 double *pz = z.fortran_vec ();
1784 Array<octave_idx_type> iz (nc, 1); 1784 Array<octave_idx_type> iz (dim_vector (nc, 1));
1785 octave_idx_type *piz = iz.fortran_vec (); 1785 octave_idx_type *piz = iz.fortran_vec ();
1786 1786
1787 F77_XFCN (dpocon, DPOCON, (F77_CONST_CHAR_ARG2 (&job, 1), 1787 F77_XFCN (dpocon, DPOCON, (F77_CONST_CHAR_ARG2 (&job, 1),
1788 nr, tmp_data, nr, anorm, 1788 nr, tmp_data, nr, anorm,
1789 rcon, pz, piz, info 1789 rcon, pz, piz, info
1829 1829
1830 if (typ == MatrixType::Full) 1830 if (typ == MatrixType::Full)
1831 { 1831 {
1832 info = 0; 1832 info = 0;
1833 1833
1834 Array<octave_idx_type> ipvt (nr, 1); 1834 Array<octave_idx_type> ipvt (dim_vector (nr, 1));
1835 octave_idx_type *pipvt = ipvt.fortran_vec (); 1835 octave_idx_type *pipvt = ipvt.fortran_vec ();
1836 1836
1837 Matrix atmp = *this; 1837 Matrix atmp = *this;
1838 double *tmp_data = atmp.fortran_vec (); 1838 double *tmp_data = atmp.fortran_vec ();
1839 if(anorm < 0.) 1839 if(anorm < 0.)
1840 anorm = atmp.abs().sum().row(static_cast<octave_idx_type>(0)).max(); 1840 anorm = atmp.abs().sum().row(static_cast<octave_idx_type>(0)).max();
1841 1841
1842 Array<double> z (4 * nc, 1); 1842 Array<double> z (dim_vector (4 * nc, 1));
1843 double *pz = z.fortran_vec (); 1843 double *pz = z.fortran_vec ();
1844 Array<octave_idx_type> iz (nc, 1); 1844 Array<octave_idx_type> iz (dim_vector (nc, 1));
1845 octave_idx_type *piz = iz.fortran_vec (); 1845 octave_idx_type *piz = iz.fortran_vec ();
1846 1846
1847 F77_XFCN (dgetrf, DGETRF, (nr, nr, tmp_data, nr, pipvt, info)); 1847 F77_XFCN (dgetrf, DGETRF, (nr, nr, tmp_data, nr, pipvt, info));
1848 1848
1849 // Throw-away extra info LAPACK gives so as to not change output. 1849 // Throw-away extra info LAPACK gives so as to not change output.
2269 2269
2270 Matrix atmp = *this; 2270 Matrix atmp = *this;
2271 double *tmp_data = atmp.fortran_vec (); 2271 double *tmp_data = atmp.fortran_vec ();
2272 2272
2273 double *pretval = retval.fortran_vec (); 2273 double *pretval = retval.fortran_vec ();
2274 Array<double> s (minmn, 1); 2274 Array<double> s (dim_vector (minmn, 1));
2275 double *ps = s.fortran_vec (); 2275 double *ps = s.fortran_vec ();
2276 2276
2277 // Ask DGELSD what the dimension of WORK should be. 2277 // Ask DGELSD what the dimension of WORK should be.
2278 octave_idx_type lwork = -1; 2278 octave_idx_type lwork = -1;
2279 2279
2280 Array<double> work (1, 1); 2280 Array<double> work (dim_vector (1, 1));
2281 2281
2282 octave_idx_type smlsiz; 2282 octave_idx_type smlsiz;
2283 F77_FUNC (xilaenv, XILAENV) (9, F77_CONST_CHAR_ARG2 ("DGELSD", 6), 2283 F77_FUNC (xilaenv, XILAENV) (9, F77_CONST_CHAR_ARG2 ("DGELSD", 6),
2284 F77_CONST_CHAR_ARG2 (" ", 1), 2284 F77_CONST_CHAR_ARG2 (" ", 1),
2285 0, 0, 0, 0, smlsiz 2285 0, 0, 0, 0, smlsiz
2307 nlvl = 0; 2307 nlvl = 0;
2308 2308
2309 octave_idx_type liwork = 3 * minmn * nlvl + 11 * minmn; 2309 octave_idx_type liwork = 3 * minmn * nlvl + 11 * minmn;
2310 if (liwork < 1) 2310 if (liwork < 1)
2311 liwork = 1; 2311 liwork = 1;
2312 Array<octave_idx_type> iwork (liwork, 1); 2312 Array<octave_idx_type> iwork (dim_vector (liwork, 1));
2313 octave_idx_type* piwork = iwork.fortran_vec (); 2313 octave_idx_type* piwork = iwork.fortran_vec ();
2314 2314
2315 F77_XFCN (dgelsd, DGELSD, (m, n, nrhs, tmp_data, m, pretval, maxmn, 2315 F77_XFCN (dgelsd, DGELSD, (m, n, nrhs, tmp_data, m, pretval, maxmn,
2316 ps, rcon, rank, work.fortran_vec (), 2316 ps, rcon, rank, work.fortran_vec (),
2317 lwork, piwork, info)); 2317 lwork, piwork, info));
2467 2467
2468 Matrix atmp = *this; 2468 Matrix atmp = *this;
2469 double *tmp_data = atmp.fortran_vec (); 2469 double *tmp_data = atmp.fortran_vec ();
2470 2470
2471 double *pretval = retval.fortran_vec (); 2471 double *pretval = retval.fortran_vec ();
2472 Array<double> s (minmn, 1); 2472 Array<double> s (dim_vector (minmn, 1));
2473 double *ps = s.fortran_vec (); 2473 double *ps = s.fortran_vec ();
2474 2474
2475 // Ask DGELSD what the dimension of WORK should be. 2475 // Ask DGELSD what the dimension of WORK should be.
2476 octave_idx_type lwork = -1; 2476 octave_idx_type lwork = -1;
2477 2477
2478 Array<double> work (1, 1); 2478 Array<double> work (dim_vector (1, 1));
2479 2479
2480 octave_idx_type smlsiz; 2480 octave_idx_type smlsiz;
2481 F77_FUNC (xilaenv, XILAENV) (9, F77_CONST_CHAR_ARG2 ("DGELSD", 6), 2481 F77_FUNC (xilaenv, XILAENV) (9, F77_CONST_CHAR_ARG2 ("DGELSD", 6),
2482 F77_CONST_CHAR_ARG2 (" ", 1), 2482 F77_CONST_CHAR_ARG2 (" ", 1),
2483 0, 0, 0, 0, smlsiz 2483 0, 0, 0, 0, smlsiz
2498 nlvl = 0; 2498 nlvl = 0;
2499 2499
2500 octave_idx_type liwork = 3 * minmn * nlvl + 11 * minmn; 2500 octave_idx_type liwork = 3 * minmn * nlvl + 11 * minmn;
2501 if (liwork < 1) 2501 if (liwork < 1)
2502 liwork = 1; 2502 liwork = 1;
2503 Array<octave_idx_type> iwork (liwork, 1); 2503 Array<octave_idx_type> iwork (dim_vector (liwork, 1));
2504 octave_idx_type* piwork = iwork.fortran_vec (); 2504 octave_idx_type* piwork = iwork.fortran_vec ();
2505 2505
2506 F77_XFCN (dgelsd, DGELSD, (m, n, nrhs, tmp_data, m, pretval, maxmn, 2506 F77_XFCN (dgelsd, DGELSD, (m, n, nrhs, tmp_data, m, pretval, maxmn,
2507 ps, rcon, rank, work.fortran_vec (), 2507 ps, rcon, rank, work.fortran_vec (),
2508 lwork, piwork, info)); 2508 lwork, piwork, info));