fpc/packages/numlib/mdt.pas

{
    $Id$
    This file is part of the Numlib package.
    Copyright (c) 1986-2000 by
     Kees van Ginneken, Wil Kortsmit and Loek van Reij of the
     Computational centre of the Eindhoven University of Technology

    FPC port Code          by Marco van de Voort ([email protected])
             documentation by Michael van Canneyt ([email protected])

    Unit with undocumented procedures. If you recognize the algoritms, please
    contact one of the above people.

    This is probably only one actual procedure/algoritm, but (faster)
    variants exist for special matrices. (like band, tridiagonal etc)

    Contrary to the other undocumented units, this unit is exported in the
    DLL.

    See the file COPYING.FPC, included in this distribution,
    for details about the copyright.

    This program is distributed in the hope that it will be useful,
    but WITHOUT ANY WARRANTY; without even the implied warranty of
    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.

 **********************************************************************}
Unit mdt;

interface
{$I DIRECT.INC}

uses typ, dsl, omv;

Procedure mdtgen(n, rwidth: ArbInt; Var alu: ArbFloat; Var p: ArbInt;
                 Var ca:ArbFloat; Var term: ArbInt);

Procedure mdtgtr(n: ArbInt; Var l, d, u, l1, d1, u1, u2: ArbFloat;
                 Var p: boolean; Var ca: ArbFloat; Var term: ArbInt);

Procedure mdtgsy(n, rwidth: ArbInt; Var a: ArbFloat; Var pp:ArbInt;
                 Var qq:boolean; Var ca:ArbFloat; Var term:ArbInt);

Procedure mdtgpd(n, rwidth: ArbInt; Var al, ca: ArbFloat; Var term: ArbInt);

Procedure mdtgba(n, lb, rb, rwa: ArbInt; Var a: ArbFloat; rwl: ArbInt;
                 Var l:ArbFloat; Var p: ArbInt; Var ca: ArbFloat; Var term:ArbInt);

Procedure mdtgpb(n, lb, rwidth: ArbInt; Var al, ca: ArbFloat;
                 Var term: ArbInt);

Procedure mdtdtr(n: ArbInt; Var l, d, u, l1, d1, u1: ArbFloat;
                 Var term:ArbInt);

implementation

Procedure mdtgen(n, rwidth: ArbInt; Var alu: ArbFloat; Var p: ArbInt;
                 Var ca:ArbFloat; Var term: ArbInt);

Var 
         indi, indk, nsr, ind, i, j, k, indexpivot : ArbInt;
      normr, sumrowi, pivot, l, normt, maxim, h, s : ArbFloat;
                                   palu, sumrow, t : ^arfloat1;
                                                pp : ^arint1;
                                          singular : boolean;
Begin
  If (n<1) Or (rwidth<1) Then
    Begin
      term := 3;
     exit
    End; {wrong input}
  palu := @alu;
  pp := @p;
  nsr := n*sizeof(ArbFloat);
  getmem(sumrow, nsr);
  getmem(t, nsr);
  normr := 0;
  singular := false ;
  For i:=1 To n Do
    Begin
      ind := (i-1)*rwidth;
      pp^[i] := i;
      sumrowi := 0;
      For j:=1 To n Do
        sumrowi := sumrowi+abs(palu^[ind+j]);
      sumrow^[i] := sumrowi;
     h := 2*random-1;
     t^[i] := sumrowi*h;
      h := abs(h);
     If normr<h Then normr := h;
      If sumrowi=0 Then
        singular := true
    End; {i}
  For k:=1 To n Do
    Begin
      maxim := 0;
     indexpivot := k;
      For i:=k To n Do
        Begin
          ind := (i-1)*rwidth;
          sumrowi := sumrow^[i];
          If sumrowi <> 0 Then
            Begin
              h := abs(palu^[ind+k])/sumrowi;
              If maxim<h Then
                Begin
                  maxim := h;
                 indexpivot := i
                End {maxim<h}
            End {sumrow <> 0}
        End; {i}
      If maxim=0 Then
        singular := true
      Else
        Begin
          If indexpivot <> k Then
            Begin
              ind := (indexpivot-1)*rwidth;
              indk := (k-1)*rwidth;
              For j:=1 To n Do
                Begin
                  h := palu^[ind+j];
                  palu^[ind+j] := palu^[indk+j];
                  palu^[indk+j] := h
                End; {j}
              h := t^[indexpivot];
             t^[indexpivot] := t^[k];
              t^[k] := h;
             pp^[k] := indexpivot;
              sumrow^[indexpivot] := sumrow^[k]
            End; {indexpivot <> k}
          pivot := palu^[(k-1)*rwidth+k];
          For i:=k+1 To n Do
            Begin
              ind := (i-1)*rwidth;
              l := palu^[ind+k]/pivot;
              palu^[ind+k] := l;
              If l <> 0 Then
                Begin
                  For j:=k+1 To n Do
                    palu^[ind+j] := palu^[ind+j]-l*palu^[(k-1)*rwidth+j];
                  If Not singular Then
                    t^[i] := t^[i]-l*t^[k]
                End {l <> 0}
            End {i}
        End {maxim <> 0}
    End; {k}
    If Not singular Then
      Begin
        normt := 0;
        For i:=n Downto 1 Do
          Begin
            indi := (i-1)*rwidth;
            s := 0;
            For j:=i+1 To n Do
              s := s+t^[j]*palu^[indi+j];
            t^[i] := (t^[i]-s)/palu^[indi+i];
            h := abs(t^[i]);
            If normt<h Then
              normt := h
          End; {i}
        ca := normt/normr
      End; {not singular}
    If singular Then
      Begin
        term := 4;
       ca := giant
      End
    Else
      term := 1;
  freemem(sumrow, nsr);
  freemem(t, nsr)
