Mercurial > hg > octave-lyh

comparison liboctave/dMatrix.cc @ 3466:9ff5622c993e

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
[project @ 2000-01-21 19:30:46 by hodelas] Fixed expm to use direct application of scaling, permutation operations in step (2) and reverse step (2).
author hodelas
date Fri, 21 Jan 2000 19:30:46 +0000
parents 87721841efd7
children a2dc6de198f9
comparison
equal deleted inserted replaced
3465:996bb7ea4507 3466:9ff5622c993e
1354 // Preconditioning step 2: balancing; code follows development 1354 // Preconditioning step 2: balancing; code follows development
1355 // in AEPBAL 1355 // in AEPBAL
1356 1356
1357 double *p_m = m.fortran_vec (); 1357 double *p_m = m.fortran_vec ();
1358 1358
1359 Array<double> scale(nc); 1359 int info, ilo, ihi,ilos,ihis, ii, jj;
1360 double *pscale = scale.fortran_vec (); 1360 Array<double> dpermute(nc);
1361 1361 Array<double> dscale(nc);
1362 int info, ilo, ihi; 1362 double *dp;
1363 1363
1364 // both scale and permute 1364 char job = 'P'; // permutation first
1365 char job = 'B'; 1365 dp = dpermute.fortran_vec();
1366 1366 F77_XFCN (dgebal, DGEBAL, (&job, nc, p_m, nc, ilo, ihi, dp, info, 1L, 1L));
1367 F77_XFCN (dgebal, DGEBAL, (&job, nc, p_m, nc, ilo, ihi, pscale, info, 1367
1368 1L, 1L)); 1368 job = 'S'; // then scaling
1369 dp = dscale.fortran_vec();
1370 F77_XFCN (dgebal, DGEBAL, (&job, nc, p_m, nc, ilos, ihis, dp, info, 1L, 1L));
1369 1371
1370 if (f77_exception_encountered) 1372 if (f77_exception_encountered)
1371 { 1373 {
1372 (*current_liboctave_error_handler) ("unrecoverable error in dgebal"); 1374 (*current_liboctave_error_handler) ("unrecoverable error in dgebal");
1373 return retval;
1374 }
1375
1376 // construct balancing matrices D, Dinv
1377
1378 Matrix dmat = Matrix (nc, nc, 0.0);
1379 Matrix dinv = Matrix (nc, nc, 0.0);
1380
1381 for (int i = 0; i < nc; i++)
1382 dmat(i,i) = dinv(i,i) = 1.0;
1383
1384 // dgebak, dside=R => dmat := D*dmat
1385 char dside = 'R';
1386 F77_XFCN (dgebak, DGEBAK, (&job, &dside, nc, ilo, ihi, pscale, nc,
1387 dmat.fortran_vec(), nc, info, 1L, 1L));
1388
1389 if (f77_exception_encountered)
1390 {
1391 (*current_liboctave_error_handler) ("unrecoverable error in dgebak");
1392 return retval;
1393 }
1394
1395 // dgebak, dside=L => dinv := dinv*D^{-1}
1396 dside = 'L';
1397 F77_XFCN (dgebak, DGEBAK, (&job, &dside, nc, ilo, ihi, pscale, nc,
1398 dinv.fortran_vec(), nc, info, 1L, 1L));
1399
1400 if (f77_exception_encountered)
1401 {
1402 (*current_liboctave_error_handler) ("unrecoverable error in dgebak");
1403 return retval; 1375 return retval;
1404 } 1376 }
1405 1377
1406 // Preconditioning step 3: scaling. 1378 // Preconditioning step 3: scaling.
1407 1379
1470 retval = retval * retval; 1442 retval = retval * retval;
1471 sqpow--; 1443 sqpow--;
1472 } 1444 }
1473 1445
1474 // Reverse preconditioning step 2: inverse balancing. 1446 // Reverse preconditioning step 2: inverse balancing.
1475 retval = dmat*retval*dinv; 1447 // apply inverse scaling to computed exponential
1476 1448 for(ii = 0; ii < nc ; ii++)
1449 for(jj=0; jj < nc ; jj++)
1450 retval(ii,jj) *= dscale(ii)/dscale(jj);
1451
1452 // construct balancing permutation vector
1453 Array<int> ipermute(nc);
1454 for(ii=0 ; ii < nc ; ii++) ipermute(ii) = ii; // identity permutation
1455
1456 // leading permutations in forward order
1457 for( ii = 0 ; ii < (ilo-1) ; ii++)
1458 {
1459 int swapidx = ( (int) dpermute(ii) ) -1;
1460 int tmp = ipermute(ii);
1461 ipermute(ii) = ipermute( swapidx );
1462 ipermute(swapidx) = tmp;
1463 }
1464
1465 // trailing permutations must be done in reverse order
1466 for( ii = nc-1 ; ii >= ihi ; ii--)
1467 {
1468 int swapidx = ( (int) dpermute(ii) ) -1;
1469 int tmp = ipermute(ii);
1470 ipermute(ii) = ipermute( swapidx );
1471 ipermute(swapidx) = tmp;
1472 }
1473
1474 // construct inverse balancing permutation vector
1475 Array<int> invpvec(nc);
1476 for( ii = 0 ; ii < nc ; ii++)
1477 invpvec(ipermute(ii)) = ii ; // Thanks to R. A. Lippert for this method
1478
1479 Matrix tmpMat = retval;
1480 for( ii = 0 ; ii < nc ; ii ++)
1481 for( jj= 0 ; jj < nc ; jj++ )
1482 retval(ii,jj) = tmpMat(invpvec(ii),invpvec(jj));
1483
1477 // Reverse preconditioning step 1: fix trace normalization. 1484 // Reverse preconditioning step 1: fix trace normalization.
1478 1485
1479 if (trshift > 0.0) 1486 if (trshift > 0.0)
1480 retval = exp (trshift) * retval; 1487 retval = exp (trshift) * retval;
1481 1488