
// private replacement of buggy strutils.isEmptyStr()
{}function mtystr(s:string):boolean;
var i, n: longint;
begin
  mtystr := false;
  n := length(s);
  for i := 1 to n do begin
    if not ( s[i] in BLANK ) then exit;
  end;
  mtystr := true
end;

// Set of functions to write a matrix to a text file
// UGLY KLUDGE to make a pretty generic double output format
// the public global var ugly enables to use this or not...
{}procedure writx(var f: text; x: double);
const
  r=4.99E-7;
  wd=0; // using tab as delimiter so fieldwidth 0 is ok
  sf=6; // was 5 check that it makes no problem
var
  y:double; p,pp,de,w,dc:longint;
  intg,n,ne:boolean;
begin
  if (abs(x)<=maxint) then intg:=abs(x-trunc(x))=0 else intg:=false;
  if intg then write(f,trunc(x):wd) else begin
    n := x<0;
    if n then y:=-x else y:=x;
    p:=0;
    while (y>=10) do begin y:=y/10; p:=succ(p) end;
    while (y<1)   do begin y:=y*10; p:=pred(p) end;
    if (p=0) then write(f,x:wd:sf-1) else begin
      ne:=p<0;
      if ne then pp:=-p else pp:=p;
      de:=2;
      while (pp>9) do begin pp:=pp div 10; de:=succ(de) end;
      if ne then de:=succ(de);
      y:=y+double(r);
      w:=wd-de;
      if (p>-(sf-2)) and (p<sf) then begin
        dc:=pred(sf-p);
        if (dc<1) then dc:=1;
        write(f,x:wd:dc) // TODO check if sf=6 make sense ?
      end else
        if n then write(f,-y:w:pred(sf),'E',p:1)
             else write(f,y:w:pred(sf),'E',p:1)
    end
  end
end;

{}procedure writx(x:double); inline;
begin
  writx(stdout,x);
end;

{ ==============================================================================
  3 variants mimic the most used prototypes of mtx functions
}
{p}procedure wtrace(m:mtx; s:string);
begin
  if not ptrace then exit;
  writeln(output,'# ',s, ' [',m.lgn,',',m.col,']');
  flush(output);
  mwr(m,output);
end;

{}procedure wtrace(a,b:mtx; s:string);
begin
  if not ptrace then exit;
  writeln(output,'# ',s,' [',a.lgn,',',a.col,'] [',b.lgn,',',b.col,']');
  flush(output);
  mwr(a,output);
end;

{}procedure wtrace(r,a,b:mtx; s:string);
begin
  if not ptrace then exit;
  writeln(output,'# ',s,' [',r.lgn,',',r.col,'] [',a.lgn,',',a.col,'] [',b.lgn,',',b.col,']');
  flush(output);
  mwr(r,output);
end;

//==============================================================================
//
//                     Dynamic memory management
//     a zoo of mtx function to perform allocation and initialisation
//
//==============================================================================

{
NOTE :
The code relies on this properties of the memory allocation function
described in fpc documentation:

** function ReAllocMem(var p: pointer;Size: PtrInt) : pointer;

Description: ReAllocMem resizes the memory pointed to by P so it has
size Size. The value of Pmay change during this operation. The contents
of the memory pointed to by P (if any) will be copied to the new
location, but may be truncated if the newly allocated memory block is
smaller in size. If a larger block is allocated, only the used memory is
initialized, extra memory will not be zeroed out. Note that P may be nil,
in that case the behaviour of ReAllocMem is equivalent to Getmem.
}

