Changeset 1092

Show
Ignore:
Timestamp:
02/27/08 22:34:27 (6 months ago)
Author:
otter
Message:

Math.Matrices.eigenValues/rank/singularValues lead to a crash, for matrices with zero dimensions.
This has been fixed

Files:
1 modified

Legend:

Unmodified
Added
Removed

  • branches/maintenance/3.0/Modelica/Math/package.mo

    r1083 r1092  
    88 
    99annotation ( 
    10      
    1110  Invisible=true, 
    1211  Icon(coordinateSystem(preserveAspectRatio=true, extent={{-100,-100},{100,100}}),  
    … …  
    413412  annotation ( 
    414413    version="0.8.1", 
    415     versionDate="2004-08-21" 
    416     , 
     414    versionDate="2004-08-21", 
    417415    Documentation(info="<HTML> 
    418416<h4>Library content</h4> 
    … …  
    791789 
    792790    annotation ( 
    793        
    794791      Documentation(info="<HTML> 
    795792<h4>Syntax</h4> 
    … …  
    829826<a href=\"Modelica://Modelica.Math.Matrices.LU_solve\">Matrices.LU_solve</a> 
    830827</HTML>")); 
     828 
    831829  protected 
    832830    Integer info; 
    … …  
    847845 
    848846    annotation ( 
    849        
    850847      Documentation(info="<HTML> 
    851848<h4>Syntax</h4> 
    … …  
    891888<a href=\"Modelica://Modelica.Math.Matrices.LU_solve2\">Matrices.LU_solve2</a> 
    892889</HTML>")); 
     890 
    893891  protected 
    894892    Integer info; 
    … …  
    909907 
    910908    annotation ( 
    911        
    912909      Documentation(info="<HTML> 
    913910<h4>Syntax</h4> 
    … …  
    935932</p> 
    936933</HTML>")); 
     934 
    937935  protected 
    938936    Integer info; 
    … …  
    956954 
    957955    annotation ( 
    958        
    959956      Documentation(info="<HTML> 
    960957<h4>Syntax</h4> 
    … …  
    990987</p> 
    991988</HTML>")); 
     989 
    992990  protected 
    993991    Integer info; 
    … …  
    13771375 
    13781376    annotation ( 
    1379        
    13801377      Documentation(info="<HTML> 
    13811378<h4>Syntax</h4> 
    … …  
    14311428</HTML> 
    14321429")); 
     1430 
    14331431  protected 
    14341432    Integer info; 
    … …  
    14361434    Boolean onlyEigenvalues = false; 
    14371435  algorithm 
     1436  if size(A,1) > 0 then 
    14381437    if onlyEigenvalues then 
    14391438       (eigenvalues[:, 1],eigenvalues[:, 2],info) := LAPACK.dgeev_eigenValues(A); 
    … …  
    14451444\"Matrices.eigenvalues\" is not possible, since the 
    14461445numerical algorithm does not converge."); 
     1446  end if; 
    14471447  end eigenValues; 
    14481448 
    … …  
    14571457 
    14581458    annotation ( 
    1459        
    14601459      Documentation(info="<HTML> 
    14611460<h4>Syntax</h4> 
    … …  
    14881487<a href=\"Modelica://Modelica.Math.Matrices.eigenValues\">Matrices.eigenValues</a> 
    14891488</HTML>")); 
     1489 
    14901490  protected 
    14911491    Integer n=size(eigenValues, 1); 
    … …  
    15131513    input Real A[:, :] "Matrix"; 
    15141514    output Real sigma[min(size(A, 1), size(A, 2))] "Singular values"; 
    1515     output Real U[size(A, 1), size(A, 1)]=zeros(size(A, 1), size(A, 1)) 
     1515    output Real U[size(A, 1), size(A, 1)]=identity(size(A, 1)) 
    15161516      "Left orthogonal matrix"; 
    1517     output Real VT[size(A, 2), size(A, 2)]=zeros(size(A, 2), size(A, 2)) 
     1517    output Real VT[size(A, 2), size(A, 2)]=identity(size(A, 2)) 
    15181518      "Transposed right orthogonal matrix "; 
    15191519 
    … …  
    15651565    Integer n=min(size(A, 1), size(A, 2)) "Number of singular values"; 
    15661566  algorithm 
    1567     (sigma,U,VT,info) := Matrices.LAPACK.dgesvd(A); 
     1567  if n>0 then 
     1568    (sigma,U,VT,info) := Modelica.Math.Matrices.LAPACK.dgesvd(A); 
    15681569    assert(info == 0, "The numerical algorithm to compute the 
    15691570singular value decomposition did not converge"); 
     1571  end if; 
    15701572  end singularValues; 
    15711573 
    … …  
    16311633      "If eps > 0, the singular values are checked against eps; otherwise eps=max(size(A))*norm(A)*Modelica.Constants.eps is used"; 
    16321634    output Integer result "Rank of matrix A"; 
     1635 
     1636    annotation (Documentation(info="<html> 
     1637   
     1638</html>")); 
    16331639  protected 
    16341640    Integer n=min(size(A, 1), size(A, 2)); 
    16351641    Integer i=n; 
    1636     Real sigma[n]=singularValues(A) "Singular values"; 
    1637     Real eps2=if eps > 0 then eps else max(size(A))*sigma[1]*Modelica.Constants.eps; 
     1642    Real sigma[n]; 
     1643    Real eps2; 
    16381644  algorithm 
    1639     result := n; 
    1640     while i > 0 loop 
    1641       if sigma[i] > eps2 then 
    1642         result := i; 
    1643         i := 0; 
    1644       end if; 
    1645       i := i - 1; 
    1646     end while; 
    1647  
    1648     annotation (Documentation(info="<html> 
    1649    
    1650 </html>")); 
     1645    result := 0; 
     1646    if n > 0 then 
     1647      sigma := Modelica.Math.Matrices.singularValues(A); 
     1648      eps2 := if eps > 0 then eps else max(size(A))*sigma[1]*Modelica.Constants.eps; 
     1649      while i > 0 loop 
     1650        if sigma[i] > eps2 then 
     1651          result := i; 
     1652          i := 0; 
     1653        end if; 
     1654        i := i - 1; 
     1655      end while; 
     1656    end if; 
    16511657  end rank; 
    16521658 
    … …  
    17161722       Implemented. 
    17171723</li> 
    1718 </html>") 
    1719       ); 
     1724</html>")); 
    17201725  algorithm 
    17211726 
    … …  
    17601765 
    17611766    annotation ( 
    1762        
    17631767      Documentation(info="<HTML> 
    17641768<p>This function computes</p> 
    … …  
    18231827</ul> 
    18241828</html>")); 
     1829 
    18251830  protected 
    18261831    parameter Integer nmax=21; 
    … …  
    19181923 
    19191924    annotation ( 
    1920        
    19211925      Documentation(info="<HTML> 
    19221926<p> 
    … …  
    20762080 
    20772081    annotation ( 
    2078        
    20792082      Documentation(info="<HTML> 
    20802083<p> 
    … …  
    21172120</ul> 
    21182121</html>")); 
     2122 
    21192123  algorithm 
    21202124    F := [A, B, zeros(na, nb); zeros(2*nb, na), zeros(2*nb, nb), [identity(nb); 
    … …  
    21492153 
    21502154      annotation ( 
    2151          
    21522155        Documentation(info="Lapack documentation 
    21532156    Purpose    
    … …  
    22262229                  have converged.    
    22272230")); 
     2231 
    22282232    external "Fortran 77" dgeev("N", "V", n, Awork, n, eigenReal, eigenImag, 
    22292233        eigenVectors, n, eigenVectors, n, work, size(work, 1), info)  
    … …  
    22502254 
    22512255      annotation ( 
    2252          
    22532256        Documentation(info="Lapack documentation 
    22542257    Purpose    
    … …  
    25292532 
    25302533      annotation ( 
    2531          
    25322534        Documentation(info="Lapack documentation 
    25332535  Purpose                                                                  
    … …  
    26222624          < 0:  if INFO = -i, the i-th argument had an illegal value       
    26232625                                                                          ")); 
     2626 
    26242627    end dgels_vec; 
    26252628 
    … …  
    26472650 
    26482651      annotation ( 
    2649          
    26502652        Documentation(info="Lapack documentation 
    26512653  Purpose                                                                
    … …  
    27382740          = 0:  successful exit                                          
    27392741          < 0:  if INFO = -i, the i-th argument had an illegal value    ")); 
     2742 
    27402743    end dgelsx_vec; 
    27412744 
    … …  
    27942797                  has been completed, but the factor U is exactly    
    27952798                  singular, so the solution could not be computed.    
    2796 ")     ); 
     2799")); 
    27972800 
    27982801    external "FORTRAN 77" dgesv(size(A, 1), size(B, 2), Awork, size(A, 1), ipiv, 
    … …  
    28152818Same as function LAPACK.dgesv, but right hand side is a vector and not a matrix. 
    28162819For details of the arguments, see documentation of dgesv. 
    2817 ")     ); 
     2820")); 
    28182821 
    28192822    external "FORTRAN 77" dgesv(size(A, 1), 1, Awork, size(A, 1), ipiv, x, size( 
    … …  
    29252928          = 0:  successful exit. 
    29262929          < 0:  if INFO = -i, the i-th argument had an illegal value. 
    2927 ")     ); 
     2930")); 
    29282931    end dgglse_vec; 
    29292932 
    … …  
    29892992   
    29902993                  completed unless i = N.    
    2991 ")     ); 
     2994")); 
    29922995 
    29932996    external "FORTRAN 77" dgtsv(size(diag, 1), size(B, 2), subdiagwork, 
    … …  
    30153018Same as function LAPACK.dgtsv, but right hand side is a vector and not a matrix. 
    30163019For details of the arguments, see documentation of dgtsv. 
    3017 ")     ); 
     3020")); 
    30183021 
    30193022    external "FORTRAN 77" dgtsv(size(diag, 1), 1, subdiagwork, diagwork, 
    … …  
    30993102Array elements marked * are not used by the routine; elements marked 
    31003103+ need not be set on entry, but are required by the routine to store 
    3101 elements of U because of fill-in resulting from the row interchanges.") 
    3102         ); 
     3104elements of U because of fill-in resulting from the row interchanges.")); 
    31033105 
    31043106    external "FORTRAN 77" dgbsv(n, kLower, kUpper, size(B, 2), Awork, size( 
    … …  
    31223124      annotation ( 
    31233125        Documentation(info="Lapack documentation:   
    3124 ")     ); 
     3126")); 
    31253127 
    31263128    external "FORTRAN 77" dgbsv(n, kLower, kUpper, 1, Awork, size(Awork, 1), 
    … …  
    32363238                  did not converge to zero. See the description of WORK    
    32373239                  above for details.    
    3238 ")     ); 
     3240")); 
    32393241 
    32403242    external "Fortran 77" dgesvd("A", "A", size(A, 1), size(A, 2), Awork, size( 
    … …  
    33513353                  did not converge to zero. See the description of WORK    
    33523354                  above for details.    
    3353 ")     ); 
     3355")); 
    33543356 
    33553357    external "Fortran 77" dgesvd("N", "N", size(A, 1), size(A, 2), Awork, size( 
    … …  
    34133415              singular, and division by zero will occur if it is used 
    34143416              to solve a system of equations. 
    3415 ")     ); 
     3417")); 
    34163418 
    34173419    external "FORTRAN 77" dgetrf(size(A, 1), size(A, 2), LU, size(A, 1), pivots, 
    … …  
    34783480        = 0:  successful exit 
    34793481        < 0:  if INFO = -i, the i-th argument had an illegal value 
    3480 ")       ); 
     3482")); 
    34813483 
    34823484    protected 
    … …  
    35463548        = 0:  successful exit 
    35473549        < 0:  if INFO = -i, the i-th argument had an illegal value 
    3548 ")     ); 
     3550")); 
    35493551 
    35503552    protected 
    … …  
    36073609        < 0:  if INFO = -i, the i-th argument had an illegal value 
    36083610        > 0:  if INFO = i, U(i,i) is exactly zero; the matrix is 
    3609               singular and its inverse could not be computed.") 
    3610         ); 
     3611              singular and its inverse could not be computed.")); 
    36113612    protected 
    36123613      Integer lwork=min(10, size(LU, 1))*size(LU, 1) "Length of work array"; 
    … …  
    36833684The matrix P is represented in jpvt as follows: If 
    36843685   jpvt(j) = i 
    3685 then the jth column of P is the ith canonical unit vector.") 
    3686         ); 
     3686then the jth column of P is the ith canonical unit vector.")); 
    36873687    protected 
    36883688      Integer info; 
    … …  
    37513751        = 0:  successful exit 
    37523752        < 0:  if INFO = -i, the i-th argument has an illegal value 
    3753 ")     ); 
     3753")); 
    37543754 
    37553755    protected 
    … …  
    38003800 
    38013801  annotation ( 
    3802      
    38033802    Icon(coordinateSystem( 
    38043803        preserveAspectRatio=true, 
    … …  
    38733872</p> 
    38743873</html>")); 
     3874 
    38753875external "C" y = sin(u); 
    38763876annotation(Library="ModelicaExternalC"); 
    … …  
    38843884 
    38853885  annotation ( 
    3886      
    38873886    Icon(coordinateSystem( 
    38883887        preserveAspectRatio=true, 
    … …  
    39533952</p> 
    39543953</html>")); 
     3954 
    39553955external "C" y = cos(u); 
    39563956annotation(Library="ModelicaExternalC"); 
    … …  
    39643964 
    39653965  annotation ( 
    3966      
    39673966    Icon(coordinateSystem( 
    39683967        preserveAspectRatio=true, 
    … …  
    40364035</p> 
    40374036</html>")); 
     4037 
    40384038external "C" y = tan(u); 
    40394039annotation(Library="ModelicaExternalC"); 
    … …  
    40474047 
    40484048  annotation ( 
    4049      
    40504049    Icon(coordinateSystem( 
    40514050        preserveAspectRatio=true, 
    … …  
    41204119</p> 
    41214120</html>")); 
     4121 
    41224122external "C" y = asin(u); 
    41234123annotation(Library="ModelicaExternalC"); 
    … …  
    41314131 
    41324132  annotation ( 
    4133      
    41344133    Icon(coordinateSystem( 
    41354134        preserveAspectRatio=true, 
    … …  
    42004199</p> 
    42014200</html>")); 
     4201 
    42024202external "C" y = acos(u); 
    42034203annotation(Library="ModelicaExternalC"); 
    … …  
    42114211 
    42124212  annotation ( 
    4213      
    42144213    Icon(coordinateSystem( 
    42154214        preserveAspectRatio=true, 
    … …  
    42744273</p> 
    42754274</html>")); 
     4275 
    42764276external "C" y = atan(u); 
    42774277annotation(Library="ModelicaExternalC"); 
    … …  
    42864286 
    42874287  annotation ( 
    4288      
    42894288    Icon(coordinateSystem( 
    42904289        preserveAspectRatio=true, 
    … …  
    43794378</HTML> 
    43804379")); 
     4380 
    43814381external "C" y = atan2(u1, u2); 
    43824382annotation(Library="ModelicaExternalC"); 
    … …  
    43944394 
    43954395  annotation ( 
    4396      
    43974396    Icon(coordinateSystem( 
    43984397        preserveAspectRatio=true, 
    … …  
    44834482</HTML> 
    44844483")); 
     4484 
    44854485protected 
    44864486  Real pi = Modelica.Constants.pi; 
    … …  
    44984498 
    44994499  annotation ( 
    4500      
    45014500    Icon(coordinateSystem( 
    45024501        preserveAspectRatio=true, 
    … …  
    45734572</p> 
    45744573</html>")); 
     4574 
    45754575external "C" y = sinh(u); 
    45764576annotation(Library="ModelicaExternalC"); 
    … …  
    45844584 
    45854585  annotation ( 
    4586      
    45874586    Icon(coordinateSystem( 
    45884587        preserveAspectRatio=true, 
    … …  
    46594658</p> 
    46604659</html>")); 
     4660 
    46614661external "C" y = cosh(u); 
    46624662annotation(Library="ModelicaExternalC"); 
    … …  
    46704670 
    46714671  annotation ( 
    4672      
    46734672    Icon(coordinateSystem( 
    46744673        preserveAspectRatio=true, 
    … …  
    47334732</p> 
    47344733</html>")); 
     4734 
    47354735external "C" y = tanh(u); 
    47364736annotation(Library="ModelicaExternalC"); 
    … …  
    47444744 
    47454745  annotation ( 
    4746      
    47474746    Icon(coordinateSystem( 
    47484747        preserveAspectRatio=true, 
    … …  
    48204819</p> 
    48214820</html>")); 
     4821 
    48224822algorithm 
    48234823  y :=Modelica.Math.log(u + sqrt(u*u + 1)); 
    … …  
    48324832 
    48334833  annotation ( 
    4834      
    48354834    Icon(coordinateSystem( 
    48364835        preserveAspectRatio=true, 
    … …  
    49304929 
    49314930  annotation ( 
    4932      
    49334931    Icon(coordinateSystem( 
    49344932        preserveAspectRatio=true, 
    … …  
    50035001</p> 
    50045002</html>")); 
     5003 
    50055004external "C" y = exp(u); 
    50065005annotation(Library="ModelicaExternalC"); 
    … …  
    50145013 
    50155014  annotation ( 
    5016      
    50175015    Icon(coordinateSystem( 
    50185016        preserveAspectRatio=true, 
    … …  
    51015099 
    51025100  annotation ( 
    5103      
    51045101    Icon(coordinateSystem( 
    51055102        preserveAspectRatio=true,