End; {mdtgen}

Procedure mdtgtr(n: ArbInt; Var l, d, u, l1, d1, u1, u2: ArbFloat;
                 Var p: boolean; Var ca: ArbFloat; Var term: ArbInt);

Var 
                         i, j, k, nmin1, sr : ArbInt;
   normr, normt, sumrowi, h, lj, di, ui, ll : ArbFloat;
                                       sing : boolean;
           pd, pu, pd1, pu1, pu2, t, sumrow : ^arfloat1;
                                    pl, pl1 : ^arfloat2;
                                         pp : ^arbool1;
Begin
  If n<1 Then
    Begin
      term := 3;
     exit
    End; {wrong input}
  pl := @l;
 pd := @d;
 pu := @u;
  pl1 := @l1;
 pd1 := @d1;
 pu1 := @u1;
 pu2 := @u2;
 pp := @p;
  sr := sizeof(ArbFloat);
  move(pl^, pl1^, (n-1)*sr);
  move(pd^, pd1^, n*sr);
  move(pu^, pu1^, (n-1)*sr);
  getmem(t, n*sr);
  getmem(sumrow, n*sr);
  normr := 0;
 sing := false;
  nmin1 := n-1;
  For i:=1 To n Do
    Begin
      pp^[i] := false;
      If i=1 Then
        sumrowi := abs(pd^[1])+abs(pu^[1])
      Else
        If i=n Then
          sumrowi := abs(pl^[n])+abs(pd^[n])
        Else
          sumrowi := abs(pl^[i])+abs(pd^[i])+abs(pu^[i]);
      sumrow^[i] := sumrowi;
     h := 2*random-1;
     t^[i] := sumrowi*h;
      h := abs(h);
      If normr<h Then
        normr := h;
      If sumrowi=0 Then
        sing := true
    End; {i}
  j := 1;
  while (j <> n) Do
    Begin
      i := j;
     j := j+1;
     lj := pl1^[j];
      If lj <> 0 Then
        Begin
          di := pd1^[i];
          If di=0 Then
            pp^[i] := true
          Else
            pp^[i] := abs(di/sumrow^[i])<abs(lj/sumrow^[j]);
          If pp^[i] Then
            Begin
              ui := pu1^[i];
             pd1^[i] := lj;
              pu1^[i] := pd1^[j];
             pl1^[j] := di/lj;
             ll := pl1^[j];
              pd1^[j] := ui-ll*pd1^[j];
              If i<nmin1 Then
                Begin
                  pu2^[i] := pu1^[j];
                  pu1^[j] := -ll*pu2^[i]
                End; {i<nmin1}
              sumrow^[j] := sumrow^[i];
              If (Not sing) Then
                Begin
                  h := t^[i];
                 t^[i] := t^[j];
                  t^[j] := h-ll*t^[i]
                End {not sing}
            End {pp^[i]}
          Else
            Begin
              pl1^[j] := lj/di;
             ll := pl1^[j];
              pd1^[j] := pd1^[j]-ll*pu1^[i];
              If i<nmin1 Then
                pu2^[i] := 0;
              If (Not sing) Then
                t^[j] := t^[j]-ll*t^[i]
            End {not pp^[i]}
        End {lj<>0}
      Else
        Begin
          If i<nmin1 Then
            pu2^[i] := 0;
          If pd1^[i]=0 Then
            sing := true
        End {lj=0}
    End; {j}
  If pd1^[n]=0 Then
    sing := true;
  If (Not sing) Then
    Begin
      normt := 0;
      t^[n] := t^[n]/pd1^[n];
      h := abs(t^[n]);
      If normt<h Then
        normt := h;
      If n > 1 Then
        Begin
          t^[nmin1] := (t^[nmin1]-pu1^[nmin1]*t^[n])/pd1^[nmin1];
          h := abs(t^[nmin1])
        End; {n > 1}
      If normt<h Then
        normt := h;
      For i:=n-2 Downto 1 Do
        Begin
          t^[i] := (t^[i]-pu1^[i]*t^[i+1]-pu2^[i]*t^[i+2])/pd1^[i];
          h := abs(t^[i]);
          If normt<h Then
            normt := h
        End; {i}
      ca := normt/normr
    End; {not sing}
  If (sing) Then
    Begin
      term := 4;
     ca := giant
    End {sing}
  Else
    term := 1;
  freemem(t, n*sr);
  freemem(sumrow, n*sr)