// Realloc mem. if needed and possible
{}function mem(var r:mtx; m:longint): boolean;
var p: Pdouble;
begin
  mem := false;
  //writeln(stderr,'**** in mem (',r.nam,') **** m=',m,'  r.e=', longint(r.e));
  assert(m>=0);
  if r.e = NIL then begin // WARNING FAILS WITH LOCAL MTX !!!
    //writeln(stderr,'***** mem ***** First alloc. *****');
    r.nam := '';
    r.lgn := 0;
    r.col := 0;
    r.t := false;
    r.dl := 0;
    r.dc := 0;
    r.max := m;
    r.dyn := true; // TODO check it should be reset in all smx functions
    r.smx := false;
    p := ReAllocMem(r.e, (m+1)*sizeof(double));
    assert( p <> NIL );
    r.e := p;
{ ifdef gsl} // create records for ***GSL*** interface
    new(r.mat);
    new(r.mat^.block);
    // TODO check this...
    with r.mat^ do begin
      size1 := 0;
      size2 := 0;
      tda := 0;
      data := @r.e[1];
      owner := false;
      block^.size := m;
      block^.data := data
    end;
    new(r.vec);
    new(r.vec^.block);
    with r.vec^ do begin
      size := 0;
      stride := 1;
      data := @r.e[1];
      owner := false;
      block^.size := m;
      block^.data := data
    end;
{ endif} // gsl
  end else assert( not r.smx ); // I guess this make sense YES cf. zer one eye
  if ( r.max < m ) then begin
    //writeln(stderr,'***** mem ***** Realloc. *****');
    assert( r.dyn );
    p := ReAllocMem(r.e, (m+1)*sizeof(double));
    assert( p <> NIL );
    r.max := m;
    r.e := p;
{ ifdef gsl}// ***GSL*** update ALL data pointers in vec^ and mat^
    with r.mat^ do begin
      data := @p[1];
      block^.data := data;
      block^.size := m
    end;
    with r.vec^ do begin
      data := @p[1];
      block^.data := data;
      block^.size := m
    end
{ endif}
  end;
  //writeln(stderr,'***** exit mem (',r.nam,') ***** r.e =', longint(r.e) );
  mem := true
end;

// return max. vector size TODO inline ? as it is private ?
{}function mem(r:mtx): longint;
begin
  mem := r.max
end;

{}procedure dim(var r:mtx; lgn, col:longint);
begin
// ***GSL*** alloc and data fields init in mem
  assert( mem(r,lgn*col) );
  assert( not r.smx );
  r.lgn := lgn;
  r.col := col;
  r.smx := false;
  r.dl := 1;
  r.dc := col;
  r.t := false;
// ***GSL***
{ ifdef gsl}
  with r.mat^ do begin
    size1 := lgn;
    size2 := col;
    tda := col;
  end;
  with r.vec^ do begin
    size := lgn*col;
    stride := 1
  end;
{ endif}
end;

{
   dim and fill with a value : used by zer one...
}
{}procedure dim(var r:mtx; lgn, col:longint; x:double); inline;
var k: longint;
begin
  dim(r,lgn,col);
  for k := 1 to lgn*col do r.e[k] := x
end;

{}procedure chk(var r:mtx; a, b:mtx);
begin
  if (a.lgn<>b.lgn) or (a.col<>b.col) then exit;
  if r.e = NIL then r.smx := false;
  if (r.smx) then
    // so GSL fields pointers are NIL
    assert( (r.lgn=a.lgn) and (r.col=a.col) )
  else
    // r is plain and dim ( calling mem ) must update ***GSL*** fields...
    dim(r,a.lgn,a.col)
end;

{}procedure chk(var r:mtx; a, b:mtx; var bdl, bdc: longint);
begin
  bdl := b.dl;
  bdc := b.dc;
  if (b.col=1) and (a.lgn=b.lgn) then
    bdl := 0
  else if (b.lgn=1) and (a.col=b.col) then
    bdc := 0
  else if (b.lgn=1) and (b.col=1) then begin
    bdc := 0;
    bdl := 0
  end else
    assert( (a.col=b.col) and (a.lgn=b.lgn) );
  if r.e = NIL then r.smx := false;
  if (r.smx) then // so far slices not mapped to GSL stuff...
    assert(  (r.lgn=a.lgn) and (r.col=a.col) )
  else
    dim(r,a.lgn,a.col)
end;

{
  same as equ(r,a) but NOT traced
}
{}procedure cpy(var r:mtx; a:mtx);
var i, j, lgn, col, ia, ir, ia1, ir1: longint;
begin
  chk(r,a,a);
  lgn := a.lgn;
  col := a.col;
  ia1 := 1;
  ir1 := 1;
  for i := 1 to lgn do begin
    ia := ia1;
    ir := ir1;
    for j := 1 to col do begin
       r.e[ir] := a.e[ia];
       ia := ia + a.dl;
       ir := ir + r.dl
    end;
    ia1 := ia1 + a.dc;
    ir1 := ir1 + r.dc
  end
end;

//=============================================================================
//             read mtx from text file with autom. mem alloc.
//=============================================================================

