Ticket #10848: trac_10848-hermitian-two-speed.patch

File trac_10848-hermitian-two-speed.patch, 5.0 KB (added by rbeezer, 9 years ago)
  • sage/matrix/matrix_double_dense.pyx

    # HG changeset patch
    # User Rob Beezer <beezer@ups.edu>
    # Date 1301091507 25200
    # Node ID da2e9029a33f25023b5f2ac62264d60461c51660
    # Parent  07a2e975404a08cddf524ef0ad6ac59217519672
    [mq]: hermitian-two-speed
    
    diff -r 07a2e975404a -r da2e9029a33f sage/matrix/matrix_double_dense.pyx
    a b  
    13341334        self.cache(key, b)
    13351335        return b
    13361336
    1337     def is_hermitian(self, tol = 1e-12):
     1337    def is_hermitian(self, tol = 1e-12, algorithm='orthonormal'):
    13381338        r"""
    13391339        Return ``True`` if the matrix is equal to its conjugate-transpose.
    13401340
     
    13441344          absolute value of the difference between two matrix entries
    13451345          for which they will still be considered equal.
    13461346
     1347        - ``algorithm`` - default: 'orthonormal' - set to 'orthonormal'
     1348          for a stable procedure and set to 'naive' for a fast
     1349          procedure.
     1350
    13471351        OUTPUT:
    13481352
    13491353        ``True`` if the matrix is square and equal to the transpose
     
    13641368            sage: A = matrix(CDF, [[ 1 + I,  1 - 6*I, -1 - I],
    13651369            ...                    [-3 - I,     -4*I,     -2],
    13661370            ...                    [-1 + I, -2 - 8*I,  2 + I]])
    1367             sage: A.is_hermitian()
     1371            sage: A.is_hermitian(algorithm='orthonormal')
    13681372            False
    13691373            sage: B = A*A.conjugate_transpose()
    1370             sage: B.is_hermitian()
     1374            sage: B.is_hermitian(algorithm='orthonormal')
    13711375            True
    13721376
    13731377        We get a unitary matrix from the SVD routine and use this
     
    13821386            ...                    [-1 + I, -2 - 8*I,  2 + I]])
    13831387            sage: U, _, _ = A.SVD()
    13841388            sage: B=U*U.conjugate_transpose()
    1385             sage: B.is_hermitian()
     1389            sage: B.is_hermitian(algorithm='orthonormal')
    13861390            True
    1387             sage: B.is_hermitian(tol=5.0e-17)
     1391            sage: B.is_hermitian(algorithm='orthonormal', tol=5.0e-17)
    13881392            True
    1389             sage: B.is_hermitian(tol=3.0e-17)
     1393            sage: B.is_hermitian(algorithm='orthonormal', tol=3.0e-17)
    13901394            False
    13911395
    1392         A matrix that is nearly Hermitian, but for a non-real
     1396        A matrix that is nearly Hermitian, but for one non-real
    13931397        diagonal entry. ::
    13941398
    13951399            sage: A = matrix(CDF, [[    2,   2-I, 1+4*I],
    13961400            ...                    [  2+I,   3+I, 2-6*I],
    13971401            ...                    [1-4*I, 2+6*I,     5]])
    1398             sage: A.is_hermitian()
     1402            sage: A.is_hermitian(algorithm='orthonormal')
    13991403            False
    14001404            sage: A[1,1] = 132
    1401             sage: A.is_hermitian()
     1405            sage: A.is_hermitian(algorithm='orthonormal')
    14021406            True
    14031407
    1404         Rectangular matrices are never Hermitian.  ::
     1408        Rectangular matrices are never Hermitian, no matter which
     1409        algorithm is requested.  ::
    14051410
    14061411            sage: A = matrix(CDF, 3, 4)
    14071412            sage: A.is_hermitian()
    14081413            False
     1414
     1415        TEST: algorithm keyword
    14091416        """
     1417        import sage.rings.complex_double
    14101418        global numpy
    14111419        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)
    14131424        b = self.fetch(key)
    14141425        if not b is None:
    14151426            return b
     
    14191430        if numpy is None:
    14201431            import numpy
    14211432        cdef Py_ssize_t i, j
     1433        cdef Matrix_double_dense T
    14221434        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:
    14271457                    break
    14281458        self.cache(key, hermitian)
    14291459        return hermitian