End; {mdtgtr}

Procedure mdtgsy(n, rwidth: ArbInt; Var a: ArbFloat; Var pp:ArbInt;
                 Var qq:boolean; Var ca:ArbFloat; Var term:ArbInt);

Var 
   i, j, kmin1, k, kplus1, kmin2, imin2, nsr, nsi, nsb, ii,
   imin1, jmin1, indexpivot, iplus1, indi, indj, indk, indp       : ArbInt;
   ra, h, absh, maxim, pivot, ct, norma, sumrowi, normt, normr, s : ArbFloat;
                               alt, l, d, t, u, v, l1, d1, u1, t1 : ^arfloat1;
                                                                p : ^arint1;
                                                                q : ^arbool1;
Begin
  If (n<1) Or (rwidth<1) Then
    Begin
      term := 3;
     exit
    End; {if}
  alt := @a;
 p := @pp;
 q := @qq;
  nsr := n*sizeof(ArbFloat);
  nsi := n*sizeof(ArbInt);
  nsb := n*sizeof(boolean);
  getmem(l, nsr);
  getmem(d, nsr);
  getmem(t, nsr);
  getmem(u, nsr);
  getmem(v, nsr);
  getmem(l1, nsr);
  getmem(d1, nsr);
  getmem(u1, nsr);
  getmem(t1, nsr);
  norma := 0;
  For i:=1 To n Do
    Begin
      indi := (i-1)*rwidth;
      p^[i] := i;
     sumrowi := 0;
      For j:=1 To i Do
        sumrowi := sumrowi+abs(alt^[indi+j]);
      For j:=i+1 To n Do
        sumrowi := sumrowi+abs(alt^[(j-1)*rwidth+i]);
      If norma<sumrowi Then
        norma := sumrowi
    End; {i}
  kmin1 := -1;
 k := 0;
 kplus1 := 1;
  while k<n Do
    Begin
      kmin2 := kmin1;
     kmin1 := k;
     k := kplus1;
     kplus1 := kplus1+1;
      indk := kmin1*rwidth;
      If k>3 Then
        Begin
          t^[2] := alt^[rwidth+2]*alt^[indk+1]+alt^[2*rwidth+2]*alt^[indk+2];
          For i:=3 To kmin2 Do
            Begin
              indi := (i-1)*rwidth;
              t^[i] := alt^[indi+i-1]*alt^[indk+i-2]+alt^[indi+i]
                       *alt^[indk+i-1]+alt^[indi+rwidth+i]*alt^[indk+i]
            End; {i}
          t^[kmin1] := alt^[indk-rwidth+kmin2]*alt^[indk+k-3]
                       +alt^[indk-rwidth+kmin1]*alt^[indk+kmin2]
                       +alt^[indk+kmin1];
          h := alt^[indk+k];
          For j:=2 To kmin1 Do
            h := h-t^[j]*alt^[indk+j-1];
          t^[k] := h;
          alt^[indk+k] := h-alt^[indk+kmin1]*alt^[indk+kmin2]
        End {k>3}
      Else
       If k=3 Then
        Begin
          t^[2] := alt^[rwidth+2]*alt^[2*rwidth+1]+alt^[2*rwidth+2];
          h := alt^[2*rwidth+3]-t^[2]*alt^[2*rwidth+1];
          t^[3] := h;
          alt^[2*rwidth+3] := h-alt^[2*rwidth+2]*alt^[2*rwidth+1]
        End  {k=3}
      Else
       If k=2 Then
        t^[2] := alt^[rwidth+2];
      maxim := 0;
      For i:=kplus1 To n Do
        Begin
          indi := (i-1)*rwidth;
          h := alt^[indi+k];
          For j:=2 To k Do
            h := h-t^[j]*alt^[indi+j-1];
          absh := abs(h);
          If maxim<absh Then
            Begin
              maxim := absh;
             indexpivot := i
            End; {if}
          alt^[indi+k] := h
        End; {i}
      If maxim <> 0 Then
        Begin
          If indexpivot>kplus1 Then
            Begin
              indp := (indexpivot-1)*rwidth;
              indk := k*rwidth;
              p^[kplus1] := indexpivot;
              For j:=1 To k Do
                Begin
                  h := alt^[indk+j];
                  alt^[indk+j] := alt^[indp+j];
                  alt^[indp+j] := h
                End; {j}
              For j:=indexpivot Downto kplus1 Do
                Begin
                  indj := (j-1)*rwidth;
                  h := alt^[indj+kplus1];
                  alt^[indj+kplus1] := alt^[indp+j];
                  alt^[indp+j] := h
                End; {j}
              For i:=indexpivot To n Do
                Begin
                  indi := (i-1)*rwidth;
                  h := alt^[indi+kplus1];
                  alt^[indi+kplus1] := alt^[indi+indexpivot];
                  alt^[indi+indexpivot] := h
                End  {i}
            End; {if}
          pivot := alt^[k*rwidth+k];
          For i:=k+2 To n Do
            alt^[(i-1)*rwidth+k] := alt^[(i-1)*rwidth+k]/pivot
        End {maxim <> 0}
    End; {k}
  d^[1] := alt^[1];
 i := 1;
  while i<n Do
    Begin
      imin1 := i;
     i := i+1;
      u^[imin1] := alt^[(i-1)*rwidth+imin1];
      l^[imin1] := u^[imin1];
     d^[i] := alt^[(i-1)*rwidth+i]
    End; {i}
  mdtgtr(n, l^[1], d^[1], u^[1], l1^[1], d1^[1], u1^[1], v^[1],
         q^[1], ct, term);
  alt^[1] := d1^[1];
 alt^[rwidth+1] := l1^[1];
  alt^[rwidth+2] := d1^[2];
 alt^[2] := u1^[1];
  imin1 := 1;
 i := 2;
  while i<n Do
    Begin
      imin2 := imin1;
     imin1 := i;
     i := i+1;
      indi := imin1*rwidth;
      alt^[indi+imin1] := l1^[imin1];
     alt^[indi+i] := d1^[i];
      alt^[(imin1-1)*rwidth+i] := u1^[imin1];
      alt^[(imin2-1)*rwidth+i] := v^[imin2]
    End; {i}
  If term=1 Then
    Begin
      normr := 0;
      For i:=1 To n Do
        Begin
          t^[i] := 2*random-1;
         h := t^[i];
          h := abs(h);
          If normr<h Then
            normr := h
        End; {i}
      i := 0;
      while i<n Do
        Begin
          imin1 := i;
         i := i+1;
         j := 1;
         h := t^[i];
          while j<imin1 Do
            Begin
              jmin1 := j;
             j := j+1;
              h := h-alt^[(i-1)*rwidth+jmin1]*t^[j]
            End; {j}
          t^[i] := h
        End; {i}
      dslgtr(n, l1^[1], d1^[1], u1^[1], v^[1], q^[1], t^[1], t1^[1], term);
      i := n+1;
     imin1 := n;
     normt := 0;
      while i>2 Do
        Begin
          iplus1 := i;
         i := imin1;
         imin1 := imin1-1;
         h := t1^[i];
          For j:=iplus1 To n Do
            h := h-alt^[(j-1)*rwidth+imin1]*t1^[j];
          t1^[i] := h;
         h := abs(h);
          If normt<h Then
            normt := h
        End; {i}
      ca := norma*normt/normr
    End {term=1}
  Else ca := giant;
  freemem(l, nsr);
  freemem(d, nsr);
  freemem(t, nsr);
  freemem(u, nsr);
  freemem(v, nsr);
  freemem(l1, nsr);
  freemem(d1, nsr);
  freemem(u1, nsr);
  freemem(t1, nsr)