// Corrected a pb when reading anonymous mtx after named ones...
// mnam must be in the arg. list... and a named mtx skipped if searching a noname
// readln until it finds either a letter in col 1 ( first word is name of matrix )
// or a number (noname mtx i.e. name = '' )
// returns true if found a match, nam is the name of the mtx or '' for a plain xy table
// n is the number of words in the first data line of the mtx/table s is this string
// which can be processed to retrieve the n numerical values
{p}function read1st(mnam:string; var f:text; var s: string; var nam: string; var n:longint):boolean;
var s0: string; ok: boolean; p:longint; x:double; code:integer;
begin
  read1st := false;
  ok := false;
  repeat
    readln(f,s);
    if mtystr(s) then begin
       if eof(f) then exit else continue;
    end;
    if pos('&',s)>0 then exit; // end of data needed for demo code
    if pos('#',s)>0 then continue;
    s0 := ExtractWordPos(1,s+' x',WordDelim,p); // hope this func. is ok ?
 // if mnam is neither '' nor '*' we can skip val() and set code=-1
    if (mnam='') or (mnam='*') then val(s0,x,code) else code:=-1;
    if  code=0 then begin
      // we assume that s is the first line of a plain noname xy table
      // and s0 the first number of the table
      nam := '';
      // *********** OK only if nam agree with mnam
      ok := (nam = mnam) or (mnam='*'); // ************
      if not ok then while not mtystr(s) do readln(f,s)
    end else begin
      // s0 must be a name
      nam := s0;
      ok := (nam = mnam) or (mnam='*');
      if ok then readln(f,s) else while not mtystr(s) do readln(f,s)
      // s must be the first data line of mtx mnam
    end;
  until ok;
  // count the number of item in this first line
  // TODO achtung n returned by WordCount() is integer !!!!!!!!
  n := WordCount(s,WordDelim);
  // assert all words are valid numbers...
  read1st := n>0;
end;

// private tool used by next lgn() avatar
{p}function inrois(x:double; r:mtx):boolean;
var i: longint;
begin
  assert( col(r)=2 );
  inrois:= false;
  for i := 1 to lgn(r) do begin
    if (x>elm(r,i,1)) and (x<=elm(r,i,2)) then begin
      inrois := true;
      break
    end
  end
end;

// If "r" is a submatrix ( field smx true ? TODO quid de trp() ) then
// returns true if product r = a*b is ok
// if "r" is a plain matrix
// redim and realloc dyn. mem if needed...
{p}procedure chkprd(var r:mtx; a, b:mtx); inline;
begin
  assert( (a.col=b.lgn) );
  if r.e = NIL then r.smx := false;
  if (r.smx) then
    assert( (r.lgn=a.lgn) and (r.col=b.col) )
  else
    dim(r,a.lgn,b.col)
end;

//=============================================================================
//            Gauss/Jordan max. pivot matrix inversion functions
//=============================================================================
// TODO MAY BE replace by a legacy gsl or netlib code ???
// or use SVD codes if required ( ill cond. mtx ).
// may be borrow some scilab library fortran code to do this...
// NOTE :
// we dont care "t" flag because inv(A') = inv(A)' and det(a) = det(a')
// so in situ inversion is made: if t then we leave it like this
// so the inverse is also affected by the transpose flag
// legacy Gauss/Jordan max pivot in situ inversion...
{}function matinv(a:mtx; var  d: double): boolean;
const
  eps = EpsDouble; // TODO check for a good value ?
var
  n, i, j, k, l: longint;
  dd, amax: extended;
  save: double;
  ik, jk: ^longint;
  //t: boolean;
begin
  matinv := false;
