Return to multprec_complex_linear_solvers.adb CVS log

Up to [local] / OpenXM_contrib / PHC / Ada / Math_Lib / Matrices

Import the second public release of PHCpack.

OKed by Jan Verschelde.

with Multprec_Complex_Numbers;           use Multprec_Complex_Numbers;

package body Multprec_Complex_Linear_Solvers is

-- AUXLILIARIES :

  function dconjg ( x : Complex_Number ) return Complex_Number is

  -- DESCRIPTION :
  --   Returns the complex conjugate of x.

    res : Complex_Number;
    re : Floating_Number := REAL_PART(x);
    im : Floating_Number := IMAG_PART(x);

  begin
    Min(im);
    res := Create(re,im);
    Clear(re); Clear(im);
    return res;
  end dconjg;

  function csign ( x,y : Complex_Number ) return Complex_Number is

  -- DESCRIPTION :
  --   Returns |x|*y/|y|.

    res : Complex_Number;
    fltacc : Floating_Number;
    cmpacc : Complex_Number;

  begin
    fltacc := AbsVal(x);
    res := Create(fltacc);
    Clear(fltacc);
    Mul(res,y);
    fltacc := AbsVal(y);
    cmpacc := Create(fltacc);
    Div(res,cmpacc);
    Clear(fltacc);
    Clear(cmpacc);
    return res;
  end csign;

  function Inverse_Abs_Sum ( z : Vector ) return Floating_Number is

  -- DESCRIPTION :
  --   Returns the reciprocal of the sum of the absolute values in z.

    res,sum,acc : Floating_Number;

  begin
    sum := Create(0);
    for i in z'range loop
      acc := AbsVal(z(i));
      Add(sum,acc);
      Clear(acc);
    end loop;
    res := Create(1);
    Div(res,sum);
    Clear(sum);
    return res;
  end Inverse_Abs_Sum;

-- TARGET ROUTINES :

  procedure Scale ( a : in out Matrix; b : in out Vector ) is

    fac : Complex_Number;

    function Maximum ( a : in Matrix; i : in integer ) return Complex_Number is

      res : integer := a'first(2);
      max : Floating_Number := AbsVal(a(i,res));
      tmp : Floating_Number;

    begin
      for j in a'first(2)+1..a'last(2) loop
        tmp := AbsVal(a(i,j));
        if tmp > max
         then max := tmp; res := j;
        end if;
      end loop;
      return a(i,res);
    end Maximum;

    procedure Divide ( a : in out Matrix; b : in out Vector;
                       i : in integer; fac : in Complex_Number ) is
    begin
      for j in a'range(2) loop
        a(i,j) := a(i,j)/fac;
      end loop;
      b(i) := b(i)/fac;
    end Divide;
  
  begin
    for i in a'range(1) loop
      fac := Maximum(a,i);
      Divide(a,b,i,fac);
    end loop;
  end Scale;