End; {mdtgsy}

Procedure mdtgpd(n, rwidth: ArbInt; Var al, ca: ArbFloat; Var term: ArbInt);

Var 
    posdef                               : boolean;
    i, j, k, kmin1, indk, indi           : ArbInt;
    h, lkk, normr, normt, sumrowi, norma : ArbFloat;
    pal, t                               : ^arfloat1;
Begin
  If (n<1) Or (rwidth<1) Then
    Begin
      term := 3;
     exit
    End; {wrong input}
  getmem(t, sizeof(ArbFloat)*n);
  pal := @al;
  normr := 0;
  posdef := true;
  norma := 0;
  For i:=1 To n Do
    Begin
      sumrowi := 0;
      For j:=1 To i Do
        sumrowi := sumrowi+abs(pal^[(i-1)*rwidth+j]);
      For j:=i+1 To n Do
        sumrowi := sumrowi+abs(pal^[(j-1)*rwidth+i]);
      If norma<sumrowi Then
        norma := sumrowi;
      t^[i] := 2*random-1;
     h := t^[i];
      h := abs(h);
      If normr<h Then
        normr := h
    End; {i}
  k := 0;
  while (k<n) and posdef Do
    Begin
      kmin1 := k;
     k := k+1;
      indk := (k-1)*rwidth;
      lkk := pal^[indk+k];
      For j:=1 To kmin1 Do
        lkk := lkk-sqr(pal^[indk+j]);
      If lkk <= 0 Then
        posdef := false
      Else
        Begin
          pal^[indk+k] := sqrt(lkk);
         lkk := pal^[indk+k];
          For i:=k+1 To n Do
            Begin
              indi := (i-1)*rwidth;
              h := pal^[indi+k];
              For j:=1 To kmin1 Do
                h := h-pal^[indk+j]*pal^[indi+j];
              pal^[indi+k] := h/lkk
            End; {i}
          h := t^[k];
          For j:=1 To kmin1 Do
            h := h-pal^[indk+j]*t^[j];
          t^[k] := h/lkk
        End {posdef}
    End; {k}
  If posdef Then
    Begin
      normt := 0;
      For i:=n Downto 1 Do
        Begin
          h := t^[i];
          For j:=i+1 To n Do
            h := h-pal^[(j-1)*rwidth+i]*t^[j];
          t^[i] := h/pal^[(i-1)*rwidth+i];
          h := abs(t^[i]);
          If normt<h Then
            normt := h
        End; {i}
      ca := norma*normt/normr
    End; {posdef}
  If posdef Then
    term := 1
  Else
    term := 2;
  freemem(t, sizeof(ArbFloat)*n);
