# HG changeset patch
# User Rob Beezer <beezer@ups.edu>
# Date 1301091507 25200
# Node ID da2e9029a33f25023b5f2ac62264d60461c51660
# Parent 07a2e975404a08cddf524ef0ad6ac59217519672
[mq]: hermitiantwospeed
diff r 07a2e975404a r da2e9029a33f sage/matrix/matrix_double_dense.pyx
a

b


1334  1334  self.cache(key, b) 
1335  1335  return b 
1336  1336  
1337   def is_hermitian(self, tol = 1e12): 
 1337  def is_hermitian(self, tol = 1e12, algorithm='orthonormal'): 
1338  1338  r""" 
1339  1339  Return ``True`` if the matrix is equal to its conjugatetranspose. 
1340  1340  
… 
… 

1344  1344  absolute value of the difference between two matrix entries 
1345  1345  for which they will still be considered equal. 
1346  1346  
 1347   ``algorithm``  default: 'orthonormal'  set to 'orthonormal' 
 1348  for a stable procedure and set to 'naive' for a fast 
 1349  procedure. 
 1350  
1347  1351  OUTPUT: 
1348  1352  
1349  1353  ``True`` if the matrix is square and equal to the transpose 
… 
… 

1364  1368  sage: A = matrix(CDF, [[ 1 + I, 1  6*I, 1  I], 
1365  1369  ... [3  I, 4*I, 2], 
1366  1370  ... [1 + I, 2  8*I, 2 + I]]) 
1367   sage: A.is_hermitian() 
 1371  sage: A.is_hermitian(algorithm='orthonormal') 
1368  1372  False 
1369  1373  sage: B = A*A.conjugate_transpose() 
1370   sage: B.is_hermitian() 
 1374  sage: B.is_hermitian(algorithm='orthonormal') 
1371  1375  True 
1372  1376  
1373  1377  We get a unitary matrix from the SVD routine and use this 
… 
… 

1382  1386  ... [1 + I, 2  8*I, 2 + I]]) 
1383  1387  sage: U, _, _ = A.SVD() 
1384  1388  sage: B=U*U.conjugate_transpose() 
1385   sage: B.is_hermitian() 
 1389  sage: B.is_hermitian(algorithm='orthonormal') 
1386  1390  True 
1387   sage: B.is_hermitian(tol=5.0e17) 
 1391  sage: B.is_hermitian(algorithm='orthonormal', tol=5.0e17) 
1388  1392  True 
1389   sage: B.is_hermitian(tol=3.0e17) 
 1393  sage: B.is_hermitian(algorithm='orthonormal', tol=3.0e17) 
1390  1394  False 
1391  1395  
1392   A matrix that is nearly Hermitian, but for a nonreal 
 1396  A matrix that is nearly Hermitian, but for one nonreal 
1393  1397  diagonal entry. :: 
1394  1398  
1395  1399  sage: A = matrix(CDF, [[ 2, 2I, 1+4*I], 
1396  1400  ... [ 2+I, 3+I, 26*I], 
1397  1401  ... [14*I, 2+6*I, 5]]) 
1398   sage: A.is_hermitian() 
 1402  sage: A.is_hermitian(algorithm='orthonormal') 
1399  1403  False 
1400  1404  sage: A[1,1] = 132 
1401   sage: A.is_hermitian() 
 1405  sage: A.is_hermitian(algorithm='orthonormal') 
1402  1406  True 
1403  1407  
1404   Rectangular matrices are never Hermitian. :: 
 1408  Rectangular matrices are never Hermitian, no matter which 
 1409  algorithm is requested. :: 
1405  1410  
1406  1411  sage: A = matrix(CDF, 3, 4) 
1407  1412  sage: A.is_hermitian() 
1408  1413  False 
 1414  
 1415  TEST: algorithm keyword 
1409  1416  """ 
 1417  import sage.rings.complex_double 
1410  1418  global numpy 
1411  1419  tol = float(tol) 
1412   key = 'hermitian_%s'%tol 
 1420  if not algorithm in ['naive', 'orthonormal']: 
 1421  raise ValueError("algorithm must be 'naive' or 'orthonormal', not {0}".format(algorithm)) 
 1422  
 1423  key = 'hermitian_{0}_{1}'.format(algorithm, tol) 
1413  1424  b = self.fetch(key) 
1414  1425  if not b is None: 
1415  1426  return b 
… 
… 

1419  1430  if numpy is None: 
1420  1431  import numpy 
1421  1432  cdef Py_ssize_t i, j 
 1433  cdef Matrix_double_dense T 
1422  1434  hermitian = True 
1423   for i from 0 < i < self._nrows: 
1424   for j from 0 <= j <= i: 
1425   if numpy.absolute(self.get_unsafe(i,j)  self.get_unsafe(j,i).conjugate()) > tol: 
1426   hermitian = False 
 1435  if algorithm == 'orthonormal': 
 1436  # Schur decomposition over CDF will be diagonal and real iff Hermitian 
 1437  _, T = self.schur(base_ring=sage.rings.complex_double.CDF) 
 1438  for i from 0 < i < T._nrows: 
 1439  for j from i+1 <= j <= T._ncols: 
 1440  if numpy.absolute(T.get_unsafe(i,j)) > tol: 
 1441  hermitian = False 
 1442  break 
 1443  if not hermitian: 
 1444  break 
 1445  if hermitian: 
 1446  for i from 0 < i < T._nrows: 
 1447  if numpy.absolute(numpy.imag(T.get_unsafe(i,i))) > tol: 
 1448  hermitian = False 
 1449  break 
 1450  elif algorithm == 'naive': 
 1451  for i from 0 < i < self._nrows: 
 1452  for j from 0 <= j <= i: 
 1453  if numpy.absolute(self.get_unsafe(i,j)  self.get_unsafe(j,i).conjugate()) > tol: 
 1454  hermitian = False 
 1455  break 
 1456  if not hermitian: 
1427  1457  break 
1428  1458  self.cache(key, hermitian) 
1429  1459  return hermitian 