// only plain matrices are usable not a slice.
// Note that inverse of transpose is transpose of inverse so
// we can ignore transposition alltogether keping dl and dc untouched
  assert( not a.smx );
  dd := 1.0;
  n := a.lgn;
  assert(a.col=n);
  GetMem(ik,(n+1)*sizeof(longint));
  GetMem(jk,(n+1)*sizeof(longint));
  for k := 1 to n do begin // k=1 to n
    amax := 0.0;
    for i := k to n do
      for j := k to n do begin
        save := a.e[(i-1)*n+j];
        if abs(save)>abs(amax) then begin
          amax := save; ik[k] := i; jk[k] := j
        end
      end;
    if abs(amax)<eps then begin
      d := 0;
      FreeMem(jk);
      FreeMem(ik);
      exit;
    end;
    i := ik[k];
    if i<>k then for j := 1 to n do begin
      save := a.e[(k-1)*n+j];
      a.e[(k-1)*n+j]:=a.e[(i-1)*n+j];
      a.e[(i-1)*n+j] := -save
    end;
    j := jk[k];
    if j<>k then for i := 1 to n do begin
      save := a.e[(i-1)*n+k];
      a.e[(i-1)*n+k]:=a.e[(i-1)*n+j];
      a.e[(i-1)*n+j]:=-save
    end;
    for i := 1 to n do
      if i<>k then
        a.e[(i-1)*n+k] := -a.e[(i-1)*n+k]/amax;
    for i := 1 to n do
      for j := 1 to n do
        if (i<>k)and(j<>k)then
          a.e[(i-1)*n+j] := a.e[(i-1)*n+j]+a.e[(i-1)*n+k]*a.e[(k-1)*n+j];
    for j := 1 to n do
      if j<>k then a.e[(k-1)*n+j]:=a.e[(k-1)*n+j]/amax;
    a.e[(k-1)*n+k] := 1.0/amax;
    dd := dd*amax
  end; // k=1 to n
  for l := 1 to n do begin // l=1 to n
    k := n-l+1;
    j := ik[k];
    if j>k then for i := 1 to n do begin
      save := a.e[(i-1)*n+k];
      a.e[(i-1)*n+k] := -a.e[(i-1)*n+j];
      a.e[(i-1)*n+j] := save
    end;
    i := jk[k];
    if i>k then for j := 1 to n do begin
      save := a.e[(k-1)*n+j];
      a.e[(k-1)*n+j] := -a.e[(i-1)*n+j];
      a.e[(i-1)*n+j] := save
    end
  end; // l=1 to n
  FreeMem(jk);
  FreeMem(ik);
  d := double(dd);
  wtrace(a,'(mat)inv');
  matinv := true
end;

//=============================================================================
//  interface with netlib fmm's spline interpolation functions: spline & seval
//                       http://www.netlib.org/fmm
//=============================================================================
// Prototypes of fortran routines spline_ and seval_
// source spline.f; compile with:  gfortran -c spline.f
// { linklib j}   // link to object in libj.a or if not installed in a lib
// {$L spline.o}  // link to local spline.o
{p}procedure spline_(var n:longint; x, y, b, c, d:pointer); cdecl; external;
{p}function seval_(var n:longint; var u:double;  x, y, b, c, d:pointer):double; cdecl; external;

// return index of interpolation segment
// f=f(u) fp=f'(u)
{p}function seval1_(var n:longint; var u:double;  x, y, b, c, d:pointer; var f, fp:double):longint; cdecl; external;

//=============================================================================
// ancillary code to implement roi rex functions and lgn() avatars
//=============================================================================

// return index of an element in the internal
// mtx.e array of double ( element 0 is allocated but
// mostly ignored for pascal/f77 like indexing over
// C style array double*
function idx(a:mtx; i, j: longint):longint;
begin
  assert( (i>0) and (i<=a.lgn) );
  assert( (j>0) and (j<=a.col) );
  idx := 1 + pred(i)*a.dc + pred(j)*a.dl;
end;

// dichotomic search of index of first xi >= u starting from x1
// This function seems to behave properly afaik
function fnd(a:mtx; u:double):longint;
var
  n: longint;
  i, j, k: longint;
  x1, xn, xi, xj, xk: double;
begin
  if mty(a) then exit;
  n := lgn(a);
  x1 := elm(a,1,1);
  xn := elm(a,n,1);
  if u<=x1 then begin
    i := 0; //
    j := 1  // all >= u
  end else if u>xn then begin
    i := n;  //
    j := n+1 // none >= u ( return n+1 TODO should it be 0 TODO )
  end else begin
    i := 1;
    j := n;
    xi := x1;
    xj := xn;
    while (j-i)>1 do begin
      k := (i+j) div 2;
      xk := elm(a,k,1);
      if u <= xk then begin
        j := k;
        xj := xk
      end else begin
        i := k;
        xi := xk
      end;
      assert(u>xi);
      assert(xj>=u)
    end
  end;
  assert((j-i)=1);
  fnd := j
end;

