comparison liboctave/CMatrix.cc @ 7789:82be108cc558

First attempt at single precision tyeps * * * corrections to qrupdate single precision routines * * * prefer demotion to single over promotion to double * * * Add single precision support to log2 function * * * Trivial PROJECT file update * * * Cache optimized hermitian/transpose methods * * * Add tests for tranpose/hermitian and ChangeLog entry for new transpose code
author David Bateman <dbateman@free.fr>
date Sun, 27 Apr 2008 22:34:17 +0200
parents 45f5faba05a2
children 5861b95e9879
comparison
equal deleted inserted replaced
7788:45f5faba05a2 7789:82be108cc558
189 // double complex arguments. They have been modified by adding an 189 // double complex arguments. They have been modified by adding an
190 // implicit double precision (a-h,o-z) statement at the beginning of 190 // implicit double precision (a-h,o-z) statement at the beginning of
191 // each subroutine. 191 // each subroutine.
192 192
193 F77_RET_T 193 F77_RET_T
194 F77_FUNC (cffti, CFFTI) (const octave_idx_type&, Complex*); 194 F77_FUNC (zffti, ZFFTI) (const octave_idx_type&, Complex*);
195 195
196 F77_RET_T 196 F77_RET_T
197 F77_FUNC (cfftf, CFFTF) (const octave_idx_type&, Complex*, Complex*); 197 F77_FUNC (zfftf, ZFFTF) (const octave_idx_type&, Complex*, Complex*);
198 198
199 F77_RET_T 199 F77_RET_T
200 F77_FUNC (cfftb, CFFTB) (const octave_idx_type&, Complex*, Complex*); 200 F77_FUNC (zfftb, ZFFTB) (const octave_idx_type&, Complex*, Complex*);
201 201
202 F77_RET_T 202 F77_RET_T
203 F77_FUNC (zlartg, ZLARTG) (const Complex&, const Complex&, 203 F77_FUNC (zlartg, ZLARTG) (const Complex&, const Complex&,
204 double&, Complex&, Complex&); 204 double&, Complex&, Complex&);
205 205
880 octave_idx_type nr_insert = nr; 880 octave_idx_type nr_insert = nr;
881 ComplexMatrix retval (nr + a.rows (), nc); 881 ComplexMatrix retval (nr + a.rows (), nc);
882 retval.insert (*this, 0, 0); 882 retval.insert (*this, 0, 0);
883 retval.insert (a, nr_insert, 0); 883 retval.insert (a, nr_insert, 0);
884 return retval; 884 return retval;
885 }
886
887 ComplexMatrix
888 ComplexMatrix::hermitian (void) const
889 {
890 octave_idx_type nr = rows ();
891 octave_idx_type nc = cols ();
892 ComplexMatrix result;
893 if (length () > 0)
894 {
895 result.resize (nc, nr);
896 for (octave_idx_type j = 0; j < nc; j++)
897 for (octave_idx_type i = 0; i < nr; i++)
898 result.elem (j, i) = conj (elem (i, j));
899 }
900 return result;
901 } 885 }
902 886
903 ComplexMatrix 887 ComplexMatrix
904 conj (const ComplexMatrix& a) 888 conj (const ComplexMatrix& a)
905 { 889 {
1354 Complex *pwsave = wsave.fortran_vec (); 1338 Complex *pwsave = wsave.fortran_vec ();
1355 1339
1356 retval = *this; 1340 retval = *this;
1357 Complex *tmp_data = retval.fortran_vec (); 1341 Complex *tmp_data = retval.fortran_vec ();
1358 1342
1359 F77_FUNC (cffti, CFFTI) (npts, pwsave); 1343 F77_FUNC (zffti, ZFFTI) (npts, pwsave);
1360 1344
1361 for (octave_idx_type j = 0; j < nsamples; j++) 1345 for (octave_idx_type j = 0; j < nsamples; j++)
1362 { 1346 {
1363 OCTAVE_QUIT; 1347 OCTAVE_QUIT;
1364 1348
1365 F77_FUNC (cfftf, CFFTF) (npts, &tmp_data[npts*j], pwsave); 1349 F77_FUNC (zfftf, ZFFTF) (npts, &tmp_data[npts*j], pwsave);
1366 } 1350 }
1367 1351
1368 return retval; 1352 return retval;
1369 } 1353 }
1370 1354
1395 Complex *pwsave = wsave.fortran_vec (); 1379 Complex *pwsave = wsave.fortran_vec ();
1396 1380
1397 retval = *this; 1381 retval = *this;
1398 Complex *tmp_data = retval.fortran_vec (); 1382 Complex *tmp_data = retval.fortran_vec ();
1399 1383
1400 F77_FUNC (cffti, CFFTI) (npts, pwsave); 1384 F77_FUNC (zffti, ZFFTI) (npts, pwsave);
1401 1385
1402 for (octave_idx_type j = 0; j < nsamples; j++) 1386 for (octave_idx_type j = 0; j < nsamples; j++)
1403 { 1387 {
1404 OCTAVE_QUIT; 1388 OCTAVE_QUIT;
1405 1389
1406 F77_FUNC (cfftb, CFFTB) (npts, &tmp_data[npts*j], pwsave); 1390 F77_FUNC (zfftb, ZFFTB) (npts, &tmp_data[npts*j], pwsave);
1407 } 1391 }
1408 1392
1409 for (octave_idx_type j = 0; j < npts*nsamples; j++) 1393 for (octave_idx_type j = 0; j < npts*nsamples; j++)
1410 tmp_data[j] = tmp_data[j] / static_cast<double> (npts); 1394 tmp_data[j] = tmp_data[j] / static_cast<double> (npts);
1411 1395
1439 Complex *pwsave = wsave.fortran_vec (); 1423 Complex *pwsave = wsave.fortran_vec ();
1440 1424
1441 retval = *this; 1425 retval = *this;
1442 Complex *tmp_data = retval.fortran_vec (); 1426 Complex *tmp_data = retval.fortran_vec ();
1443 1427
1444 F77_FUNC (cffti, CFFTI) (npts, pwsave); 1428 F77_FUNC (zffti, ZFFTI) (npts, pwsave);
1445 1429
1446 for (octave_idx_type j = 0; j < nsamples; j++) 1430 for (octave_idx_type j = 0; j < nsamples; j++)
1447 { 1431 {
1448 OCTAVE_QUIT; 1432 OCTAVE_QUIT;
1449 1433
1450 F77_FUNC (cfftf, CFFTF) (npts, &tmp_data[npts*j], pwsave); 1434 F77_FUNC (zfftf, ZFFTF) (npts, &tmp_data[npts*j], pwsave);
1451 } 1435 }
1452 1436
1453 npts = nc; 1437 npts = nc;
1454 nsamples = nr; 1438 nsamples = nr;
1455 nn = 4*npts+15; 1439 nn = 4*npts+15;
1458 pwsave = wsave.fortran_vec (); 1442 pwsave = wsave.fortran_vec ();
1459 1443
1460 Array<Complex> tmp (npts); 1444 Array<Complex> tmp (npts);
1461 Complex *prow = tmp.fortran_vec (); 1445 Complex *prow = tmp.fortran_vec ();
1462 1446
1463 F77_FUNC (cffti, CFFTI) (npts, pwsave); 1447 F77_FUNC (zffti, ZFFTI) (npts, pwsave);
1464 1448
1465 for (octave_idx_type j = 0; j < nsamples; j++) 1449 for (octave_idx_type j = 0; j < nsamples; j++)
1466 { 1450 {
1467 OCTAVE_QUIT; 1451 OCTAVE_QUIT;
1468 1452
1469 for (octave_idx_type i = 0; i < npts; i++) 1453 for (octave_idx_type i = 0; i < npts; i++)
1470 prow[i] = tmp_data[i*nr + j]; 1454 prow[i] = tmp_data[i*nr + j];
1471 1455
1472 F77_FUNC (cfftf, CFFTF) (npts, prow, pwsave); 1456 F77_FUNC (zfftf, ZFFTF) (npts, prow, pwsave);
1473 1457
1474 for (octave_idx_type i = 0; i < npts; i++) 1458 for (octave_idx_type i = 0; i < npts; i++)
1475 tmp_data[i*nr + j] = prow[i]; 1459 tmp_data[i*nr + j] = prow[i];
1476 } 1460 }
1477 1461
1505 Complex *pwsave = wsave.fortran_vec (); 1489 Complex *pwsave = wsave.fortran_vec ();
1506 1490
1507 retval = *this; 1491 retval = *this;
1508 Complex *tmp_data = retval.fortran_vec (); 1492 Complex *tmp_data = retval.fortran_vec ();
1509 1493
1510 F77_FUNC (cffti, CFFTI) (npts, pwsave); 1494 F77_FUNC (zffti, ZFFTI) (npts, pwsave);
1511 1495
1512 for (octave_idx_type j = 0; j < nsamples; j++) 1496 for (octave_idx_type j = 0; j < nsamples; j++)
1513 { 1497 {
1514 OCTAVE_QUIT; 1498 OCTAVE_QUIT;
1515 1499
1516 F77_FUNC (cfftb, CFFTB) (npts, &tmp_data[npts*j], pwsave); 1500 F77_FUNC (zfftb, ZFFTB) (npts, &tmp_data[npts*j], pwsave);
1517 } 1501 }
1518 1502
1519 for (octave_idx_type j = 0; j < npts*nsamples; j++) 1503 for (octave_idx_type j = 0; j < npts*nsamples; j++)
1520 tmp_data[j] = tmp_data[j] / static_cast<double> (npts); 1504 tmp_data[j] = tmp_data[j] / static_cast<double> (npts);
1521 1505
1527 pwsave = wsave.fortran_vec (); 1511 pwsave = wsave.fortran_vec ();
1528 1512
1529 Array<Complex> tmp (npts); 1513 Array<Complex> tmp (npts);
1530 Complex *prow = tmp.fortran_vec (); 1514 Complex *prow = tmp.fortran_vec ();
1531 1515
1532 F77_FUNC (cffti, CFFTI) (npts, pwsave); 1516 F77_FUNC (zffti, ZFFTI) (npts, pwsave);
1533 1517
1534 for (octave_idx_type j = 0; j < nsamples; j++) 1518 for (octave_idx_type j = 0; j < nsamples; j++)
1535 { 1519 {
1536 OCTAVE_QUIT; 1520 OCTAVE_QUIT;
1537 1521
1538 for (octave_idx_type i = 0; i < npts; i++) 1522 for (octave_idx_type i = 0; i < npts; i++)
1539 prow[i] = tmp_data[i*nr + j]; 1523 prow[i] = tmp_data[i*nr + j];
1540 1524
1541 F77_FUNC (cfftb, CFFTB) (npts, prow, pwsave); 1525 F77_FUNC (zfftb, ZFFTB) (npts, prow, pwsave);
1542 1526
1543 for (octave_idx_type i = 0; i < npts; i++) 1527 for (octave_idx_type i = 0; i < npts; i++)
1544 tmp_data[i*nr + j] = prow[i] / static_cast<double> (npts); 1528 tmp_data[i*nr + j] = prow[i] / static_cast<double> (npts);
1545 } 1529 }
1546 1530