End; {mdtgpd}

Procedure mdtgba(n, lb, rb, rwa: ArbInt; Var a: ArbFloat; rwl: ArbInt;
                 Var l:ArbFloat; Var p: ArbInt; Var ca: ArbFloat; Var term:ArbInt);

Var 
  sr, i, j, k, ipivot, m, lbj, lbi, ubi, ls,
             ii, jj, ll, s, js, jl, ubj       : ArbInt;
  ra, normr, sumrowi, pivot, normt, maxim, h  : ArbFloat;
          pl, au, sumrow, t, row              : ^arfloat1;
                                           pp : ^arint1;

Begin
  If (n<1) Or (lb<0) Or (rb<0) Or (lb>n-1) Or (rb>n-1) Or (rwl<0) Or (rwa<1) Then
    Begin
      term := 3;
     exit
    End; {term=3}
  sr := sizeof(ArbFloat);
  au := @a;
 pl := @l;
 pp := @p;
  ll := lb+rb+1;
  ls := ll*sr;
  getmem(sumrow, n*sr);
  getmem(t, n*sr);
  getmem(row, ls);
  lbi := n-rb+1;
 lbj := 0;
  jj := 1;
  For i:=lb Downto 1 Do
    Begin
      move(au^[i+jj], au^[jj], (ll-i)*sr);
      fillchar(au^[jj+ll-i], i*sr, 0);
      jj := jj+rwa
    End; {i}
  jj := (n-rb)*rwa+ll;
  For i:=1 To rb Do
    Begin
      fillchar(au^[jj], i*sr, 0);
      jj := jj+rwa-1
    End; {i}
  normr := 0;
 term := 1;
  ii := 1;
  For i:=1 To n Do
    Begin
      pp^[i] := i;
      sumrowi := omvn1v(au^[ii], ll);
      ii := ii+rwa;
      sumrow^[i] := sumrowi;
      h := 2*random-1;
     t^[i] := sumrowi*h;
      h := abs(h);
      If normr<h Then
        normr := h;
      If sumrowi=0 Then
        term := 4
    End; {i}
  ubi := lb;
 jj := 1;
  For k:=1 To n Do
    Begin
      maxim := 0;
     ipivot := k;
     ii := jj;
      If ubi<n Then
        ubi := ubi+1;
      For i:=k To ubi Do
        Begin
          sumrowi := sumrow^[i];
          If sumrowi <> 0 Then
            Begin
              h := abs(au^[ii])/sumrowi;
              ii := ii+rwa;
              If maxim<h Then
                Begin
                  maxim := h;
                 ipivot := i
                End {maxim<h}
            End {sumrowi <> 0}
        End; {i}
      If maxim=0 Then
        Begin
          lbj := 1;
         ubj := ubi-k;
          For j:=lbj To ubj Do
            pl^[(k-1)*rwl+j] := 0;
          For i:=k+1 To ubi Do
            Begin
              ii := (i-1)*rwa;
              For j:=2 To ll Do
                au^[ii+j-1] := au^[ii+j];
              au^[ii+ll] := 0
            End; {i}
          term := 4
        End {maxim=0}
      Else
        Begin
          If ipivot <> k Then
            Begin
              ii := (ipivot-1)*rwa+1;
              move(au^[ii], row^, ls);
              move(au^[jj], au^[ii], ls);
              move(row^, au^[jj], ls);
              h := t^[ipivot];
             t^[ipivot] := t^[k];
             t^[k] := h;
              pp^[k] := ipivot;
              sumrow^[ipivot] := sumrow^[k]
            End; {ipivot <> k}
          pivot := au^[jj];
         jl := 0;
         ii := jj;
          For i:=k+1 To ubi Do
            Begin
              jl := jl+1;
              ii := ii+rwa;
              h := au^[ii]/pivot;
              pl^[(k-1)*rwl+jl] := h;
              For j:=0 To ll-2 Do
                au^[ii+j] := au^[ii+j+1]-h*au^[jj+j+1];
              au^[ii+ll-1] := 0;
              If term=1 Then
                t^[i] := t^[i]-h*t^[k]
            End {i}
        End; {maxim <> 0}
      jj := jj+rwa
    End; {k}
  If term=1 Then
    Begin
      normt := 0;
     ubj := -lb-1;
      jj := n*rwa+1;
      For i:=n Downto 1 Do
        Begin
          jj := jj-rwa;
          If ubj<rb Then
            ubj := ubj+1;
          h := t^[i];
          For j:=1 To ubj+lb Do
            h := h-au^[jj+j]*t^[i+j];
          t^[i] := h/au^[jj];
          h := abs(t^[i]);
          If normt<h Then
            normt := h
        End; {i}
      ca := normt/normr
    End {term=1}
  Else
   ca := giant;
  freemem(sumrow, n*sr);
  freemem(t, n*sr);
  freemem(row, ls)
