Mercurial > hg > octave-lyh

comparison liboctave/dMatrix.cc @ 8336:9813c07ca946

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
make det take advantage of matrix type
author Jaroslav Hajek <highegg@gmail.com>
date Wed, 19 Nov 2008 15:26:39 +0100
parents 64cf956a109c
children e02242c54c49
comparison
equal deleted inserted replaced
8335:64cf956a109c 8336:9813c07ca946
49 #include "mx-base.h" 49 #include "mx-base.h"
50 #include "mx-m-dm.h" 50 #include "mx-m-dm.h"
51 #include "mx-dm-m.h" 51 #include "mx-dm-m.h"
52 #include "mx-inlines.cc" 52 #include "mx-inlines.cc"
53 #include "oct-cmplx.h" 53 #include "oct-cmplx.h"
54 #include "oct-norm.h"
54 #include "quit.h" 55 #include "quit.h"
55 56
56 #if defined (HAVE_FFTW3) 57 #if defined (HAVE_FFTW3)
57 #include "oct-fftw.h" 58 #include "oct-fftw.h"
58 #endif 59 #endif
1238 } 1239 }
1239 1240
1240 DET 1241 DET
1241 Matrix::determinant (octave_idx_type& info, double& rcon, int calc_cond) const 1242 Matrix::determinant (octave_idx_type& info, double& rcon, int calc_cond) const
1242 { 1243 {
1243 DET retval; 1244 MatrixType mattype (*this);
1245 return determinant (mattype, info, rcon, calc_cond);
1246 }
1247
1248 DET
1249 Matrix::determinant (MatrixType& mattype,
1250 octave_idx_type& info, double& rcon, int calc_cond) const
1251 {
1252 DET retval (1.0);
1244 1253
1245 octave_idx_type nr = rows (); 1254 octave_idx_type nr = rows ();
1246 octave_idx_type nc = cols (); 1255 octave_idx_type nc = cols ();
1247 1256
1248 if (nr == 0 || nc == 0) 1257 if (nr != nc)
1249 { 1258 (*current_liboctave_error_handler) ("matrix must be square");
1250 retval = DET (1.0, 0);
1251 }
1252 else 1259 else
1253 { 1260 {
1254 Array<octave_idx_type> ipvt (nr); 1261 int typ = mattype.type ();
1255 octave_idx_type *pipvt = ipvt.fortran_vec (); 1262
1256 1263 if (typ == MatrixType::Lower || typ == MatrixType::Upper)
1257 Matrix atmp = *this; 1264 {
1258 double *tmp_data = atmp.fortran_vec (); 1265 for (octave_idx_type i = 0; i < nc; i++)
1259 1266 retval *= elem (i,i);
1260 info = 0; 1267 }
1261 1268 else if (typ == MatrixType::Hermitian)
1262 // Calculate the norm of the matrix, for later use. 1269 {
1263 double anorm = 0; 1270 Matrix atmp = *this;
1264 if (calc_cond) 1271 double *tmp_data = atmp.fortran_vec ();
1265 anorm = atmp.abs().sum().row(static_cast<octave_idx_type>(0)).max(); 1272
1266 1273 info = 0;
1267 F77_XFCN (dgetrf, DGETRF, (nr, nr, tmp_data, nr, pipvt, info)); 1274 double anorm = 0;
1268 1275 if (calc_cond) anorm = xnorm (*this, 1);
1269 // Throw-away extra info LAPACK gives so as to not change output. 1276
1270 rcon = 0.0; 1277
1271 if (info != 0) 1278 char job = 'L';
1272 { 1279 F77_XFCN (dpotrf, DPOTRF, (F77_CONST_CHAR_ARG2 (&job, 1), nr,
1273 info = -1; 1280 tmp_data, nr, info
1274 retval = DET (); 1281 F77_CHAR_ARG_LEN (1)));
1275 } 1282
1276 else 1283 if (info != 0)
1277 { 1284 {
1278 if (calc_cond) 1285 rcon = 0.0;
1279 { 1286 mattype.mark_as_unsymmetric ();
1280 // Now calc the condition number for non-singular matrix. 1287 typ = MatrixType::Full;
1281 char job = '1'; 1288 }
1282 Array<double> z (4 * nc); 1289 else
1283 double *pz = z.fortran_vec (); 1290 {
1284 Array<octave_idx_type> iz (nc); 1291 Array<double> z (3 * nc);
1285 octave_idx_type *piz = iz.fortran_vec (); 1292 double *pz = z.fortran_vec ();
1286 1293 Array<octave_idx_type> iz (nc);
1287 F77_XFCN (dgecon, DGECON, (F77_CONST_CHAR_ARG2 (&job, 1), 1294 octave_idx_type *piz = iz.fortran_vec ();
1288 nc, tmp_data, nr, anorm, 1295
1289 rcon, pz, piz, info 1296 F77_XFCN (dpocon, DPOCON, (F77_CONST_CHAR_ARG2 (&job, 1),
1290 F77_CHAR_ARG_LEN (1))); 1297 nr, tmp_data, nr, anorm,
1291 } 1298 rcon, pz, piz, info
1292 1299 F77_CHAR_ARG_LEN (1)));
1293 if (info != 0) 1300
1294 { 1301 if (info != 0)
1295 info = -1; 1302 rcon = 0.0;
1296 retval = DET (); 1303
1297 } 1304 for (octave_idx_type i = 0; i < nc; i++)
1298 else 1305 retval *= elem (i,i);
1299 { 1306
1300 retval = DET (1.0); 1307 retval = retval.square ();
1301 1308 }
1302 for (octave_idx_type i = 0; i < nc; i++) 1309 }
1303 { 1310 else if (typ != MatrixType::Full)
1304 double c = atmp(i,i); 1311 (*current_liboctave_error_handler) ("det: invalid dense matrix type");
1305 retval *= (ipvt(i) != (i+1)) ? -c : c; 1312
1313 if (typ == MatrixType::Full)
1314 {
1315 Array<octave_idx_type> ipvt (nr);
1316 octave_idx_type *pipvt = ipvt.fortran_vec ();
1317
1318 Matrix atmp = *this;
1319 double *tmp_data = atmp.fortran_vec ();
1320
1321 info = 0;
1322
1323 // Calculate the norm of the matrix, for later use.
1324 double anorm = 0;
1325 if (calc_cond) anorm = xnorm (*this, 1);
1326
1327 F77_XFCN (dgetrf, DGETRF, (nr, nr, tmp_data, nr, pipvt, info));
1328
1329 // Throw-away extra info LAPACK gives so as to not change output.
1330 rcon = 0.0;
1331 if (info != 0)
1332 {
1333 info = -1;
1334 retval = DET ();
1335 }
1336 else
1337 {
1338 if (calc_cond)
1339 {
1340 // Now calc the condition number for non-singular matrix.
1341 char job = '1';
1342 Array<double> z (4 * nc);
1343 double *pz = z.fortran_vec ();
1344 Array<octave_idx_type> iz (nc);
1345 octave_idx_type *piz = iz.fortran_vec ();
1346
1347 F77_XFCN (dgecon, DGECON, (F77_CONST_CHAR_ARG2 (&job, 1),
1348 nc, tmp_data, nr, anorm,
1349 rcon, pz, piz, info
1350 F77_CHAR_ARG_LEN (1)));
1306 } 1351 }
1307 } 1352
1308 } 1353 if (info != 0)
1354 {
1355 info = -1;
1356 retval = DET ();
1357 }
1358 else
1359 {
1360 for (octave_idx_type i = 0; i < nc; i++)
1361 {
1362 double c = atmp(i,i);
1363 retval *= (ipvt(i) != (i+1)) ? -c : c;
1364 }
1365 }
1366 }
1367 }
1309 } 1368 }
1310 1369
1311 return retval; 1370 return retval;
1312 } 1371 }
1313 1372