Программа DECOMP

[Ketsyki]

Здравствйте. Пытаюсь оценить обусловленность матрицы при изменении входных параметров.
Изменяю значения только одного элемента матрицы и получаю обусловленности
при -2 - $2 \cdot 10^5$
при -2,9 - $3 \cdot 10^7$
и т.д. (т.е. обусловленность растет) Но потом получается такая штука:
при -2,9999 - $1 \cdot 10^{14}$
А при -2,999999 - $4 \cdot 10^{13}$
Может такое быть? Программа написана на паскале. Может, разрядной сетки не хватает?

[Ketsyki]

Вот сама программа Decomp

Используется синтаксис Pascal

Const    ndim  = 10;
          ndimS = 10;
          ndimC = 10;
          maxnfe   = 3000;

type
  float = extended;

type  floatvector   = array [1..ndim]           of float;
      floatmatrix   = array [1..ndim,1..ndim]   of float;
      rvecn         = array [1..6*ndim+3]       of float;
      ivec5         = array [1..5]              of integer;
      ivector       = array [1..ndim]           of integer;
      floatvectorS  = array [1..ndimS]          of float;
      floatvectorC  = array [1..ndimC]          of float;
      floatmatrixCS = array [1..ndimS,1..ndimC] of float;

procedure DECOMP(
                   n  : integer;
             var   A  : floatmatrix ;
             var cond : float   ;
             var ipvt : ivector;
             var work : floatvector
                                   );
  VAR

        nm1,i,j,k,kp1,kb,km1,m  : integer;
        ek,t,anorm,ynorm,znorm  : float;
begin
        ipvt[n]:=1;
        if n = 1  then  begin
                          cond:=1.0;
                          if a[1,1] <> 0.0 then exit
                                           else begin
                                                  cond:=1.0e+32;
                                                  exit
                                                end
                        end;
        nm1:=n-1;
        anorm:=0.0;
        for j:=1 to n do begin
                              t:=0.0;
                              for i:=1 to n do  t:=t+abs(a[i,j]);
                              if t > anorm then anorm:=t;
                         end;
{do 35} for k:=1 to  nm1 do begin
            kp1:=k + 1;
            m:=k;
            for i:=kp1 to n do
                         if abs(a[i,k]) > abs(a[m,k]) then m:=i;
            ipvt[k]:=m;
            if m <> k then ipvt[n]:=-ipvt[n];
            t:=a[m,k];
            a[m,k]:=a[k,k];
            a[k,k]:=t;
            if t <>  0.0 then begin
{20}        for i:=kp1 to n do a[i,k]:=-a[i,k]/t;
{do 30}     for j:=kp1 to n do   begin
                                   t:=a[m,j];
                                   a[m,j]:=a[k,j];
                                   a[k,j]:=t;
                                   if t <> 0.0 then for i:=kp1 to n do
{25}                                          a[i,j]:=a[i,j] + a[i,k]*t;
{30}                              end;
{35}                           end;
                            end;
{do 50}for k:=1 to n do begin
                                   t:=0.0;
                                   if k <> 1 then begin
                                                    km1:=k-1;
                                                    for i:=1 to km1 do
                                                 t:=t + a[i,k] * work[i];
                                                  end;
{45}                               ek:=1.0;
             if t < 0.0 then ek:=-1.0;
             if a[k,k] = 0.0 then begin
                                    cond:=1.0E+32;
                                    exit
                                  end;
             work[k]:=-(ek+t)  / a[k,k];
{50}                    end;
{do 60}for kb:=1 to nm1 do begin
                k:=n-kb;
                t:=0.0;
                kp1:=k + 1;
                for i:=kp1 to n do t:=t + a[i,k] * work[k];
                work[k]:=t;
                m:=ipvt[k];
                if m <> k then begin
                                 t:=work[m];
                                 work[m]:=work[k];
                                 work[k]:=t
                               end;
{60}                       end;
        ynorm:=0.0;
        for i:=1 to n do  ynorm:=ynorm + abs(work[i]);
        solve(n,a,work,ipvt);
        znorm:=0.0;
        for i:=1 to n do znorm:=znorm + abs(work[i]);
        cond:=anorm*znorm/ynorm;
        if cond < 1.0   then cond:=1.0;
        exit
END;