-- TARGET ROUTINES :

  procedure lufac ( a : in out Matrix; n : in integer;
                    ipvt : out Standard_Natural_Vectors.Vector;
                    info : out natural ) is

    kp1,l,nm1 : integer;
    smax,fltacc : Floating_Number;
    temp,cmpacc : Complex_Number;

  begin
    info := 0;
    nm1 := n - 1;
    if nm1 >= 1
     then for k in 1..nm1 loop
            kp1 := k + 1;
            l := k;
            smax := AbsVal(a(k,k));                  -- find the pivot index l
            for i in kp1..n loop
              fltacc := AbsVal(a(i,k));
              if fltacc > smax
               then l := i;
                    Copy(fltacc,smax);
              end if;
              Clear(fltacc);
            end loop;
            ipvt(k) := l;
            if Equal(smax,0.0)       -- this column is already triangularized
             then info := k;
             else if l /= k                       -- interchange if necessary
                   then Copy(a(l,k),temp);
                        Copy(a(k,k),a(l,k));
                        Copy(temp,a(k,k)); Clear(temp);
                  end if;
                  cmpacc := Create(1);                 -- compute multipliers
                  Div(cmpacc,a(k,k));
                  Min(cmpacc);
                  for i in kp1..n loop
                    Mul(a(i,k),cmpacc);
                  end loop;
                  Clear(cmpacc);
                  for j in kp1..n loop                     -- row elimination
                    Copy(a(l,j),temp);
                    if l /= k
                     then Copy(a(k,j),a(l,j));
                          Copy(temp,a(k,j));
                    end if;
                    for i in kp1..n loop
                      cmpacc := temp*a(i,k);
                      Add(a(i,j),cmpacc);
                      Clear(cmpacc);
                    end loop;
                  end loop;
                  Clear(temp);
            end if;
            Clear(smax);
         end loop;
    end if;
    ipvt(n) := n;
    fltacc := AbsVal(a(n,n));
    if Equal(fltacc,0.0)
     then info := n;
    end if;
    Clear(fltacc);
  end lufac;

  procedure lufco ( a : in out Matrix; n : in integer;
                    ipvt : out Standard_Natural_Vectors.Vector;
                    rcond : out Floating_Number ) is

  -- NOTE :
  --   rcond = 1/(norm(a)*(estimate of norm(inverse(a))))
  --   estimate = norm(z)/norm(y) where a*z = y and ctrans(a)*y = e.
  --   ctrans(a) is the conjugate transpose of a.
  --   The components of e are chosen to cause maximum local
  --   growth in teh elements of w where ctrans(u)*w = e.
  --   The vectors are frequently rescaled to avoid overflow.

    z : Multprec_Complex_Vectors.Vector(1..n);
    info,kb,kp1,l : integer;
    s,sm,sum,anorm,ynorm,fltacc1,fltacc2 : Floating_Number;
    ek,t,wk,wkm,cmpacc1,cmpacc2 : Complex_Number;
    ipvtt : Standard_Natural_Vectors.Vector(1..n);

  begin
    anorm := Create(0);                              -- compute 1-norm of a
    for j in 1..n loop
      sum := Create(0);
      for i in 1..n loop
        fltacc1 := AbsVal(a(i,j));
        Add(sum,fltacc1);
        Clear(fltacc1);
      end loop;
      if sum > anorm
       then Copy(sum,anorm);
      end if; 
      Clear(sum);
    end loop;
    lufac(a,n,ipvtt,info);                                        -- factor
    ipvt := ipvtt;
    ek := Create(1);                               -- solve ctrans(u)*w = e
    for j in 1..n loop
      z(j) := Create(0);
    end loop;
    for k in 1..n loop
      fltacc1 := AbsVal(z(k));
      if not Equal(fltacc1,0.0)
       then cmpacc1 := -z(k);
            Clear(ek);
            ek := csign(ek,cmpacc1);
            Clear(cmpacc1);
      end if;
      Clear(fltacc1);
      cmpacc1 := ek - z(k);
      fltacc1 := AbsVal(cmpacc1);
      Clear(cmpacc1);
      fltacc2 := AbsVal(a(k,k));
      if fltacc1 > fltacc2
       then s := fltacc2/fltacc1;
            cmpacc1 := Create(s);
            Mul(z,cmpacc1);
            Mul(ek,cmpacc1);
            Clear(cmpacc1);
            Clear(s);
      end if;
      Clear(fltacc1); Clear(fltacc2);
      Clear(wk);  wk := ek - z(k);
      Clear(wkm); wkm := ek + z(k);
      Min(wkm);
      s := AbsVal(wk);
      sm := AbsVal(wkm);
      fltacc1 := AbsVal(a(k,k));
      if Equal(fltacc1,0.0)
       then Clear(wk);  wk := Create(1);
            Clear(wkm); wkm := Create(1);
       else cmpacc1 := dconjg(a(k,k));
            Div(wk,cmpacc1);
            Div(wkm,cmpacc1);
            Clear(cmpacc1);
      end if;
      Clear(fltacc1);
      kp1 := k + 1;
      if kp1 <= n
       then for j in kp1..n loop
              cmpacc2 := dconjg(a(k,j));
              cmpacc1 := wkm*cmpacc2;
              Add(cmpacc1,z(j));
              fltacc1 := AbsVal(cmpacc1);
              Add(sm,fltacc1);
              Clear(fltacc1); 
              Clear(cmpacc1);
              cmpacc1 := wk*cmpacc2;
              Add(z(j),cmpacc1);
              Clear(cmpacc1);
              Clear(cmpacc2);
              fltacc1 := AbsVal(z(j));
              Add(s,fltacc1);
              Clear(fltacc1);
            end loop;
            if s < sm
             then t := wkm - wk;
                  Copy(wkm,wk);
                  for j in kp1..n loop
                    cmpacc1 := t*dconjg(a(k,j));
                    Add(z(j),cmpacc1);
                    Clear(cmpacc1);
                  end loop;
                  Clear(t);
            end if;
      end if;
      Copy(wk,z(k));
    end loop;
    s := Inverse_Abs_Sum(z);
    cmpacc1 := Create(s);
    Mul(z,cmpacc1);
    Clear(cmpacc1);
    Clear(s);
    for k in 1..n loop                           -- solve ctrans(l)*y = w
      kb := n+1-k;
      if kb < n
       then t := Create(0);
            for i in (kb+1)..n loop
              cmpacc1 := dconjg(a(i,kb))*z(i);
              Add(t,cmpacc1);
              Clear(cmpacc1);
            end loop;
            Add(z(kb),t);
      end if;
      fltacc1 := AbsVal(z(kb));
      if fltacc1 > 1.0 
       then s := Create(1);
            Div(s,fltacc1);
            cmpacc1 := Create(s);
            Mul(z,cmpacc1);
            Clear(cmpacc1);
            Clear(s);
      end if;
      Clear(fltacc1);
      l := ipvtt(kb);
      if l /= kb
       then Copy(z(l),t);
            Copy(z(kb),z(l));
            Copy(t,z(kb));
            Clear(t);
      end if;
    end loop;
    s := Inverse_Abs_Sum(z);
    cmpacc1 := Create(s);
    Clear(s);
    Mul(z,cmpacc1);
    Clear(cmpacc1);
    ynorm := Create(1);
    for k in 1..n loop                                    -- solve l*v = y
      l := ipvtt(k);
      if l /= k
       then Copy(z(l),t);
            Copy(z(k),z(l));
            Copy(t,z(k));
       else Copy(z(l),t);
      end if;
      if k < n
       then for i in (k+1)..n loop
              cmpacc1 := t*a(i,k);
              Add(z(i),cmpacc1);
              Clear(cmpacc1);
            end loop;
      end if;
      Clear(t);
      fltacc1 := AbsVal(z(k));
      if fltacc1 > 1.0
       then s := Create(1);
            Div(s,fltacc1);
            cmpacc1 := Create(s);
            Mul(z,cmpacc1);
            Clear(cmpacc1);
            Mul(ynorm,s);
            Clear(s);
      end if;
      Clear(fltacc1);
    end loop;
    s := Inverse_Abs_Sum(z);
    cmpacc1 := Create(s);
    Mul(z,cmpacc1);
    Clear(cmpacc1);
    Mul(ynorm,s);
    Clear(s);
    for k in 1..n loop                                    -- solve u*z = v
      kb := n+1-k;
      fltacc1 := AbsVal(z(kb));
      fltacc2 := AbsVal(a(kb,kb));
      if fltacc1 > fltacc2
       then s := fltacc2/fltacc1;
            cmpacc1 := Create(s);
            Mul(z,cmpacc1);
            Mul(ynorm,s);
            Clear(s);
      end if;
      if Equal(fltacc2,0.0)
       then Clear(z(kb));
            z(kb) := Create(1);
       else Div(z(kb),a(kb,kb));
      end if;
      Clear(fltacc1);
      Clear(fltacc2);
      t := -z(kb);
      for i in 1..(kb-1) loop
        cmpacc1 := t*a(i,kb);
        Add(z(i),cmpacc1);
        Clear(cmpacc1);
      end loop;
      Clear(t);
    end loop;
    s := Inverse_Abs_Sum(z);                             -- make znorm = 1.0
    cmpacc1 := Create(s);
    Mul(z,cmpacc1);
    Clear(cmpacc1);
    Mul(ynorm,s);
    if Equal(anorm,0.0)
     then rcond := Create(0);
     else rcond := ynorm/anorm;
    end if;
    Clear(ynorm); Clear(s); Clear(sm);
    Clear(ek); Clear(t); Clear(wk); Clear(wkm);
    Clear(z);
  end lufco;

  procedure lusolve ( a : in Matrix; n : in integer;
                      ipvt : in Standard_Natural_Vectors.Vector;
                      b : in out Vector ) is

    l,nm1,kb : integer;
    temp,acc : Complex_Number;
 
  begin
    nm1 := n-1;
    if nm1 >= 1                                             -- solve l*y = b
     then for k in 1..nm1 loop
            l := ipvt(k);
            Copy(b(l),temp);
            if l /= k 
             then Copy(b(k),b(l));
                  Copy(temp,b(k));
            end if;
            for i in (k+1)..n loop
              acc := temp*a(i,k);
              Add(b(i),acc);
              Clear(acc);
            end loop;
            Clear(temp);
          end loop;
    end if;
    for k in 1..n loop                                     -- solve u*x = y
      kb := n+1-k;
      Div(b(kb),a(kb,kb));
      temp := -b(kb);
      for j in 1..(kb-1) loop
        acc := temp*a(j,kb);
        Add(b(j),acc);
        Clear(acc);
      end loop;
      Clear(temp);
    end loop;
  end lusolve;

  procedure Triangulate ( a : in out Matrix; n,m : in integer ) is

    max,cbs : Floating_Number;
    temp : Complex_Number;
    pivot,k,kcolumn : integer;

  begin
    k := 1;
    kcolumn := 1;
    while (k <= n) and (kcolumn <= m) loop
      max := Create(0);                                        -- find pivot
      pivot := 0;
      for l in k..n loop
        cbs := AbsVal(a(l,kcolumn));
        if cbs > max
         then max := cbs;
              pivot := l;
        end if;
      end loop;
      if pivot = 0
       then kcolumn := kcolumn + 1;
       else if pivot /= k                       -- interchange if necessary
             then for i in 1..m loop
                    temp := a(pivot,i);
                    a(pivot,i) := a(k,i);
                    a(k,i) := temp;
                  end loop;
            end if;
            for j in (kcolumn+1)..m loop                   -- triangulate a
              a(k,j) := a(k,j) / a(k,kcolumn);
            end loop;
            a(k,kcolumn) := Create(1);
            for i in (k+1)..n loop
              for j in (kcolumn+1)..m loop
                a(i,j) := a(i,j) - a(i,kcolumn) * a(k,j);
              end loop;
              a(i,kcolumn) := Create(0);
            end loop;
            k := k + 1;
            kcolumn := kcolumn + 1;
      end if;
    end loop;
  end Triangulate;

  procedure Diagonalize ( a : in out Matrix; n,m : in integer ) is

    max : Floating_Number;
    temp : Complex_Number;
    pivot,k,kcolumn : integer;

  begin
    k := 1;
    kcolumn := 1;
    while (k <= n) and (kcolumn <= m) loop
      max := Create(0);                                          -- find pivot
      for l in k..n loop
        if AbsVal(a(l,kcolumn)) > max
         then max := AbsVal(a(l,kcolumn));
              pivot := l;
        end if;
      end loop;
      if Equal(max,0.0)
       then kcolumn := kcolumn + 1;
       else if pivot /= k                        -- interchange if necessary
             then for i in 1..m loop
                    temp := a(pivot,i);
                    a(pivot,i) := a(k,i);
                    a(k,i) := temp;
                  end loop;
            end if;
            for j in (kcolumn+1)..m loop                    -- diagonalize a
              a(k,j) := a(k,j) / a(k,kcolumn);
            end loop;
            a(k,kcolumn) := Create(1);
            for i in 1..(k-1) loop
              for j in (kcolumn+1)..m loop
                a(i,j) := a(i,j) - a(i,kcolumn) * a(k,j);
              end loop;
            end loop;
            for i in (k+1)..n loop
              for j in (kcolumn+1)..m loop
                a(i,j) := a(i,j) - a(i,kcolumn) * a(k,j);
              end loop;
            end loop;
            for j in 1..(k-1) loop
              a(j,kcolumn) := Create(0);
            end loop;
            for j in (k+1)..n loop
              a(j,kcolumn) := Create(0);
            end loop;
            k := k + 1;
            kcolumn := kcolumn + 1;
      end if;
    end loop;
  end Diagonalize;

end Multprec_Complex_Linear_Solvers;