End; {mdtgba}

Procedure mdtgpb(n, lb, rwidth: ArbInt; Var al, ca: ArbFloat;
                 Var term: ArbInt);

Var 
    posdef                                           : boolean;
    i, j, k, r, p, q, ll, llmin1, jmin1, indi        : ArbInt;
    h, normr, normt, sumrowi, alim, norma            : ArbFloat;
    pal, t                                           : ^arfloat1;

    Procedure decomp(i, r: ArbInt);

    Var 
        k, ii, ir : ArbInt;
    Begin
      ii := (i-1)*rwidth;
      ir := (r-1)*rwidth;
      h := pal^[ii+j];
     q := ll-j+p;
      For k:=p To jmin1 Do
        Begin
          h := h-pal^[ii+k]*pal^[ir+q];
         q := q+1
        End; {k}
      If j<ll Then
        pal^[ii+j] := h/pal^[ir+ll]
    End; {decomp}

    Procedure lmin1t(i: ArbInt);

    Var 
        k:ArbInt;
    Begin
      h := t^[i];
     q := i;
      For k:=llmin1 Downto p Do
        Begin
          q := q-1;
         h := h-pal^[indi+k]*t^[q]
        End; {k}
      t^[i] := h/alim
    End; {lmin1t}

Begin
  If (lb<0) Or (lb>n-1) Or (n<1) Or (rwidth<1) Then
    Begin
      term := 3;
     exit
    End; {wrong input}
  pal := @al;
  getmem(t, n*sizeof(ArbFloat));
  ll := lb+1;
 normr := 0;
 p := ll+1;
 norma := 0;
  For i:=1 To n Do
    Begin
      If p>1 Then
        p := p-1;
      indi := (i-1)*rwidth+p;
      sumrowi := omvn1v(pal^[indi], ll-p+1);
      r := i;
     j := ll;
      while (r<n) and (j>1) Do
        Begin
          r := r+1;
         j := j-1;
          sumrowi := sumrowi+abs(pal^[(r-1)*rwidth+j])
        End; {r,j}
      If norma<sumrowi Then
        norma := sumrowi;
      h := 2*random-1;
     t^[i] := h;
      h := abs(h);
      If normr<h Then
        normr := h
    End; {i}
    llmin1 := ll-1;
    p := ll+1;
    i := 0;
    posdef := true ;
    while (i<n) and posdef Do
      Begin
        i := i+1;
        indi := (i-1)*rwidth;
        If p>1 Then
          p := p-1;
        r := i-ll+p;
       j := p-1;
        while j<llmin1 Do
          Begin
            jmin1 := j;
           j := j+1;
            decomp(i, r);
           r := r+1
          End; {j}
        jmin1 := llmin1;
       j := ll;
       decomp(i, i);
        If h <= 0 Then
          posdef := false
        Else
          Begin
            alim := sqrt(h);
           pal^[indi+ll] := alim;
            lmin1t(i)
          End
      End; {i}
    If posdef Then
      Begin
        normt := 0;
       p := ll+1;
        For i:=n Downto 1 Do
          Begin
            If p>1 Then
              p := p-1;
            q := i;
           h := t^[i];
            For k:=llmin1 Downto p Do
              Begin
                q := q+1;
               h := h-pal^[(q-1)*rwidth+k]*t^[q]
              End; {k}
            t^[i] := h/pal^[(i-1)*rwidth+ll];
            h := abs(t^[i]);
            If normt<h Then
              normt := h
          End; {i}
        ca := norma*normt/normr
      End; {posdef}
    If posdef Then
      term := 1
    Else
      term := 2;
  freemem(t, n*sizeof(ArbFloat));
