Ticket #13660: trac_13660-auto_convert_matrix_to_RDF_CDF.patch

File trac_13660-auto_convert_matrix_to_RDF_CDF.patch, 9.6 KB (added by Punarbasu Purkayastha, 10 years ago)

apply to devel/sage

  • sage/matrix/matrix2.pyx

    # HG changeset patch
    # User Punarbasu Purkayastha <ppurka@gmail.com>
    # Date 1351358170 -28800
    # Node ID ab35d2fe15ce282154db5bfa198c10563fea45e5
    # Parent  993f139bb53ea105f344e589debc3e94ae1f17f4
    auto convert matrices over RR and CC to RDF and CDF for eigenvalue/eigenvector computations
    
    diff --git a/sage/matrix/matrix2.pyx b/sage/matrix/matrix2.pyx
    a b  
    50955095            From: Number Field in a with defining polynomial y^4 - 2*y^3 - 507*y^2 + 4988*y - 8744
    50965096            To:   Algebraic Real Field
    50975097            Defn: a |--> 16.35066086057957?)
    5098            
    5099            
     5098
     5099
    51005100        Notice the effect of the extend option.
    5101        
     5101
    51025102        ::
    5103        
     5103
    51045104            sage: M=matrix(QQ,[[0,-1,0],[1,0,0],[0,0,2]])
    51055105            sage: M.eigenvalues()
    51065106            [2, -1*I, 1*I]
    51075107            sage: M.eigenvalues(extend=False)
    51085108            [2]
    51095109
     5110        Computation of the eigenvalues over the field ``RR`` or ``CC``
     5111        defaults to computing over ``RDF`` or ``CDF``, respectively,
     5112        instead. A warning is provided about this change.::
     5113
     5114            sage: M = matrix([[1.2, 2],[2, 3]])
     5115            sage: M.eigenvalues()
     5116            doctest:1: UserWarning: Computing the eigenvalues over the RDF field since the generic algorithm for the inexact ring RR may give incorrect results due to numerical precision issues.
     5117            [-0.0931712199461, 4.29317121995]
     5118            sage: M = matrix(CC, [[1.2, I],[2, 3]])
     5119            sage: M.eigenvalues()
     5120            doctest:1: UserWarning: Computing the eigenvalues over the CDF field since the generic algorithm for the inexact ring CC may give incorrect results due to numerical precision issues.
     5121            [0.881845698329 - 0.820914065343*I, 3.31815430167 + 0.820914065343*I]
     5122
    51105123        """
    51115124        x = self.fetch('eigenvalues')
     5125        self_ring = self.base_ring()
    51125126        if not x is None:
    51135127            if not extend:
    51145128                x1=Sequence([])
    51155129                for i in x:
    5116                     if i in self.base_ring():
     5130                    if i in self_ring:
    51175131                        x1.append(i)
    51185132                x=x1
    51195133            return x
    51205134
    5121         if not self.base_ring().is_exact():
     5135        if not self_ring.is_exact():
    51225136            from warnings import warn
    5123             warn("Using generic algorithm for an inexact ring, which will probably give incorrect results due to numerical precision issues.")
     5137            from sage.rings.real_mpfr import RealField_class
     5138            from sage.rings.complex_field import ComplexField_class
     5139            if isinstance(self_ring, RealField_class):
     5140                warn("Computing the eigenvalues over the RDF field since "
     5141                     "the generic algorithm for the inexact ring RR may give "
     5142                     "incorrect results due to numerical precision issues.")
     5143                # Over RDF we get both real and complex eigenvalues
     5144                V = self.change_ring(RDF).eigenvalues()
     5145                self.cache('eigenvalues', V)
     5146                return V
     5147            elif isinstance(self_ring, ComplexField_class):
     5148                warn("Computing the eigenvalues over the CDF field since "
     5149                     "the generic algorithm for the inexact ring CC may give "
     5150                     "incorrect results due to numerical precision issues.")
     5151                # Over CDF we get both real and complex eigenvalues
     5152                V = self.change_ring(CDF).eigenvalues()
     5153                self.cache('eigenvalues', V)
     5154                return V
     5155            else:
     5156                warn("Using generic algorithm for an inexact ring, which "
     5157                     "will probably give incorrect results due to numerical "
     5158                     "precision issues.")
    51245159
    51255160        from sage.rings.qqbar import QQbar
    5126         G = self.fcp()   # factored charpoly of self. 
     5161        G = self.fcp()   # factored charpoly of self.
    51275162        V = []
    51285163        i=0
    51295164        for h, e in G:
     
    51365171                    try:
    51375172                        alpha = F.gen(0).galois_conjugates(QQbar)
    51385173                    except AttributeError, TypeError:
    5139                         raise NotImplementedError, "eigenvalues() is not implemented for matrices with eigenvalues that are not in the fraction field of the base ring or in QQbar"
     5174                        raise NotImplementedError("eigenvalues() is not "
     5175                            "implemented for matrices with eigenvalues that "
     5176                            "are not in the fraction field of the base ring "
     5177                            "or in QQbar")
    51405178                    V.extend(alpha*e)
    51415179            i+=1
    51425180        V = Sequence(V)
     
    51455183        if not extend:
    51465184            V1=Sequence([])
    51475185            for i in V:
    5148                 if i in self.base_ring():
     5186                if i in self_ring:
    51495187                    V1.append(i)
    51505188            V=V1
    51515189        return V
     
    51835221            sage: delta = eval*evec - evec*A
    51845222            sage: abs(abs(delta)) < 1e-10
    51855223            True
    5186        
     5224
    51875225        Notice the difference between considering ring extensions or not.
    5188        
     5226
    51895227        ::
    5190            
     5228
    51915229            sage: M=matrix(QQ,[[0,-1,0],[1,0,0],[0,0,2]])
    51925230            sage: M.eigenvectors_left()
    51935231            [(2, [
     
    51985236            (0, 0, 1)
    51995237            ], 1)]
    52005238
     5239        Computation of the eigenvectors over the field ``RR`` or ``CC``
     5240        defaults to computing over ``RDF`` or ``CDF``, respectively,
     5241        instead. A warning is provided about this change.::
     5242
     5243            sage: M = matrix([[1.2, 2],[2, 3]])
     5244            sage: M.eigenvectors_left()
     5245            doctest:1: UserWarning: Computing the eigenvectors over the RDF field since the generic algorithm for the inexact ring RR may give incorrect results due to numerical precision issues.
     5246            [(-0.0931712199461, [(-0.839751355262, 0.542971142268)], 1),
     5247            (4.29317121995, [(-0.542971142268, -0.839751355262)], 1)]
     5248            sage: M = matrix(CC, [[1.2, I],[2, 3]])
     5249            sage: M.eigenvectors_left()
     5250            doctest:1: UserWarning: Computing the eigenvectors over the CDF field since the generic algorithm for the inexact ring CC may give incorrect results due to numerical precision issues.
     5251            [(0.881845698329 - 0.820914065343*I, [(0.915245825866, -0.145594698293 + 0.37566908585*I)], 1),
     5252            (3.31815430167 + 0.820914065343*I, [(0.616145932705 + 0.238794153033*I, 0.750560818381)], 1)]
     5253
    52015254        """
    52025255        x = self.fetch('eigenvectors_left')
    52035256        if not x is None:
    52045257            return x
    52055258
    5206         if not self.base_ring().is_exact():
     5259        self_ring = self.base_ring()
     5260        if not self_ring.is_exact():
    52075261            from warnings import warn
    5208             warn("Using generic algorithm for an inexact ring, which may result in garbage from numerical precision issues.")
     5262            from sage.rings.real_mpfr import RealField_class
     5263            from sage.rings.complex_field import ComplexField_class
     5264            if isinstance(self_ring, RealField_class):
     5265                warn("Computing the eigenvectors over the RDF field since "
     5266                     "the generic algorithm for the inexact ring RR may give "
     5267                     "incorrect results due to numerical precision issues.")
     5268                # Over RDF we get both real and complex eigenvalues
     5269                evec_eval_list = self.change_ring(RDF).eigenvectors_left()
     5270                self.cache('eigenvectors_left', evec_eval_list)
     5271                return evec_eval_list
     5272            elif isinstance(self_ring, ComplexField_class):
     5273                warn("Computing the eigenvectors over the CDF field since "
     5274                     "the generic algorithm for the inexact ring CC may give "
     5275                     "incorrect results due to numerical precision issues.")
     5276                # Over CDF we get both real and complex eigenvalues
     5277                evec_eval_list = self.change_ring(CDF).eigenvectors_left()
     5278                self.cache('eigenvectors_left', evec_eval_list)
     5279                return evec_eval_list
     5280            else:
     5281                warn("Using generic algorithm for an inexact ring, which "
     5282                     "may result in garbage from numerical "
     5283                     "precision issues.")
    52095284
    52105285        V = []
    52115286        from sage.rings.qqbar import QQbar
     
    52145289        evec_list=[]
    52155290        n = self._nrows
    52165291        evec_eval_list = []
    5217         F = self.base_ring().fraction_field()
     5292        F = self_ring.fraction_field()
    52185293        for ev in eigenspaces:
    52195294            eigval = ev[0]
    52205295            eigbasis = ev[1].basis()
    52215296            eigmult = ev[2]
    5222             if eigval in self.base_ring() or extend:
     5297            if eigval in self_ring or extend:
    52235298                if eigval.parent().fraction_field() == F:
    52245299                    evec_eval_list.append((eigval, eigbasis, eigmult))
    52255300                else:
    52265301                    try:
    52275302                        eigval_conj = eigval.galois_conjugates(QQbar)
    52285303                    except AttributeError, TypeError:
    5229                         raise NotImplementedError, "eigenvectors are not implemented for matrices with eigenvalues that are not in the fraction field of the base ring or in QQbar"
    5230                        
     5304                        raise NotImplementedError("eigenvectors are not "
     5305                            "implemented for matrices with eigenvalues that "
     5306                            "are not in the fraction field of the base ring "
     5307                            "or in QQbar")
     5308
    52315309                    for e in eigval_conj:
    52325310                        m = hom(eigval.parent(), e.parent(), e)
    52335311                        space = (e.parent())**n