Mercurial > hg > octave-lyh

comparison liboctave/dMatrix.cc @ 10350:12884915a8e4

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
merge MArray classes & improve Array interface
author Jaroslav Hajek <highegg@gmail.com>
date Sat, 23 Jan 2010 21:41:03 +0100
parents 07ebe522dac2
children a3635bc1ea19
comparison
equal deleted inserted replaced
10349:d4d13389c957 10350:12884915a8e4
212 } 212 }
213 213
214 // Matrix class. 214 // Matrix class.
215 215
216 Matrix::Matrix (const RowVector& rv) 216 Matrix::Matrix (const RowVector& rv)
217 : MArray2<double> (Array2<double> (rv, 1, rv.length ())) 217 : MArray<double> (rv)
218 { 218 {
219 } 219 }
220 220
221 Matrix::Matrix (const ColumnVector& cv) 221 Matrix::Matrix (const ColumnVector& cv)
222 : MArray2<double> (Array2<double> (cv, cv.length (), 1)) 222 : MArray<double> (cv)
223 { 223 {
224 } 224 }
225 225
226 Matrix::Matrix (const DiagMatrix& a) 226 Matrix::Matrix (const DiagMatrix& a)
227 : MArray2<double> (a.rows (), a.cols (), 0.0) 227 : MArray<double> (a.rows (), a.cols (), 0.0)
228 { 228 {
229 for (octave_idx_type i = 0; i < a.length (); i++) 229 for (octave_idx_type i = 0; i < a.length (); i++)
230 elem (i, i) = a.elem (i, i); 230 elem (i, i) = a.elem (i, i);
231 } 231 }
232 232
233 Matrix::Matrix (const PermMatrix& a) 233 Matrix::Matrix (const PermMatrix& a)
234 : MArray2<double> (a.rows (), a.cols (), 0.0) 234 : MArray<double> (a.rows (), a.cols (), 0.0)
235 { 235 {
236 const Array<octave_idx_type> ia (a.pvec ()); 236 const Array<octave_idx_type> ia (a.pvec ());
237 octave_idx_type len = a.rows (); 237 octave_idx_type len = a.rows ();
238 if (a.is_col_perm ()) 238 if (a.is_col_perm ())
239 for (octave_idx_type i = 0; i < len; i++) 239 for (octave_idx_type i = 0; i < len; i++)
245 245
246 // FIXME -- could we use a templated mixed-type copy function 246 // FIXME -- could we use a templated mixed-type copy function
247 // here? 247 // here?
248 248
249 Matrix::Matrix (const boolMatrix& a) 249 Matrix::Matrix (const boolMatrix& a)
250 : MArray2<double> (a.rows (), a.cols ()) 250 : MArray<double> (a)
251 { 251 {
252 for (octave_idx_type i = 0; i < a.rows (); i++)
253 for (octave_idx_type j = 0; j < a.cols (); j++)
254 elem (i, j) = a.elem (i, j);
255 } 252 }
256 253
257 Matrix::Matrix (const charMatrix& a) 254 Matrix::Matrix (const charMatrix& a)
258 : MArray2<double> (a.rows (), a.cols ()) 255 : MArray<double> (a.rows (), a.cols ())
259 { 256 {
260 for (octave_idx_type i = 0; i < a.rows (); i++) 257 for (octave_idx_type i = 0; i < a.rows (); i++)
261 for (octave_idx_type j = 0; j < a.cols (); j++) 258 for (octave_idx_type j = 0; j < a.cols (); j++)
262 elem (i, j) = static_cast<unsigned char> (a.elem (i, j)); 259 elem (i, j) = static_cast<unsigned char> (a.elem (i, j));
263 } 260 }
294 } 291 }
295 292
296 Matrix& 293 Matrix&
297 Matrix::insert (const Matrix& a, octave_idx_type r, octave_idx_type c) 294 Matrix::insert (const Matrix& a, octave_idx_type r, octave_idx_type c)
298 { 295 {
299 Array2<double>::insert (a, r, c); 296 Array<double>::insert (a, r, c);
300 return *this; 297 return *this;
301 } 298 }
302 299
303 Matrix& 300 Matrix&
304 Matrix::insert (const RowVector& a, octave_idx_type r, octave_idx_type c) 301 Matrix::insert (const RowVector& a, octave_idx_type r, octave_idx_type c)
612 // extract row or column i. 609 // extract row or column i.
613 610
614 RowVector 611 RowVector
615 Matrix::row (octave_idx_type i) const 612 Matrix::row (octave_idx_type i) const
616 { 613 {
617 return MArray<double> (index (idx_vector (i), idx_vector::colon)); 614 return index (idx_vector (i), idx_vector::colon);
618 } 615 }
619 616
620 ColumnVector 617 ColumnVector
621 Matrix::column (octave_idx_type i) const 618 Matrix::column (octave_idx_type i) const
622 { 619 {
623 return MArray<double> (index (idx_vector::colon, idx_vector (i))); 620 return index (idx_vector::colon, idx_vector (i));
624 } 621 }
625 622
626 Matrix 623 Matrix
627 Matrix::inverse (void) const 624 Matrix::inverse (void) const
628 { 625 {
732 729
733 if (nr != nc || nr == 0 || nc == 0) 730 if (nr != nc || nr == 0 || nc == 0)
734 (*current_liboctave_error_handler) ("inverse requires square matrix"); 731 (*current_liboctave_error_handler) ("inverse requires square matrix");
735 else 732 else
736 { 733 {
737 Array<octave_idx_type> ipvt (nr); 734 Array<octave_idx_type> ipvt (nr, 1);
738 octave_idx_type *pipvt = ipvt.fortran_vec (); 735 octave_idx_type *pipvt = ipvt.fortran_vec ();
739 736
740 retval = *this; 737 retval = *this;
741 double *tmp_data = retval.fortran_vec (); 738 double *tmp_data = retval.fortran_vec ();
742 739
743 Array<double> z(1); 740 Array<double> z(1, 1);
744 octave_idx_type lwork = -1; 741 octave_idx_type lwork = -1;
745 742
746 // Query the optimum work array size. 743 // Query the optimum work array size.
747 F77_XFCN (dgetri, DGETRI, (nc, tmp_data, nr, pipvt, 744 F77_XFCN (dgetri, DGETRI, (nc, tmp_data, nr, pipvt,
748 z.fortran_vec (), lwork, info)); 745 z.fortran_vec (), lwork, info));
749 746
750 lwork = static_cast<octave_idx_type> (z(0)); 747 lwork = static_cast<octave_idx_type> (z(0));
751 lwork = (lwork < 2 *nc ? 2*nc : lwork); 748 lwork = (lwork < 2 *nc ? 2*nc : lwork);
752 z.resize (lwork); 749 z.resize (lwork, 1);
753 double *pz = z.fortran_vec (); 750 double *pz = z.fortran_vec ();
754 751
755 info = 0; 752 info = 0;
756 753
757 // Calculate the norm of the matrix, for later use. 754 // Calculate the norm of the matrix, for later use.
769 { 766 {
770 octave_idx_type dgecon_info = 0; 767 octave_idx_type dgecon_info = 0;
771 768
772 // Now calculate the condition number for non-singular matrix. 769 // Now calculate the condition number for non-singular matrix.
773 char job = '1'; 770 char job = '1';
774 Array<octave_idx_type> iz (nc); 771 Array<octave_idx_type> iz (nc, 1);
775 octave_idx_type *piz = iz.fortran_vec (); 772 octave_idx_type *piz = iz.fortran_vec ();
776 F77_XFCN (dgecon, DGECON, (F77_CONST_CHAR_ARG2 (&job, 1), 773 F77_XFCN (dgecon, DGECON, (F77_CONST_CHAR_ARG2 (&job, 1),
777 nc, tmp_data, nr, anorm, 774 nc, tmp_data, nr, anorm,
778 rcon, pz, piz, dgecon_info 775 rcon, pz, piz, dgecon_info
779 F77_CHAR_ARG_LEN (1))); 776 F77_CHAR_ARG_LEN (1)));
1004 nsamples = nc; 1001 nsamples = nc;
1005 } 1002 }
1006 1003
1007 octave_idx_type nn = 4*npts+15; 1004 octave_idx_type nn = 4*npts+15;
1008 1005
1009 Array<Complex> wsave (nn); 1006 Array<Complex> wsave (nn, 1);
1010 Complex *pwsave = wsave.fortran_vec (); 1007 Complex *pwsave = wsave.fortran_vec ();
1011 1008
1012 retval = ComplexMatrix (*this); 1009 retval = ComplexMatrix (*this);
1013 Complex *tmp_data = retval.fortran_vec (); 1010 Complex *tmp_data = retval.fortran_vec ();
1014 1011
1045 nsamples = nc; 1042 nsamples = nc;
1046 } 1043 }
1047 1044
1048 octave_idx_type nn = 4*npts+15; 1045 octave_idx_type nn = 4*npts+15;
1049 1046
1050 Array<Complex> wsave (nn); 1047 Array<Complex> wsave (nn, 1);
1051 Complex *pwsave = wsave.fortran_vec (); 1048 Complex *pwsave = wsave.fortran_vec ();
1052 1049
1053 retval = ComplexMatrix (*this); 1050 retval = ComplexMatrix (*this);
1054 Complex *tmp_data = retval.fortran_vec (); 1051 Complex *tmp_data = retval.fortran_vec ();
1055 1052
1089 nsamples = nc; 1086 nsamples = nc;
1090 } 1087 }
1091 1088
1092 octave_idx_type nn = 4*npts+15; 1089 octave_idx_type nn = 4*npts+15;
1093 1090
1094 Array<Complex> wsave (nn); 1091 Array<Complex> wsave (nn, 1);
1095 Complex *pwsave = wsave.fortran_vec (); 1092 Complex *pwsave = wsave.fortran_vec ();
1096 1093
1097 retval = ComplexMatrix (*this); 1094 retval = ComplexMatrix (*this);
1098 Complex *tmp_data = retval.fortran_vec (); 1095 Complex *tmp_data = retval.fortran_vec ();
1099 1096
1108 1105
1109 npts = nc; 1106 npts = nc;
1110 nsamples = nr; 1107 nsamples = nr;
1111 nn = 4*npts+15; 1108 nn = 4*npts+15;
1112 1109
1113 wsave.resize (nn); 1110 wsave.resize (nn, 1);
1114 pwsave = wsave.fortran_vec (); 1111 pwsave = wsave.fortran_vec ();
1115 1112
1116 Array<Complex> tmp (npts); 1113 Array<Complex> tmp (npts, 1);
1117 Complex *prow = tmp.fortran_vec (); 1114 Complex *prow = tmp.fortran_vec ();
1118 1115
1119 F77_FUNC (zffti, ZFFTI) (npts, pwsave); 1116 F77_FUNC (zffti, ZFFTI) (npts, pwsave);
1120 1117
1121 for (octave_idx_type j = 0; j < nsamples; j++) 1118 for (octave_idx_type j = 0; j < nsamples; j++)
1155 nsamples = nc; 1152 nsamples = nc;
1156 } 1153 }
1157 1154
1158 octave_idx_type nn = 4*npts+15; 1155 octave_idx_type nn = 4*npts+15;
1159 1156
1160 Array<Complex> wsave (nn); 1157 Array<Complex> wsave (nn, 1);
1161 Complex *pwsave = wsave.fortran_vec (); 1158 Complex *pwsave = wsave.fortran_vec ();
1162 1159
1163 retval = ComplexMatrix (*this); 1160 retval = ComplexMatrix (*this);
1164 Complex *tmp_data = retval.fortran_vec (); 1161 Complex *tmp_data = retval.fortran_vec ();
1165 1162
1177 1174
1178 npts = nc; 1175 npts = nc;
1179 nsamples = nr; 1176 nsamples = nr;
1180 nn = 4*npts+15; 1177 nn = 4*npts+15;
1181 1178
1182 wsave.resize (nn); 1179 wsave.resize (nn, 1);
1183 pwsave = wsave.fortran_vec (); 1180 pwsave = wsave.fortran_vec ();
1184 1181
1185 Array<Complex> tmp (npts); 1182 Array<Complex> tmp (npts, 1);
1186 Complex *prow = tmp.fortran_vec (); 1183 Complex *prow = tmp.fortran_vec ();
1187 1184
1188 F77_FUNC (zffti, ZFFTI) (npts, pwsave); 1185 F77_FUNC (zffti, ZFFTI) (npts, pwsave);
1189 1186
1190 for (octave_idx_type j = 0; j < nsamples; j++) 1187 for (octave_idx_type j = 0; j < nsamples; j++)
1277 mattype.mark_as_unsymmetric (); 1274 mattype.mark_as_unsymmetric ();
1278 typ = MatrixType::Full; 1275 typ = MatrixType::Full;
1279 } 1276 }
1280 else 1277 else
1281 { 1278 {
1282 Array<double> z (3 * nc); 1279 Array<double> z (3 * nc, 1);
1283 double *pz = z.fortran_vec (); 1280 double *pz = z.fortran_vec ();
1284 Array<octave_idx_type> iz (nc); 1281 Array<octave_idx_type> iz (nc, 1);
1285 octave_idx_type *piz = iz.fortran_vec (); 1282 octave_idx_type *piz = iz.fortran_vec ();
1286 1283
1287 F77_XFCN (dpocon, DPOCON, (F77_CONST_CHAR_ARG2 (&job, 1), 1284 F77_XFCN (dpocon, DPOCON, (F77_CONST_CHAR_ARG2 (&job, 1),
1288 nr, tmp_data, nr, anorm, 1285 nr, tmp_data, nr, anorm,
1289 rcon, pz, piz, info 1286 rcon, pz, piz, info
1301 else if (typ != MatrixType::Full) 1298 else if (typ != MatrixType::Full)
1302 (*current_liboctave_error_handler) ("det: invalid dense matrix type"); 1299 (*current_liboctave_error_handler) ("det: invalid dense matrix type");
1303 1300
1304 if (typ == MatrixType::Full) 1301 if (typ == MatrixType::Full)
1305 { 1302 {
1306 Array<octave_idx_type> ipvt (nr); 1303 Array<octave_idx_type> ipvt (nr, 1);
1307 octave_idx_type *pipvt = ipvt.fortran_vec (); 1304 octave_idx_type *pipvt = ipvt.fortran_vec ();
1308 1305
1309 Matrix atmp = *this; 1306 Matrix atmp = *this;
1310 double *tmp_data = atmp.fortran_vec (); 1307 double *tmp_data = atmp.fortran_vec ();
1311 1308
1328 { 1325 {
1329 if (calc_cond) 1326 if (calc_cond)
1330 { 1327 {
1331 // Now calc the condition number for non-singular matrix. 1328 // Now calc the condition number for non-singular matrix.
1332 char job = '1'; 1329 char job = '1';
1333 Array<double> z (4 * nc); 1330 Array<double> z (4 * nc, 1);
1334 double *pz = z.fortran_vec (); 1331 double *pz = z.fortran_vec ();
1335 Array<octave_idx_type> iz (nc); 1332 Array<octave_idx_type> iz (nc, 1);
1336 octave_idx_type *piz = iz.fortran_vec (); 1333 octave_idx_type *piz = iz.fortran_vec ();
1337 1334
1338 F77_XFCN (dgecon, DGECON, (F77_CONST_CHAR_ARG2 (&job, 1), 1335 F77_XFCN (dgecon, DGECON, (F77_CONST_CHAR_ARG2 (&job, 1),
1339 nc, tmp_data, nr, anorm, 1336 nc, tmp_data, nr, anorm,
1340 rcon, pz, piz, info 1337 rcon, pz, piz, info
1393 octave_idx_type info = 0; 1390 octave_idx_type info = 0;
1394 char norm = '1'; 1391 char norm = '1';
1395 char uplo = 'U'; 1392 char uplo = 'U';
1396 char dia = 'N'; 1393 char dia = 'N';
1397 1394
1398 Array<double> z (3 * nc); 1395 Array<double> z (3 * nc, 1);
1399 double *pz = z.fortran_vec (); 1396 double *pz = z.fortran_vec ();
1400 Array<octave_idx_type> iz (nc); 1397 Array<octave_idx_type> iz (nc, 1);
1401 octave_idx_type *piz = iz.fortran_vec (); 1398 octave_idx_type *piz = iz.fortran_vec ();
1402 1399
1403 F77_XFCN (dtrcon, DTRCON, (F77_CONST_CHAR_ARG2 (&norm, 1), 1400 F77_XFCN (dtrcon, DTRCON, (F77_CONST_CHAR_ARG2 (&norm, 1),
1404 F77_CONST_CHAR_ARG2 (&uplo, 1), 1401 F77_CONST_CHAR_ARG2 (&uplo, 1),
1405 F77_CONST_CHAR_ARG2 (&dia, 1), 1402 F77_CONST_CHAR_ARG2 (&dia, 1),
1421 octave_idx_type info = 0; 1418 octave_idx_type info = 0;
1422 char norm = '1'; 1419 char norm = '1';
1423 char uplo = 'L'; 1420 char uplo = 'L';
1424 char dia = 'N'; 1421 char dia = 'N';
1425 1422
1426 Array<double> z (3 * nc); 1423 Array<double> z (3 * nc, 1);
1427 double *pz = z.fortran_vec (); 1424 double *pz = z.fortran_vec ();
1428 Array<octave_idx_type> iz (nc); 1425 Array<octave_idx_type> iz (nc, 1);
1429 octave_idx_type *piz = iz.fortran_vec (); 1426 octave_idx_type *piz = iz.fortran_vec ();
1430 1427
1431 F77_XFCN (dtrcon, DTRCON, (F77_CONST_CHAR_ARG2 (&norm, 1), 1428 F77_XFCN (dtrcon, DTRCON, (F77_CONST_CHAR_ARG2 (&norm, 1),
1432 F77_CONST_CHAR_ARG2 (&uplo, 1), 1429 F77_CONST_CHAR_ARG2 (&uplo, 1),
1433 F77_CONST_CHAR_ARG2 (&dia, 1), 1430 F77_CONST_CHAR_ARG2 (&dia, 1),
1466 mattype.mark_as_unsymmetric (); 1463 mattype.mark_as_unsymmetric ();
1467 typ = MatrixType::Full; 1464 typ = MatrixType::Full;
1468 } 1465 }
1469 else 1466 else
1470 { 1467 {
1471 Array<double> z (3 * nc); 1468 Array<double> z (3 * nc, 1);
1472 double *pz = z.fortran_vec (); 1469 double *pz = z.fortran_vec ();
1473 Array<octave_idx_type> iz (nc); 1470 Array<octave_idx_type> iz (nc, 1);
1474 octave_idx_type *piz = iz.fortran_vec (); 1471 octave_idx_type *piz = iz.fortran_vec ();
1475 1472
1476 F77_XFCN (dpocon, DPOCON, (F77_CONST_CHAR_ARG2 (&job, 1), 1473 F77_XFCN (dpocon, DPOCON, (F77_CONST_CHAR_ARG2 (&job, 1),
1477 nr, tmp_data, nr, anorm, 1474 nr, tmp_data, nr, anorm,
1478 rcon, pz, piz, info 1475 rcon, pz, piz, info
1485 1482
1486 if (typ == MatrixType::Full) 1483 if (typ == MatrixType::Full)
1487 { 1484 {
1488 octave_idx_type info = 0; 1485 octave_idx_type info = 0;
1489 1486
1490 Array<octave_idx_type> ipvt (nr); 1487 Array<octave_idx_type> ipvt (nr, 1);
1491 octave_idx_type *pipvt = ipvt.fortran_vec (); 1488 octave_idx_type *pipvt = ipvt.fortran_vec ();
1492 1489
1493 if(anorm < 0.) 1490 if(anorm < 0.)
1494 anorm = atmp.abs().sum(). 1491 anorm = atmp.abs().sum().
1495 row(static_cast<octave_idx_type>(0)).max(); 1492 row(static_cast<octave_idx_type>(0)).max();
1496 1493
1497 Array<double> z (4 * nc); 1494 Array<double> z (4 * nc, 1);
1498 double *pz = z.fortran_vec (); 1495 double *pz = z.fortran_vec ();
1499 Array<octave_idx_type> iz (nc); 1496 Array<octave_idx_type> iz (nc, 1);
1500 octave_idx_type *piz = iz.fortran_vec (); 1497 octave_idx_type *piz = iz.fortran_vec ();
1501 1498
1502 F77_XFCN (dgetrf, DGETRF, (nr, nr, tmp_data, nr, pipvt, info)); 1499 F77_XFCN (dgetrf, DGETRF, (nr, nr, tmp_data, nr, pipvt, info));
1503 1500
1504 if (info != 0) 1501 if (info != 0)
1565 { 1562 {
1566 char norm = '1'; 1563 char norm = '1';
1567 char uplo = 'U'; 1564 char uplo = 'U';
1568 char dia = 'N'; 1565 char dia = 'N';
1569 1566
1570 Array<double> z (3 * nc); 1567 Array<double> z (3 * nc, 1);
1571 double *pz = z.fortran_vec (); 1568 double *pz = z.fortran_vec ();
1572 Array<octave_idx_type> iz (nc); 1569 Array<octave_idx_type> iz (nc, 1);
1573 octave_idx_type *piz = iz.fortran_vec (); 1570 octave_idx_type *piz = iz.fortran_vec ();
1574 1571
1575 F77_XFCN (dtrcon, DTRCON, (F77_CONST_CHAR_ARG2 (&norm, 1), 1572 F77_XFCN (dtrcon, DTRCON, (F77_CONST_CHAR_ARG2 (&norm, 1),
1576 F77_CONST_CHAR_ARG2 (&uplo, 1), 1573 F77_CONST_CHAR_ARG2 (&uplo, 1),
1577 F77_CONST_CHAR_ARG2 (&dia, 1), 1574 F77_CONST_CHAR_ARG2 (&dia, 1),
1665 { 1662 {
1666 char norm = '1'; 1663 char norm = '1';
1667 char uplo = 'L'; 1664 char uplo = 'L';
1668 char dia = 'N'; 1665 char dia = 'N';
1669 1666
1670 Array<double> z (3 * nc); 1667 Array<double> z (3 * nc, 1);
1671 double *pz = z.fortran_vec (); 1668 double *pz = z.fortran_vec ();
1672 Array<octave_idx_type> iz (nc); 1669 Array<octave_idx_type> iz (nc, 1);
1673 octave_idx_type *piz = iz.fortran_vec (); 1670 octave_idx_type *piz = iz.fortran_vec ();
1674 1671
1675 F77_XFCN (dtrcon, DTRCON, (F77_CONST_CHAR_ARG2 (&norm, 1), 1672 F77_XFCN (dtrcon, DTRCON, (F77_CONST_CHAR_ARG2 (&norm, 1),
1676 F77_CONST_CHAR_ARG2 (&uplo, 1), 1673 F77_CONST_CHAR_ARG2 (&uplo, 1),
1677 F77_CONST_CHAR_ARG2 (&dia, 1), 1674 F77_CONST_CHAR_ARG2 (&dia, 1),
1771 } 1768 }
1772 else 1769 else
1773 { 1770 {
1774 if (calc_cond) 1771 if (calc_cond)
1775 { 1772 {
1776 Array<double> z (3 * nc); 1773 Array<double> z (3 * nc, 1);
1777 double *pz = z.fortran_vec (); 1774 double *pz = z.fortran_vec ();
1778 Array<octave_idx_type> iz (nc); 1775 Array<octave_idx_type> iz (nc, 1);
1779 octave_idx_type *piz = iz.fortran_vec (); 1776 octave_idx_type *piz = iz.fortran_vec ();
1780 1777
1781 F77_XFCN (dpocon, DPOCON, (F77_CONST_CHAR_ARG2 (&job, 1), 1778 F77_XFCN (dpocon, DPOCON, (F77_CONST_CHAR_ARG2 (&job, 1),
1782 nr, tmp_data, nr, anorm, 1779 nr, tmp_data, nr, anorm,
1783 rcon, pz, piz, info 1780 rcon, pz, piz, info
1823 1820
1824 if (typ == MatrixType::Full) 1821 if (typ == MatrixType::Full)
1825 { 1822 {
1826 info = 0; 1823 info = 0;
1827 1824
1828 Array<octave_idx_type> ipvt (nr); 1825 Array<octave_idx_type> ipvt (nr, 1);
1829 octave_idx_type *pipvt = ipvt.fortran_vec (); 1826 octave_idx_type *pipvt = ipvt.fortran_vec ();
1830 1827
1831 Matrix atmp = *this; 1828 Matrix atmp = *this;
1832 double *tmp_data = atmp.fortran_vec (); 1829 double *tmp_data = atmp.fortran_vec ();
1833 if(anorm < 0.) 1830 if(anorm < 0.)
1834 anorm = atmp.abs().sum().row(static_cast<octave_idx_type>(0)).max(); 1831 anorm = atmp.abs().sum().row(static_cast<octave_idx_type>(0)).max();
1835 1832
1836 Array<double> z (4 * nc); 1833 Array<double> z (4 * nc, 1);
1837 double *pz = z.fortran_vec (); 1834 double *pz = z.fortran_vec ();
1838 Array<octave_idx_type> iz (nc); 1835 Array<octave_idx_type> iz (nc, 1);
1839 octave_idx_type *piz = iz.fortran_vec (); 1836 octave_idx_type *piz = iz.fortran_vec ();
1840 1837
1841 F77_XFCN (dgetrf, DGETRF, (nr, nr, tmp_data, nr, pipvt, info)); 1838 F77_XFCN (dgetrf, DGETRF, (nr, nr, tmp_data, nr, pipvt, info));
1842 1839
1843 // Throw-away extra info LAPACK gives so as to not change output. 1840 // Throw-away extra info LAPACK gives so as to not change output.
2262 2259
2263 Matrix atmp = *this; 2260 Matrix atmp = *this;
2264 double *tmp_data = atmp.fortran_vec (); 2261 double *tmp_data = atmp.fortran_vec ();
2265 2262
2266 double *pretval = retval.fortran_vec (); 2263 double *pretval = retval.fortran_vec ();
2267 Array<double> s (minmn); 2264 Array<double> s (minmn, 1);
2268 double *ps = s.fortran_vec (); 2265 double *ps = s.fortran_vec ();
2269 2266
2270 // Ask DGELSD what the dimension of WORK should be. 2267 // Ask DGELSD what the dimension of WORK should be.
2271 octave_idx_type lwork = -1; 2268 octave_idx_type lwork = -1;
2272 2269
2273 Array<double> work (1); 2270 Array<double> work (1, 1);
2274 2271
2275 octave_idx_type smlsiz; 2272 octave_idx_type smlsiz;
2276 F77_FUNC (xilaenv, XILAENV) (9, F77_CONST_CHAR_ARG2 ("DGELSD", 6), 2273 F77_FUNC (xilaenv, XILAENV) (9, F77_CONST_CHAR_ARG2 ("DGELSD", 6),
2277 F77_CONST_CHAR_ARG2 (" ", 1), 2274 F77_CONST_CHAR_ARG2 (" ", 1),
2278 0, 0, 0, 0, smlsiz 2275 0, 0, 0, 0, smlsiz
2300 nlvl = 0; 2297 nlvl = 0;
2301 2298
2302 octave_idx_type liwork = 3 * minmn * nlvl + 11 * minmn; 2299 octave_idx_type liwork = 3 * minmn * nlvl + 11 * minmn;
2303 if (liwork < 1) 2300 if (liwork < 1)
2304 liwork = 1; 2301 liwork = 1;
2305 Array<octave_idx_type> iwork (liwork); 2302 Array<octave_idx_type> iwork (liwork, 1);
2306 octave_idx_type* piwork = iwork.fortran_vec (); 2303 octave_idx_type* piwork = iwork.fortran_vec ();
2307 2304
2308 F77_XFCN (dgelsd, DGELSD, (m, n, nrhs, tmp_data, m, pretval, maxmn, 2305 F77_XFCN (dgelsd, DGELSD, (m, n, nrhs, tmp_data, m, pretval, maxmn,
2309 ps, rcon, rank, work.fortran_vec (), 2306 ps, rcon, rank, work.fortran_vec (),
2310 lwork, piwork, info)); 2307 lwork, piwork, info));
2345 if (work(0) < lworkaround) 2342 if (work(0) < lworkaround)
2346 work(0) = lworkaround; 2343 work(0) = lworkaround;
2347 } 2344 }
2348 2345
2349 lwork = static_cast<octave_idx_type> (work(0)); 2346 lwork = static_cast<octave_idx_type> (work(0));
2350 work.resize (lwork); 2347 work.resize (lwork, 1);
2351 2348
2352 F77_XFCN (dgelsd, DGELSD, (m, n, nrhs, tmp_data, m, pretval, 2349 F77_XFCN (dgelsd, DGELSD, (m, n, nrhs, tmp_data, m, pretval,
2353 maxmn, ps, rcon, rank, 2350 maxmn, ps, rcon, rank,
2354 work.fortran_vec (), lwork, 2351 work.fortran_vec (), lwork,
2355 piwork, info)); 2352 piwork, info));
2463 2460
2464 Matrix atmp = *this; 2461 Matrix atmp = *this;
2465 double *tmp_data = atmp.fortran_vec (); 2462 double *tmp_data = atmp.fortran_vec ();
2466 2463
2467 double *pretval = retval.fortran_vec (); 2464 double *pretval = retval.fortran_vec ();
2468 Array<double> s (minmn); 2465 Array<double> s (minmn, 1);
2469 double *ps = s.fortran_vec (); 2466 double *ps = s.fortran_vec ();
2470 2467
2471 // Ask DGELSD what the dimension of WORK should be. 2468 // Ask DGELSD what the dimension of WORK should be.
2472 octave_idx_type lwork = -1; 2469 octave_idx_type lwork = -1;
2473 2470
2474 Array<double> work (1); 2471 Array<double> work (1, 1);
2475 2472
2476 octave_idx_type smlsiz; 2473 octave_idx_type smlsiz;
2477 F77_FUNC (xilaenv, XILAENV) (9, F77_CONST_CHAR_ARG2 ("DGELSD", 6), 2474 F77_FUNC (xilaenv, XILAENV) (9, F77_CONST_CHAR_ARG2 ("DGELSD", 6),
2478 F77_CONST_CHAR_ARG2 (" ", 1), 2475 F77_CONST_CHAR_ARG2 (" ", 1),
2479 0, 0, 0, 0, smlsiz 2476 0, 0, 0, 0, smlsiz
2494 nlvl = 0; 2491 nlvl = 0;
2495 2492
2496 octave_idx_type liwork = 3 * minmn * nlvl + 11 * minmn; 2493 octave_idx_type liwork = 3 * minmn * nlvl + 11 * minmn;
2497 if (liwork < 1) 2494 if (liwork < 1)
2498 liwork = 1; 2495 liwork = 1;
2499 Array<octave_idx_type> iwork (liwork); 2496 Array<octave_idx_type> iwork (liwork, 1);
2500 octave_idx_type* piwork = iwork.fortran_vec (); 2497 octave_idx_type* piwork = iwork.fortran_vec ();
2501 2498
2502 F77_XFCN (dgelsd, DGELSD, (m, n, nrhs, tmp_data, m, pretval, maxmn, 2499 F77_XFCN (dgelsd, DGELSD, (m, n, nrhs, tmp_data, m, pretval, maxmn,
2503 ps, rcon, rank, work.fortran_vec (), 2500 ps, rcon, rank, work.fortran_vec (),
2504 lwork, piwork, info)); 2501 lwork, piwork, info));
2505 2502
2506 lwork = static_cast<octave_idx_type> (work(0)); 2503 lwork = static_cast<octave_idx_type> (work(0));
2507 work.resize (lwork); 2504 work.resize (lwork, 1);
2508 2505
2509 F77_XFCN (dgelsd, DGELSD, (m, n, nrhs, tmp_data, m, pretval, 2506 F77_XFCN (dgelsd, DGELSD, (m, n, nrhs, tmp_data, m, pretval,
2510 maxmn, ps, rcon, rank, 2507 maxmn, ps, rcon, rank,
2511 work.fortran_vec (), lwork, 2508 work.fortran_vec (), lwork,
2512 piwork, info)); 2509 piwork, info));
2829 } 2826 }
2830 2827
2831 Matrix 2828 Matrix
2832 Matrix::diag (octave_idx_type k) const 2829 Matrix::diag (octave_idx_type k) const
2833 { 2830 {
2834 return MArray2<double>::diag (k); 2831 return MArray<double>::diag (k);
2835 } 2832 }
2836 2833
2837 ColumnVector 2834 ColumnVector
2838 Matrix::row_min (void) const 2835 Matrix::row_min (void) const
2839 { 2836 {
2850 octave_idx_type nc = cols (); 2847 octave_idx_type nc = cols ();
2851 2848
2852 if (nr > 0 && nc > 0) 2849 if (nr > 0 && nc > 0)
2853 { 2850 {
2854 result.resize (nr); 2851 result.resize (nr);
2855 idx_arg.resize (nr); 2852 idx_arg.resize (nr, 1);
2856 2853
2857 for (octave_idx_type i = 0; i < nr; i++) 2854 for (octave_idx_type i = 0; i < nr; i++)
2858 { 2855 {
2859 octave_idx_type idx_j; 2856 octave_idx_type idx_j;
2860 2857
2905 octave_idx_type nc = cols (); 2902 octave_idx_type nc = cols ();
2906 2903
2907 if (nr > 0 && nc > 0) 2904 if (nr > 0 && nc > 0)
2908 { 2905 {
2909 result.resize (nr); 2906 result.resize (nr);
2910 idx_arg.resize (nr); 2907 idx_arg.resize (nr, 1);
2911 2908
2912 for (octave_idx_type i = 0; i < nr; i++) 2909 for (octave_idx_type i = 0; i < nr; i++)
2913 { 2910 {
2914 octave_idx_type idx_j; 2911 octave_idx_type idx_j;
2915 2912
2960 octave_idx_type nc = cols (); 2957 octave_idx_type nc = cols ();
2961 2958
2962 if (nr > 0 && nc > 0) 2959 if (nr > 0 && nc > 0)
2963 { 2960 {
2964 result.resize (nc); 2961 result.resize (nc);
2965 idx_arg.resize (nc); 2962 idx_arg.resize (1, nc);
2966 2963
2967 for (octave_idx_type j = 0; j < nc; j++) 2964 for (octave_idx_type j = 0; j < nc; j++)
2968 { 2965 {
2969 octave_idx_type idx_i; 2966 octave_idx_type idx_i;
2970 2967
3015 octave_idx_type nc = cols (); 3012 octave_idx_type nc = cols ();
3016 3013
3017 if (nr > 0 && nc > 0) 3014 if (nr > 0 && nc > 0)
3018 { 3015 {
3019 result.resize (nc); 3016 result.resize (nc);
3020 idx_arg.resize (nc); 3017 idx_arg.resize (1, nc);
3021 3018
3022 for (octave_idx_type j = 0; j < nc; j++) 3019 for (octave_idx_type j = 0; j < nc; j++)
3023 { 3020 {
3024 octave_idx_type idx_i; 3021 octave_idx_type idx_i;
3025 3022