End; {mdtgpb}

Procedure mdtdtr(n: ArbInt; Var l, d, u, l1, d1, u1: ArbFloat;
                 Var term:ArbInt);

Var
                      i, j, s : ArbInt;
                       lj, di : ArbFloat;
             pd, pu, pd1, pu1 : ^arfloat1;
                      pl, pl1 : ^arfloat2;

Begin
  If n<1 Then
    Begin
      term := 3;
     exit
    End; {wrong input}
  pl := @l;
  pd := @d;
  pu := @u;
  pl1 := @l1;
  pd1 := @d1;
  pu1 := @u1;
  s := sizeof(ArbFloat);
  move(pl^, pl1^, (n-1)*s);
  move(pd^, pd1^, n*s);
  move(pu^, pu1^, (n-1)*s);
  j := 1;
  di := pd1^[j];
  If di=0 Then
    term := 2
  Else
    term := 1;
  while (term=1) and (j <> n) Do
    Begin
     i := j;
     j := j+1;
     lj := pl1^[j]/di;
     pl1^[j] := lj;
     di := pd1^[j]-lj*pu1^[i];
     pd1^[j] := di;
     If di=0 Then
      term := 2
    End {j}
End; {mdtdtr}

Begin
  randseed := 12345
End.

{
  $Log$
  Revision 1.1  2000-01-24 22:08:58  marco
   * initial version


}