
{$H+}          // ansi strings
{$C+}          // assertions on
{$inline on}   // used by opr() add family functions and other "wrapper"
{$L spline.o}  // link to local spline.o
{ linklib j}   // link to /usr/local/lib/libj.a

{
  @abstract(Unit to manage abstract matrix data in a Scilab/Matlab like way)

  TODO use html code here and link to Ooo presentation in html or pdf
  Introduction

  This unit provides a matrix data type and functions to perform input/output
  and most simple matrix operations. * link to concept list *

  The main goal is to achieve the compacity of scilab/matlab like code
  transparent dynamic memory management; simple unified IO NO messy
  formatting... just a pretty and ugly format ! Almost NO explicit loops
  through transparent dimension management;

  Some definitions: plain mtx versus slice

  function and procedure overloading

  exceptions handling with assert (or returned value... a bit messy)

  Understanding softdimensions dim. expansion and reduction by operator

  For advanced users, some hooks are provided to use legacy
  scientific libraries such as fortran subroutines from netlib or the
  numerical constants, vector and matrix from the Gnu Scientific library.

  TODO add pointers to urls with presentations and typical examples or
  tutorials and also memo for advanced usage or notes for developers

  TODO see also similar projects such as eigen matpack the GSL and where this
  idea comes from ( technicalities from Braquelaire C programming... the old
  book about abstract data types when object programming was not yet born ! )

  What are the main advantages over other softwares and limitations.
  Some words about R approach to i/o i.e. NO IO just keep everething together
  from one session to another in a given context/folder... persistence or
  simple I/O versus the spreadsheet way...

  Simple IO easy interaction with shell scripts of all kinds and spreadsheets*
  as cut/paste works fine.

  dark corners... potential memory leaks why use/not use the the fpc option
  to check memory leaks. The message: local mtx should only be assigned slices
  to avoid memory alloc in functions or procedures ( which will NOT be
  deallocated. One may use private mtx in implementation of user's units using
  mtx09.

  references: abstract programming Braquelaire's book GSL manual
  may be pointers to wikipedia function/procedure overloading and or the
  official freePascal documentation.

}
unit mtx09;

interface

{$ifdef MATH}
   uses
      strutils,
      math;
{$else}
   uses
      strutils;
   const
      Infinity = 1e300;
{$endif}

const
  ElmDelim=[' ',#9,',']; //< element delimiter within an ascii representation of a matrix
  RowDelim=[';',#10];    //< row delimiter
  EpsDouble = 1e-30;     //< TODO use HSL or GSL funct. or cst.
  WordDelim=[#0..' '];   //< assert ascii char set TODO check where used
  BLANK=[#32,#9];        //< spaces and tabs
  TAB=#9;                //< ascii tab character used as delim in output

type
  C_int = LongBool;   //< 4 bytes to match C int type
  size_t = longint;   //< TODO check if this is the best C<>fpc match
  Pchar  = ^char;     //< this type may be declared elsewhere
  Pdouble  = ^double; //< this type is already declared in fpc...
  // opr_t = procedure(var b:double; a:double); //< C style b += a


  //  interface with the GSL standard (double) vector and matrix types
  block=record
    size : size_t;  //< number of elements in the block
    data : ^double; //< pointer to the first element of the block
  end;
  Pblock = ^block;

  // interface with the GSL vector of doubles type
  vector = record
    size: size_t;
    stride: size_t;
    data: Pdouble;
    block: Pblock;
    owner: C_int;
  end;
  Pvector = ^vector;

  // interface with the GSL; standard matrix of double
  matrix = record
    size1: size_t;
    size2: size_t;
    tda: size_t;
    data: Pdouble;
    block: Pblock;
    owner: C_int;
  end;
  Pmatrix = ^matrix;

  {
      @abstract( mtx: abstract matrix data type )
      TODO a pointer to the relevant discussion in the GSL and elsewhere
      and explain in detail the data structure

      All elements of the matrix are stored in a memory area pointed to by
      mtx.e; (e[0] is not used as Fortran and Scilab math indexing from
      1 to n rather than C style) e[1] is the first element of the matrix i.e.
      (1,1). The fields e.dl and e.dc are the offset to add to the index of
      e[i] to access the next element within the same row or column.
      This enables to use several mtx records pointing to the same memory area
      but with different indexing schemes. For instance a transposed matrix
      can be done by swaping dl and dc values. This simple setup enables to
      handle compact sub matrices in a fairly efficient and simple way.
      See for instance the documentation of the GNU Scientific Library.
  }
  mtx = record
    smx: boolean; //< slice flag: @true if the record is a slice; @false plain matrix.
    nam: string; //< matrix name; may be an empty string
    lgn: longint; //< number of rows (ligne in french)
    col: longint; //< number of columns
    max: longint; //< maximum number within the current allocated block
    dl: longint; //< index offset for adjacent elements of the same row
    dc: longint; //< index offset for adjacent elements of the same column
    t: boolean; //< transposition flag
    e: ^double; //< pointer to the top left (first) element of the matrix
    dyn: boolean; //< true if matrix is stored in dynamic memory and can be reallocated
    mat: Pmatrix; //< stubb to use GSL matrix of doubles
    vec: Pvector; //< stubb to use GSL vectors of doubles
  end;

// include automatically build interface of functions and procedures
{$include mtx09.inc}

implementation

const
  addnullelm = 0.0; //< see opr: target in Makefile
  subnullelm = 0.0; //< !! DO NOT change names of these constants DO NOT !!
  mulnullelm = 1.0; //< see awk script generating sub mul dvd on add model
  dvdnullelm = 1.0; //< needed in scilab clone add(r,a,true)...

var
  x11: mtx;                //< manage double mtx casting
  cat0: mtx;               //< global mtx workspace of cat/stk
  elm0: mtx;               //< worspace of elm(mtx,string) avatar
  ptrace: boolean = false; //< wrt() at end of functions
  pretty: boolean=true;    //< pretty (smart) or ugly exponential output
//  format: integer = 6;   // TODO choose output format pretty (default)
                           // or full exp. format. see mwr() and writx()
  apdmode: boolean=true;   // mwr(a,'file') will append by default

{$include private.inc} // private funct./proc. NOT IN INTERFACE

//==============================================================================
//     avatar of reset rewrite (akin to those in legacy turbo pascal)
//==============================================================================

// TODO fmt(n) to control width of pretty output

{**
@abstract(@html(
<pre><span style="font-weight: bold; color: rgb(0, 0, 0);">procedure</span><span
 style="color: rgb(0, 0, 0);"> fmt(d:</span><span
 style="color: rgb(128, 0, 0);">longint</span><span
 style="color: rgb(0, 0, 0);">); </span><span
 style="font-weight: bold; color: rgb(0, 0, 128);">inline</span><span
 style="color: rgb(0, 0, 0);">;</span>
</pre>
Control output format of numbers: pretty or ugly (but more accurate)
<dl>
  <dt>d=1</dt>
  <dd>pretty output (default),</dd>
  <dt>d=0</dt>
  <dd>ugly but accurate output.</dd>
</dl>
<a href="ex/fmt_.pas.html">Example</a>
))}
procedure fmt(d:longint); inline;
begin
  pretty := d>0
end;

{**
@abstract(@html(
<pre><span style="font-weight: bold; color: rgb(0, 0, 0);">function</span><span
 style="color: rgb(0, 0, 0);"> reset(</span><span
 style="font-weight: bold; color: rgb(0, 0, 0);">var</span><span
 style="color: rgb(0, 0, 0);"> f: </span><span
 style="color: rgb(128, 0, 0);">text</span><span
 style="color: rgb(0, 0, 0);">; s:</span><span
 style="color: rgb(128, 0, 0);">string</span><span
 style="color: rgb(0, 0, 0);">): </span><span
 style="color: rgb(128, 0, 0);">boolean</span><span
 style="color: rgb(0, 0, 0);">; </span><span
 style="font-weight: bold; color: rgb(0, 0, 128);">inline</span><span
 style="color: rgb(0, 0, 0);">;</span>
</pre>
  Open text stream "f" to read from file named "s".<br>
  Returns false if it fails.<br>
  See also opn().<br>
  <a href="ex/stream.html">Example</a>
))}
function reset(var f: text; s:string):boolean; inline;
begin
 {$I-}
 assign(f,s);
 system.reset(f);
 reset := ioresult=0
 {$I+}
end;

{**
@abstract(@html(
<pre><span style="font-weight: bold; color: rgb(0, 0, 0);">procedure</span><span
 style="color: rgb(0, 0, 0);"> opn(s:</span><span
 style="color: rgb(128, 0, 0);">string</span><span
 style="color: rgb(0, 0, 0);">; </span><span
 style="font-weight: bold; color: rgb(0, 0, 0);">var</span><span
 style="color: rgb(0, 0, 0);"> f:</span><span
 style="color: rgb(128, 0, 0);">text</span><span
 style="color: rgb(0, 0, 0);">); </span><span
 style="font-weight: bold; color: rgb(0, 0, 128);">inline</span><span
 style="color: rgb(0, 0, 0);">;</span>
</pre>
Open text stream "f" to read from file named "s".
Raises an error if it fails.<br>
NOTE: see also reset()<br>
<a href="ex/stream0.pas.html">Example</a>
))}
procedure opn(s:string; var f:text); inline;
begin
  assign(f,s);
  system.reset(f);
end;

{**
@abstract(@html(
<pre><span style="font-weight: bold; color: rgb(0, 0, 0);">procedure</span><span
 style="color: rgb(0, 0, 0);"> apn(s:</span><span
 style="color: rgb(128, 0, 0);">string</span><span
 style="color: rgb(0, 0, 0);">; </span><span
 style="font-weight: bold; color: rgb(0, 0, 0);">var</span><span
 style="color: rgb(0, 0, 0);"> f:</span><span
 style="color: rgb(128, 0, 0);">text</span><span
 style="color: rgb(0, 0, 0);">); </span><span
 style="font-weight: bold; color: rgb(0, 0, 128);">inline</span><span
 style="color: rgb(0, 0, 0);">;</span>
</pre>
  Open text stream "f" to append to file named "s".<br>
  See also: rwr() opn() clo()<br>
  <a href="ex/stream0.pas.html">Example</a>
))}
procedure apn(s:string; var f:text); inline;
begin
  append(f,s)
end;

{**
@abstract(@html(
Close text stream "f". An alias of system.close().
Raises an error if it fails.<br>
See also: opn() rwr() mrd() mwr()<br>
<a href="ex/stream0.pas.html">Example</a>
))}
procedure clo(var f:text); inline;
begin
  system.close(f)
end;

{**
@abstract(@html(
<pre><span style="font-weight: bold; color: rgb(0, 0, 0);">procedure</span><span
 style="color: rgb(0, 0, 0);"> reset(</span><span
 style="font-weight: bold; color: rgb(0, 0, 0);">var</span><span
 style="color: rgb(0, 0, 0);"> f: </span><span
 style="color: rgb(128, 0, 0);">text</span><span
 style="color: rgb(0, 0, 0);">); </span><span
 style="font-weight: bold; color: rgb(0, 0, 128);">inline</span><span
 style="color: rgb(0, 0, 0);">;</span>
</pre>
Pristine freepascal reset: the freepascal standard
way to open text streams:<br>
<a href="ex/stream2.pas.html">Example</a>
))}
procedure reset(var f: text); inline;
begin
  system.reset(f)
end;

{**
@abstract(@html(
<pre><span style="font-weight: bold; color: rgb(0, 0, 0);">function</span><span
 style="color: rgb(0, 0, 0);"> rewrite(</span><span
 style="font-weight: bold; color: rgb(0, 0, 0);">var</span><span
 style="color: rgb(0, 0, 0);"> f: </span><span
 style="color: rgb(128, 0, 0);">text</span><span
 style="color: rgb(0, 0, 0);">; s:</span><span
 style="color: rgb(128, 0, 0);">string</span><span
 style="color: rgb(0, 0, 0);">): </span><span
 style="color: rgb(128, 0, 0);">boolean</span><span
 style="color: rgb(0, 0, 0);">; </span><span
 style="font-weight: bold; color: rgb(0, 0, 128);">inline</span><span
 style="color: rgb(0, 0, 0);">;</span>
</pre>
  Open text stream "f" for writing to file "s".
  Return false if rewrite failed.
  One may wrap the call with assert() for a simple exception management:
  <tt>assert(rewrite(fo,onam))</tt><br>
  <a href="ex/stream.html">Example</a>
))}
function rewrite(var f: text; s:string): boolean; inline;
begin
  {$I-}
  assign(f,s);
  system.rewrite(f);
  rewrite := ioresult=0
  {$I+}
end;

{**
@abstract(@html(
<pre><span style="font-weight: bold; color: rgb(0, 0, 0);">procedure</span><span
 style="color: rgb(0, 0, 0);"> rewrite(</span><span
 style="font-weight: bold; color: rgb(0, 0, 0);">var</span><span
 style="color: rgb(0, 0, 0);"> f: </span><span
 style="color: rgb(128, 0, 0);">text</span><span
 style="color: rgb(0, 0, 0);">); </span><span
 style="font-weight: bold; color: rgb(0, 0, 128);">inline</span><span
 style="color: rgb(0, 0, 0);">;</span>
</pre>
Pristine freepascal rewrite.<br>
The freepascal standard way tp open text streams:<br>
<a href="ex/stream2.pas.html">Example</a>
))}
procedure rewrite(var f: text); inline;
begin
  system.rewrite(f)
end;

{**
@abstract(@html(
<pre><span style="font-weight: bold; color: rgb(0, 0, 0);">function</span><span
 style="color: rgb(0, 0, 0);"> append(</span><span
 style="font-weight: bold; color: rgb(0, 0, 0);">var</span><span
 style="color: rgb(0, 0, 0);"> f:</span><span
 style="color: rgb(128, 0, 0);">text</span><span
 style="color: rgb(0, 0, 0);">; s:</span><span
 style="color: rgb(128, 0, 0);">string</span><span
 style="color: rgb(0, 0, 0);">):</span><span
 style="color: rgb(128, 0, 0);">boolean</span><span
 style="color: rgb(0, 0, 0);">;</span>
</pre>
   This function enables to append to a file with a single instruction.
   It doesn't throw an error if the file named "s"
   doesn't exist unlike pristine freePascal append()
   but does a rewrite instead.<br>
   <a href="ex/append1.pas.html">Example</a><br>
   <a href="ex/append2.pas.html">See also freePascal append</a>
))}
function append(var f:text; s:string):boolean;
begin
  append := false;
  if s='' then begin
    f := stdout;
    append := true;
    exit;
  end;
  assign(f,s);
{$I-}
  system.append(f);
  if ioresult<>0 then system.rewrite(f);
  if ioresult<>0 then exit;
{$I+}
  append := true;
end;

{**
  @abstract(@html(
<pre><span style="font-weight: bold; color: rgb(0, 0, 0);">procedure</span><span
 style="color: rgb(0, 0, 0);"> append(</span><span
 style="font-weight: bold; color: rgb(0, 0, 0);">var</span><span
 style="color: rgb(0, 0, 0);"> f: </span><span
 style="color: rgb(128, 0, 0);">text</span><span
 style="color: rgb(0, 0, 0);">); </span><span
 style="font-weight: bold; color: rgb(0, 0, 128);">inline</span><span
 style="color: rgb(0, 0, 0);">;</span>
</pre>
 Pristine freepascal append()<br>
 <a href="ex/append2.pas.html">Example</a><br>
 <a href="ex/append1.pas.html">See also the other avatar</a>
))}
procedure append(var f: text); inline;
begin
  system.append(f)
end;

//==============================================================================
//
//                     formatting and output functions
//
//==============================================================================

{**
  @abstract(@html(
<pre><span style="font-weight: bold; color: rgb(0, 0, 0);">procedure</span><span
 style="color: rgb(0, 0, 0);"> nam(</span><span
 style="font-weight: bold; color: rgb(0, 0, 0);">var</span><span
 style="color: rgb(0, 0, 0);"> p:mtx; name:</span><span
 style="color: rgb(128, 0, 0);">string</span><span
 style="color: rgb(0, 0, 0);">) </span><span
 style="font-weight: bold; color: rgb(0, 0, 128);">inline</span><span
 style="color: rgb(0, 0, 0);">;</span>
</pre>
  assign a <tt>name</tt> to <tt>p</tt><br>
  <a href="ex/nam_a.pas.html">Example</a>
))}
procedure nam(var a:mtx; name:string) inline;
begin
  a.nam := name
end;

{**
@abstract(@html(
<pre><span style="font-weight: bold; color: rgb(0, 0, 0);">function</span><span
 style="color: rgb(0, 0, 0);"> nam(p:mtx):</span><span
 style="color: rgb(128, 0, 0);">string</span><span
 style="color: rgb(0, 0, 0);">; </span><span
 style="font-weight: bold; color: rgb(0, 0, 128);">inline</span><span
 style="color: rgb(0, 0, 0);">;</span>
</pre>
  return the name of "p"<br>
  <a href="ex/nam_a.pas.html">Example</a>
))}
function nam(a:mtx):string; inline;
begin
  nam := a.nam;
end;

{**
@abstract(@html(
<pre><span style="font-weight: bold; color: rgb(0, 0, 0);">function</span><span
 style="color: rgb(0, 0, 0);"> mwr(a:mtx; </span><span
 style="font-weight: bold; color: rgb(0, 0, 0);">var</span><span
 style="color: rgb(0, 0, 0);"> f:</span><span
 style="color: rgb(128, 0, 0);">text</span><span
 style="color: rgb(0, 0, 0);">):</span><span
 style="color: rgb(128, 0, 0);">boolean</span><span
 style="color: rgb(0, 0, 0);">;</span>
</pre>
  Write matrix "a" in text stream "f"<br>
  The text stream "f" MUST be opened by append or rewrite
  or apn() or rwr() short name equivalents.<br>
  <a href="ex/mwr_af.pas.html">Example</a>
))}
function mwr(a:mtx; var f:text):boolean;
var
   i, j, ia, ia1: longint;
begin
  if a.nam<>'' then writeln(f,a.nam);
  ia1 := 1;
  for i := 1 to a.lgn do begin
    ia := ia1;
    if pretty then writx(f,a.e[ia]) else write(f,a.e[ia]);
    ia := ia + a.dl;
    for j := 2 to a.col do begin
      write(f,#9); // TODO use constant TAB
      // NOTE ascii char set assumed elsewhere
      // using exotic char set may require a lot of modifications
      if pretty then writx(f,a.e[ia]) else write(f,a.e[ia]);
      ia := ia + a.dl;
    end;
    ia1 := ia1 + a.dc;
    writeln(f);
  end;
  writeln(f);
  flush(f); // needed to ensure proper order of things in output
  mwr := true
end;

{**
@abstract(@html(
<pre><span style="font-weight: bold; color: rgb(0, 0, 0);">function</span><span
 style="color: rgb(0, 0, 0);"> mwr(a:mtx; s:</span><span
 style="color: rgb(128, 0, 0);">string</span><span
 style="color: rgb(0, 0, 0);">): </span><span
 style="color: rgb(128, 0, 0);">boolean</span><span
 style="color: rgb(0, 0, 0);">;</span>
</pre>
   This function output (append) matrix "a" to the text
   file named "s".<br> If "s" is an empty string
   then <em>stdout</em> is used.<br>
   rwr(s) to rewrite from an empty file.<br>
   or change append mode to false with rwr(true) : overwrite or
   rwr(false) ( default behaviour of mwr(a,s) ) to append
   matrix "a" to file "s"<br>
<a href="ex/mwr1.pas.html">Example</a>
))}
function mwr(a:mtx; s:string): boolean;
var
  f: text;
begin
  mwr := false;
  if apdmode then append(f,s) else rewrite(f,s);
  mwr(a,f);
  flush(f);
  close(f);
  mwr := true
end;

{**
@abstract(@html(
<pre><span style="font-weight: bold; color: rgb(0, 0, 0);">function</span><span
 style="color: rgb(0, 0, 0);"> mwr(a:mtx): </span><span
 style="color: rgb(128, 0, 0);">boolean</span><span
 style="color: rgb(0, 0, 0);">; </span><span
 style="font-weight: bold; color: rgb(0, 0, 128);">inline</span><span
 style="color: rgb(0, 0, 0);">;</span>
</pre>
  Write mtx "a" to stdout. <br>
<a href="ex/mwr0.pas.html">Example</a>
))}
function mwr(a:mtx): boolean; inline;
begin
  mwr := mwr(a,stdout);
end;

{**
@abstract(@html(
<pre><span style="font-weight: bold; color: rgb(0, 0, 0);">function</span><span
 style="color: rgb(0, 0, 0);"> mwr(x:</span><span
 style="color: rgb(128, 0, 0);">double</span><span
 style="color: rgb(0, 0, 0);">): </span><span
 style="color: rgb(128, 0, 0);">boolean</span><span
 style="color: rgb(0, 0, 0);">; </span><span
 style="font-weight: bold; color: rgb(0, 0, 128);">inline</span><span
 style="color: rgb(0, 0, 0);">;</span>
</pre>
Write a formatted number to stdout (as a 1,1 mtx) with a trailing blank line.<br>
<a href="ex/mwr_x.pas.html">Example</a>
))}
function mwr(x:double): boolean; inline;
begin
  dim(x11,1,1,x);
  mwr := mwr(x11,stdout)
end;

{**
@abstract(@html(
<pre><span style="font-weight: bold; color: rgb(0, 0, 0);">function</span><span
 style="color: rgb(0, 0, 0);"> mwr(x:</span><span
 style="color: rgb(128, 0, 0);">double</span><span
 style="color: rgb(0, 0, 0);">; </span><span
 style="font-weight: bold; color: rgb(0, 0, 0);">var</span><span
 style="color: rgb(0, 0, 0);"> f:</span><span
 style="color: rgb(128, 0, 0);">text</span><span
 style="color: rgb(0, 0, 0);">): </span><span
 style="color: rgb(128, 0, 0);">boolean</span><span
 style="color: rgb(0, 0, 0);">; </span><span
 style="font-weight: bold; color: rgb(0, 0, 128);">inline</span><span
 style="color: rgb(0, 0, 0);">;</span>
</pre>
Write a formatted number to an opened text stream (as a 1,1 mtx) with
a trailing blank line.<br>
<a href="ex/mwr_xf.pas.html">Example</a>
))}
function mwr(x:double; var f:text): boolean; inline;
begin
  dim(x11,1,1,x);
  mwr := mwr(x11,f)
end;

{**
@abstract(@html(
<pre><span style="font-weight: bold; color: rgb(0, 0, 0);">procedure</span><span
 style="color: rgb(0, 0, 0);"> mwr(s:</span><span
 style="color: rgb(128, 0, 0);">string</span><span
 style="color: rgb(0, 0, 0);">; </span><span
 style="font-weight: bold; color: rgb(0, 0, 0);">var</span><span
 style="color: rgb(0, 0, 0);"> f:</span><span
 style="color: rgb(128, 0, 0);">text</span><span
 style="color: rgb(0, 0, 0);">); </span><span
 style="font-weight: bold; color: rgb(0, 0, 128);">inline</span><span
 style="color: rgb(0, 0, 0);">;</span>
</pre>
Print string "s" into the open text stream "f".<br>
<a href="ex/mwrsf.pas.html">Example</a>
))}
procedure mwr(s:string; var f:text); inline;
begin
  writeln(f,s);
  flush(f)
end;

{**
@abstract(@html(
<pre><span style="font-weight: bold; color: rgb(0, 0, 0);">procedure</span><span
 style="color: rgb(0, 0, 0);"> mwr(s:</span><span
 style="color: rgb(128, 0, 0);">string</span><span
 style="color: rgb(0, 0, 0);">); </span><span
 style="font-weight: bold; color: rgb(0, 0, 128);">inline</span><span
 style="color: rgb(0, 0, 0);">;</span>
</pre>
Print the string "s" into stdout.<br>
<a href="ex/mwrs.pas.html">Example</a>
))}
procedure mwr(s:string); inline;
begin
  writeln(output,s);
  flush(output)
end;

{**
@abstract(@html(
<pre><span style="font-weight: bold; color: rgb(0, 0, 0);">procedure</span><span
 style="color: rgb(0, 0, 0);"> mwr(p:</span><span
 style="color: rgb(128, 0, 0);">boolean</span><span
 style="color: rgb(0, 0, 0);">, </span><span
 style="font-weight: bold; color: rgb(0, 0, 0);">var</span><span
 style="color: rgb(0, 0, 0);"> f:</span><span
 style="color: rgb(128, 0, 0);">text</span><span
 style="color: rgb(0, 0, 0);">); </span><span
 style="font-weight: bold; color: rgb(0, 0, 128);">inline</span><span
 style="color: rgb(0, 0, 0);">;</span>
</pre>
Print the boolean "p" into the opened text stream "f" as "TRUE" or "FALSE".<br>
<a href="ex/mwr_pf.pas.html">Example</a>
))}
procedure mwr(p:boolean; var f:text); inline;
begin
  if p then
    writeln(f,'TRUE')
  else
    writeln(f,'FALSE');
  flush(f)
end;

{**
@abstract(@html(
<pre><span style="font-weight: bold; color: rgb(0, 0, 0);">procedure</span><span
 style="color: rgb(0, 0, 0);"> mwr(p:</span><span
 style="color: rgb(128, 0, 0);">boolean</span><span
 style="color: rgb(0, 0, 0);">); </span><span
 style="font-weight: bold; color: rgb(0, 0, 128);">inline</span><span
 style="color: rgb(0, 0, 0);">;</span>
</pre>
Print the boolean "p" into stdout as 'TRUE' or 'FALSE'<br>
<a href="ex/mwr_p.pas.html">Example</a>
))}
procedure mwr(p:boolean); inline;
begin
  if p then
    writeln(output,'TRUE')
  else
    writeln(output,'FALSE');
  flush(output)
end;

{**
@abstract(@html(
<pre><span style="font-weight: bold; color: rgb(0, 0, 0);">procedure</span><span
 style="color: rgb(0, 0, 0);"> rwr(nam:</span><span
 style="color: rgb(128, 0, 0);">string</span><span
 style="color: rgb(0, 0, 0);">); </span><span
 style="font-weight: bold; color: rgb(0, 0, 128);">inline</span><span
 style="color: rgb(0, 0, 0);">;</span>
</pre>
This procedure will empty the file "nam"
or create it if doesn't exist.<br>
<a href="ex/rwr0.pas.html">Example</a>
))}
procedure rwr(nam:string); inline;
var f:text;
begin
  assign(f,nam);
  system.rewrite(f);
  close(f)
end;

{**
@abstract(@html(
<pre><span style="font-weight: bold; color: rgb(0, 0, 0);">procedure</span><span
 style="color: rgb(0, 0, 0);"> rwr(p:</span><span
 style="color: rgb(128, 0, 0);">boolean</span><span
 style="color: rgb(0, 0, 0);">); </span><span
 style="font-weight: bold; color: rgb(0, 0, 128);">inline</span><span
 style="color: rgb(0, 0, 0);">;</span>
</pre>
Modify append mode of mwr(a:mtx; s:string) procedure. Default is append,
rwr(true) will modify this behaviour to write a single mtx in a blank file
rather than appending to whatever pre-existing file. This flag will prevail
for all mwr(mtx,string) avatar. An alternative is to clear a file
with rwr(file_name) preserving the default append mode.<br>
<a href="ex/fmm.pas.html">Example</a>
))}
procedure rwr(p:boolean); inline;
begin
  apdmode := not p
end;

{**
@abstract(@html(
<pre><span style="font-weight: bold; color: rgb(0, 0, 0);">procedure</span><span
 style="color: rgb(0, 0, 0);"> rwr(nam:</span><span
 style="color: rgb(128, 0, 0);">string</span><span
 style="color: rgb(0, 0, 0);">; </span><span
 style="font-weight: bold; color: rgb(0, 0, 0);">var</span><span
 style="color: rgb(0, 0, 0);"> f:</span><span
 style="color: rgb(128, 0, 0);">text</span><span
 style="color: rgb(0, 0, 0);">); </span><span
 style="font-weight: bold; color: rgb(0, 0, 128);">inline</span><span
 style="color: rgb(0, 0, 0);">;</span>
</pre>
This procedure will empty the file "nam"
or create it if doesn't exist.<br>
<a href="ex/rwr1.pas.html">Example</a>
))}
procedure rwr(nam:string; var f:text); inline;
begin
  assign(f,nam);
  system.rewrite(f)
end;

{**
  @abstract(@html(
<pre><span style="font-weight: bold; color: rgb(0, 0, 0);">procedure</span><span
 style="color: rgb(0, 0, 0);"> tra(p: </span><span
 style="color: rgb(128, 0, 0);">boolean</span><span
 style="color: rgb(0, 0, 0);">); </span><span
 style="font-weight: bold; color: rgb(0, 0, 128);">inline</span><span
 style="color: rgb(0, 0, 0);">;</span>
</pre>
    Control trace mode; p=true/false : trace on/off
    When trace mode is on most functions will produce an information
    in standard output about the function name the dimensions of
    the const. mtx arguments and output the mtx result of the operation.<br>
    <a href="ex/tra_0.pas.html">Example</a>
))}
procedure tra(p: boolean); inline;
begin
  ptrace := p
end;

{**
  @abstract(@html(
<pre><span style="font-weight: bold; color: rgb(0, 0, 0);">procedure</span><span
 style="color: rgb(0, 0, 0);"> tra; </span><span
 style="font-weight: bold; color: rgb(0, 0, 128);">inline</span><span
 style="color: rgb(0, 0, 0);">;</span>
</pre>
    toggle trace mode on/off.
    When trace mode is on most functions will produce an information
    in standard output about the function name and the dimensions of
    the mtx arguments.<br>
    <a href="ex/tra_0.pas.html">Example</a>
))}
procedure tra; inline;
begin
  ptrace := not ptrace
end;

{**
@abstract(@html(
<pre><span style="font-weight: bold; color: rgb(0, 0, 0);">function</span><span
  style="color: rgb(0, 0, 0);"> ptr(a:mtx):Pdouble;</span>
</pre>
  This function enables to use plain mtx in an external fortran procedure call.<br>
  <a href="ex/ptr0.pas.html">Example</a><br>
  <a href="ex/sub1.f.html">The fortran code used</a>
))}
function ptr(a:mtx):Pdouble;
begin
  assert(a.e<>NIL);
  ptr := @a.e[1]
end;

{**
@abstract(@html(
<pre><span style="font-weight: bold; color: rgb(0, 0, 0);">procedure</span><span
 style="color: rgb(0, 0, 0);"> mym(</span><span
 style="font-weight: bold; color: rgb(0, 0, 0);">var</span><span
 style="color: rgb(0, 0, 0);"> r:mtx; </span><span
 style="font-weight: bold; color: rgb(0, 0, 0);">var</span><span
 style="color: rgb(0, 0, 0);"> x:</span><span
 style="color: rgb(128, 0, 0);">double</span><span
 style="color: rgb(0, 0, 0);">; lgn, col:</span><span
 style="color: rgb(128, 0, 0);">longint</span><span
 style="color: rgb(0, 0, 0);">);</span>
</pre>
    Memory map: use static memory in mtx objects.<br>
Matrix "r" is created using the memory starting at the location of the variable
"x" with the given numbers or rows and columns (in arg. lgn and col).
The number of elements of "r" must be equal or smaller than the static array
used. The second argument "x" is a variable argument of type double typically
the first element of the array: <tt>x[1]</tt>.
This is usefull to alias a static pascal array of double with an mtx object.
Some name convention may be used to highlight the connection between the mtx
and the static array used.<br>
    <a href="ex/mym0.pas.html">Example</a>
))}
procedure mym(var r:mtx; var x:double; lgn, col:longint);
var
  p: ^double;
  m1: longint = -1; // kludge to avoid compiler range warning below
begin
  if r.e=NIL then begin
    // allocate GSL compatibility records if needed
    new(r.mat);
    new(r.mat^.block);
    new(r.vec);
    new(r.vec^.block);
  end;
  p := @x;
  r.lgn := lgn;
  r.col := col;
  r.dl := 1;
  r.dc := col;
  r.dyn := false;
  r.max := lgn*col;
  r.smx := false; //
  r.e := @p[m1]; // this raises a range check warning with a constant
  with r.mat^ do begin
    size1 := lgn;
    size2 := col;
    tda := col;
    data := @r.e[1];
    owner := false;
    block^.size := lgn*col;
    block^.data := data
  end;
  with r.vec^ do begin
    size := lgn*col;
    stride := 1;
    data := @r.e[1];
    owner := false;
    block^.size := lgn*col;
    block^.data := data
  end
end;

{**
@abstract(@html(
<pre><span style="font-weight: bold; color: rgb(0, 0, 0);">function</span><span
 style="color: rgb(0, 0, 0);"> mty(a:mtx):</span><span
 style="color: rgb(128, 0, 0);">boolean</span><span
 style="color: rgb(0, 0, 0);">; </span><span
 style="font-weight: bold; color: rgb(0, 0, 128);">inline</span><span
 style="color: rgb(0, 0, 0);">;</span>
</pre>
  Return true if "a" is empty or not defined/allocated.<br>
  <a href="ex/mty0.pas.html">Example</a>
))}
function mty(a:mtx):boolean; inline;
begin
  mty := true;
  if a.e=nil then exit;
  mty := (a.col=0) or (a.lgn=0)
end;

//=============================================================================
//             assign and retrieve elements in a matrix
//=============================================================================

// here the most simple lgn col functions should be macro or inline and
// are needed by elm... should really be macro / or inlined...

{**
@abstract(@html(
<pre><span style="font-weight: bold; color: rgb(0, 0, 0);">function</span><span
 style="color: rgb(0, 0, 0);"> col(a:mtx): </span><span
 style="color: rgb(128, 0, 0);">longint</span><span
 style="color: rgb(0, 0, 0);">; </span><span
 style="font-weight: bold; color: rgb(0, 0, 128);">inline</span><span
 style="color: rgb(0, 0, 0);">;</span>
</pre>
  Return the number of columns of "a".<br>
  <a href="ex/axs0.pas.html">Example</a>
))}
function col(a:mtx): longint; inline;
begin
  col := a.col;
end;

{**
@abstract(@html(
<pre><span style="font-weight: bold; color: rgb(0, 0, 0);">function</span><span
 style="color: rgb(0, 0, 0);"> lgn(a:mtx): </span><span
 style="color: rgb(128, 0, 0);">longint</span><span
 style="color: rgb(0, 0, 0);">; </span><span
 style="font-weight: bold; color: rgb(0, 0, 128);">inline</span><span
 style="color: rgb(0, 0, 0);">;</span>
</pre>
Return the number of rows (french: lignes) of a.<br>
  <a href="ex/axs0.pas.html">Example</a>
))}
function lgn(a:mtx): longint; inline;
begin
  lgn := a.lgn;
end;

{**
@abstract(@html(
<pre><span style="font-weight: bold; color: rgb(0, 0, 0);">function</span><span
 style="color: rgb(0, 0, 0);"> elm(a:mtx; i, j: </span><span
 style="color: rgb(128, 0, 0);">longint</span><span
 style="color: rgb(0, 0, 0);">):</span><span
 style="color: rgb(128, 0, 0);">double</span><span
 style="color: rgb(0, 0, 0);">;</span>
</pre>
 Return a(i,j)<br>
 <a href="ex/axs0.pas.html">Example</a>
))}
function elm(a:mtx; i, j: longint):double;
var k: longint;
begin
  assert( (i>0) and (i<=a.lgn) );
  assert( (j>0) and (j<=a.col) );
  k := 1 + pred(i)*a.dc + pred(j)*a.dl;
  elm := a.e[k]
end;

{**
@abstract(@html(
<pre><span style="font-weight: bold; color: rgb(0, 0, 0);">procedure</span><span
 style="color: rgb(0, 0, 0);"> elm(a:mtx; i, j: </span><span
 style="color: rgb(128, 0, 0);">longint</span><span
 style="color: rgb(0, 0, 0);">; x:</span><span
 style="color: rgb(128, 0, 0);">double</span><span
 style="color: rgb(0, 0, 0);">);</span>
</pre>
Assign an element of a:  a(i,j) = x<br>
 <a href="ex/axs0.pas.html">Example</a>
))}
procedure elm(a:mtx; i, j: longint; x:double);
var k: longint;
begin
  assert( (i>0) and (i<=a.lgn) );
  assert( (j>0) and (j<=a.col) );
  k := 1 + pred(i)*a.dc + pred(j)*a.dl;
  a.e[k] := x
end;

{**
@abstract(@html(
<pre><span style="font-weight: bold; color: rgb(0, 0, 0);">function</span><span
 style="color: rgb(0, 0, 0);"> elm(a:mtx; b:mtx):</span><span
 style="color: rgb(128, 0, 0);">boolean</span><span
 style="color: rgb(0, 0, 0);">;</span>
</pre>
  <tt>a(i,j) = b(i,j)</tt> for all <tt>i,j</tt><br>
  <bf>b</bf> has soft dimensions i.e. rows and/or columns will be duplicated so that
  it matches the dimensions of <bf>a</bf>, if it make sense:
  i.e. b may be (1,1) (1,n) or (m,1) if a is (m,n).<br>
  <a href="ex/elm1.pas.html">Example 1</a><br>
  <a href="ex/elm2.pas.html">Example 2</a>
))}
function elm(a:mtx; b:mtx):boolean;
// TODO TODO TODO TODO TODO TODO TODO TODO TODO TODO TODO TODO TODO
// TODO check this again... should xpd() be used for softdim...
// probably YES TODO TODO TODO TODO TODO TODO TODO TODO TODO TODO
// TODO SHOULD IT BE TRACED OR NOT ????? TODO
var
  i, j, lgn1, col1, ia, ib, ia1, ib1: longint;
  bdl, bdc: longint;
begin
  elm := false;
  // "soft" dim of arg. b to match a either row, col or even scalar
  if lgn(b) = 1 then begin
    assert( (col(b)=col(a)) or (col(b)=1) );
    bdc := 0;
    bdl := b.dl;
    if col(b) = 1 then begin
      bdl := 0
    end
  end else if col(b) = 1 then begin
    assert( lgn(b)=lgn(a) );
    bdl := 0;
    bdc := b.dc;
  end else begin
    assert( lgn(a)=lgn(b) );
    assert( col(a)=col(b) );
    bdc := b.dc;
    bdl := b.dl
  end;
  lgn1 := a.lgn;
  col1 := a.col;
  ia1 := 1;
  ib1 := 1;
  for i := 1 to lgn1 do begin
    ia := ia1;
    ib := ib1;
    for j := 1 to col1 do begin
       a.e[ia] := b.e[ib];
       ia := ia + a.dl;
       ib := ib + bdl;
    end;
    ia1 := ia1 + a.dc;
    ib1 := ib1 + bdc;
  end;
  elm := true
end;

{**
  @abstract(@html(
<pre><span style="font-weight: bold; color: rgb(0, 0, 0);">procedure</span><span
 style="color: rgb(0, 0, 0);"> elm(b:mtx; x:</span><span
 style="color: rgb(128, 0, 0);">double</span><span
 style="color: rgb(0, 0, 0);">);</span>
</pre>
  Fill "b" with value "x";<br>
  "b" must have dimensions. Intended usage:
  to fill a slice with some special value.
  e.g. <tt> elm( col(a,1), 0.0) </tt><br>
  <a href="ex/elm2.pas.html">Example</a>
))}
procedure elm(b:mtx; x:double);
// TODO trace ? TODO
var i, j, lgn, col, ib, ib1: longint;
begin
  lgn := b.lgn;
  col := b.col;
  ib1 := 1;
  for i := 1 to lgn do begin
    ib := ib1;
    for j := 1 to col do begin
       b.e[ib] := x;
       ib := ib + b.dl
    end;
    ib1 := ib1 + b.dc
  end
end;

{**
@abstract(@html(
<pre><span style="font-weight: bold; color: rgb(0, 0, 0);">procedure</span><span
 style="color: rgb(0, 0, 0);"> elm(a:mtx; s:</span><span
 style="color: rgb(128, 0, 0);">string</span><span
 style="color: rgb(0, 0, 0);">);</span>
</pre>
assign elments of "a" with the values given in "s";
"a" can be a slice. The dimensions of "s" and "a" must match
in the softdim sense.<br>
<a href="ex/elm3.pas.html">Example</a>
))}
procedure elm(a:mtx; s:string);
var t: boolean;
begin
  t := ptrace;
  ptrace := false;
  equ(elm0,s);
  elm(a,elm0);
  ptrace := t;
  wtrace(a,'elm(mtx,string)')
end;

{**
@abstract(@html(
<pre><span style="font-weight: bold; color: rgb(0, 0, 0);">function</span><span
 style="color: rgb(0, 0, 0);"> mrd(</span><span
 style="font-weight: bold; color: rgb(0, 0, 0);">var</span><span
 style="color: rgb(0, 0, 0);"> r: mtx; </span><span
 style="font-weight: bold; color: rgb(0, 0, 0);">var</span><span
 style="color: rgb(0, 0, 0);"> f: </span><span
 style="color: rgb(128, 0, 0);">text</span><span
 style="color: rgb(0, 0, 0);">; mnam:</span><span
 style="color: rgb(128, 0, 0);">string</span><span
 style="color: rgb(0, 0, 0);">): </span><span
 style="color: rgb(128, 0, 0);">boolean</span><span
 style="color: rgb(0, 0, 0);">;</span>
</pre>
read in r a matrix named mnam from the opened text stream f.<br>
return true if read was successfull.<br>
<a href="ex/mrd_rfm.pas.html">Example</a><br>
))}
function mrd(var r: mtx; var f: text; mnam:string): boolean;
// NOTE : it MUST be a boolean function to allow while mrd(m)...
var
  s, s0, nam: string;
  i, j, k, m, n: longint;
  x: double;
  code: integer;
begin
  mrd := false;
  nam := '';
  if not read1st(mnam,f,s,nam,n) then exit;
  assert( n>0 );
  // first default alloc is n ( 1 line then later 2n 4n 8n... as needed )
  assert( mem(r,n) ); // sanity checks made here including "not r.smx"
  assert( not (r.e=NIL) ); // hopefully
  r.nam := nam;
  k := 0;
  m := 0;
  j := n;
  while true do begin
    if pos('#',s)=0 then begin
        inc(m);
        //if  m > ( mx div n ) then exit;
        for i := 1 to n do begin
            s0 := ExtractWord(i,s,WordDelim);
            assert(s0<>'');
            //writeln(stderr,'---------------------> s0= "',s0,'"');
            val(s0,x,code);
            //writeln(stderr,'---------------------> x= ',x:15:5);
            //assert(code=0);
            if code<>0 then exit;
            inc(k);
            if k > r.max then begin
              //writeln(stderr,'--- mrd: k, r.max, r.e: ',k,', ',r.max,', ',longint(r.e));
              assert( mem(r, r.max*2) );
              //writeln(stderr,'After realloc...  r.e=',longint(r.e),'  r.max=',r.max);
            end;
            r.e[k] := x;
        end;
    end;
    {$I-}
    readln(f,s);
    {$I+}
    j := WordCount(s,WordDelim);
    if j=0 then break;
    if j<>n then exit;
  end;
  r.dl := 1;
  r.dc := n;
  r.col := n;
  r.lgn := m;
{ ifdef gsl} // ***GSL*** interface records mat^ and vec^...
  with r.mat^ do begin
    size1 := m;
    size2 := n;
    tda := n;
  end;
  with r.vec^ do begin
    // we assume a plain mtx and a flat vector view in GSL
    size := m*n;
    stride := 1;
  end;
{ endif} // gsl ============================================
  wtrace(r,'mrd');
  mrd := true
end;

{**
@abstract(@html(
<pre><span style="font-weight: bold; color: rgb(0, 0, 0);">function</span><span
 style="color: rgb(0, 0, 0);"> mrd(</span><span
 style="font-weight: bold; color: rgb(0, 0, 0);">var</span><span
 style="color: rgb(0, 0, 0);"> r: mtx; </span><span
 style="font-weight: bold; color: rgb(0, 0, 0);">var</span><span
 style="color: rgb(0, 0, 0);"> f: </span><span
 style="color: rgb(128, 0, 0);">text</span><span
 style="color: rgb(0, 0, 0);">): </span><span
 style="color: rgb(128, 0, 0);">boolean</span><span
 style="color: rgb(0, 0, 0);">; </span><span
 style="font-weight: bold; color: rgb(0, 0, 128);">inline</span><span
 style="color: rgb(0, 0, 0);">;</span>
</pre>
    Read any matrix from text stream f in r; return true if successful.<br>
<a href="ex/mrd_rf.pas.html">Example</a>
))}
function mrd(var r: mtx; var f: text): boolean; inline;
begin
  mrd := mrd(r,f,'*')
end;

{**
@abstract(@html(
<pre><span style="font-weight: bold; color: rgb(0, 0, 0);">function</span><span
 style="color: rgb(0, 0, 0);"> mrd(</span><span
 style="font-weight: bold; color: rgb(0, 0, 0);">var</span><span
 style="color: rgb(0, 0, 0);"> r:mtx):</span><span
 style="color: rgb(128, 0, 0);">boolean</span><span
 style="color: rgb(0, 0, 0);">;</span>
</pre>
    Read any matrix from stdin in r; return true if successful.<br>
<a href="ex/mrd_r.pas.html">Example</a><br>
<a href="ex/mrd_r0.pas.html">Example</a>
))}
function mrd(var r:mtx):boolean;
begin
  mrd := mrd(r,input)
end;

{**
@abstract(@html(
<pre><span style="font-weight: bold; color: rgb(0, 0, 0);">function</span><span
 style="color: rgb(0, 0, 0);"> mrd(</span><span
 style="font-weight: bold; color: rgb(0, 0, 0);">var</span><span
 style="color: rgb(0, 0, 0);"> r:mtx; sf:</span><span
 style="color: rgb(128, 0, 0);">string</span><span
 style="color: rgb(0, 0, 0);">; sm:</span><span
 style="color: rgb(128, 0, 0);">string</span><span
 style="color: rgb(0, 0, 0);">):</span><span
 style="color: rgb(128, 0, 0);">boolean</span><span
 style="color: rgb(0, 0, 0);">;</span>
</pre>
   read "r" named "sm" from file "sf": return true if successful.<br>
   sm='*' : any matrix<br>
   sm=''  : first anonymous matrix in file<br>
  <a href="ex/mrd_rss.pas.html">Example</a>
))}
function mrd(var r:mtx; sf:string; sm:string):boolean;
var f: text;
begin
  mrd := false;   // TODO consistancy with exceptions handling elsewhere
  opn(sf,f);
  mrd := mrd(r,f,sm);
  clo(f)          // TODO uses system function
end;

{**
@abstract(@html(
<pre><span style="font-weight: bold; color: rgb(0, 0, 0);">function</span><span
 style="color: rgb(0, 0, 0);"> mrd(</span><span
 style="font-weight: bold; color: rgb(0, 0, 0);">var</span><span
 style="color: rgb(0, 0, 0);"> r:mtx; sf:</span><span
 style="color: rgb(128, 0, 0);">string</span><span
 style="color: rgb(0, 0, 0);">):</span><span
 style="color: rgb(128, 0, 0);">boolean</span><span
 style="color: rgb(0, 0, 0);">; </span><span
 style="font-weight: bold; color: rgb(0, 0, 128);">inline</span><span
 style="color: rgb(0, 0, 0);">;</span>
</pre>
   read "r" (first named or anonymous matrix) in the file "sf";
   return true if successful.<br>
<a href="ex/mrd_rs.pas.html">Example</a>
))}
function mrd(var r:mtx; sf:string):boolean; inline;
begin
  mrd := mrd(r,sf,'*')
end;

//==============================================================================
//======================     slice functions                          ==========
//==============================================================================

{ NOTE TODO move in private section............
  Return a compact slice of matrix a
}
function smx(a:mtx; i1,j1:longint; lgn, col:longint):mtx;
// This is the main slicing function other functions use it
// almost the same as sli() but no trace
var k: longint;
begin
  smx.smx := true;
  smx.max := 0; // ->freemem do nothing ? although .smx should be used
  smx.dyn := false;
  smx.lgn := lgn;
  smx.col := col;
  // calc. offset for this alias
  k := 1 + a.dc*(i1-1) + a.dl*(j1-1);
  // writeln(stderr,'k, e[k]: ',k:5,a.e[k]:10);
  // smx.e  := @a.e[k] - sizeof(double); // 0 based vectors but only use e[1..max]
  smx.e  := @a.e[k-1]; // 0 based vectors but only use e[1..max] TODO CHECK ?
  smx.t  := a.t;
  smx.dl := a.dl;
  smx.dc := a.dc;
  smx.nam := ''
end;

{ TODO move in private section
  @return( true if a is a slice i.e. not a plain mtx )
}
function smx(a:mtx):boolean;
begin
  smx := a.smx
end;

{**
@abstract(@html(
<pre><span style="font-weight: bold; color: rgb(0, 0, 0);">function</span><span
 style="color: rgb(0, 0, 0);"> col(a:mtx; j1: </span><span
 style="color: rgb(128, 0, 0);">longint</span><span
 style="color: rgb(0, 0, 0);">): mtx;</span>
</pre>
  Return a slice of a: <tt>a(:,j)</tt><br>
  <a href="ex/col1.pas.html">Example</a>
))}
function col(a:mtx; j: longint): mtx;
// TODO check if it can be inline TODO
begin
  col := smx(a, 1, j, a.lgn, 1)
end;

{**
@abstract(@html(
<pre><span style="font-weight: bold; color: rgb(0, 0, 0);">function</span><span
 style="color: rgb(0, 0, 0);"> col(a:mtx; j1, j2: </span><span
 style="color: rgb(128, 0, 0);">longint</span><span
 style="color: rgb(0, 0, 0);">): mtx;</span>
</pre>
Return a slice of a: <tt>a(:,j1:j2)</tt><br>
  <a href="ex/col1.pas.html">Example</a>
))}
function col(a:mtx; j1, j2: longint): mtx;
// TODO check if can be inline TODO
begin
  col := smx(a, 1, j1, a.lgn, j2-j1+1)
end;

{**
@abstract(@html(
<pre><span style="font-weight: bold; color: rgb(0, 0, 0);">function</span><span
 style="color: rgb(0, 0, 0);"> lgn(a:mtx; i1: </span><span
 style="color: rgb(128, 0, 0);">longint</span><span
 style="color: rgb(0, 0, 0);">): mtx;</span>
</pre>
Returns row "i1" of "a" (a slice of "a")<br>
<a href="ex/lgn0.pas.html">Example</a>
))}
function lgn(a:mtx; i1: longint): mtx;
// TODO should be renamed row for consistency with english
// TODO should be inline ?
begin
  lgn := smx(a, i1, 1, 1, a.col)
end;

{**
@abstract(@html(
<pre><span style="font-weight: bold; color: rgb(0, 0, 0);">function</span><span
 style="color: rgb(0, 0, 0);"> lgn(a:mtx; i1,i2: </span><span
 style="color: rgb(128, 0, 0);">longint</span><span
 style="color: rgb(0, 0, 0);">): mtx;</span>
</pre>
Return a slice of "a", adjacent rows from i1 to i2.<br>
Scilab: <tt>a(i1:i2,:)</tt><br>
<a href="ex/lgn1.pas.html">Example</a>
))}
function lgn(a:mtx; i1,i2: longint): mtx;
begin
  lgn := smx(a, i1, 1, i2-i1+1, a.col)
end;

//===================== ROI avatars of lgn() ==================================

{**
@abstract(@html(
<pre><span style="font-weight: bold; color: rgb(0, 0, 0);">function</span><span
 style="color: rgb(0, 0, 0);"> lgn(a:mtx; low, high:</span><span
 style="color: rgb(128, 0, 0);">double</span><span
 style="color: rgb(0, 0, 0);">):mtx;</span>
</pre>
Return a row SLICE of a, defined by a low high bracket on the
values of the first colum assumed to be monotonously increasing.
a MUST be plain (not a slice).<br>
<a href="ex/lgn2.pas.html">Example</a>
))}
function lgn(a:mtx; low, high:double):mtx;
var
  ia, i, i1, i2: longint;
begin
  assert( low<=high );
  assert( not a.smx );
  // first col of a should be ordered by increasing values
  ia := 1;
  i1 := 1;
  i2 := 0;
  for i := 1 to a.lgn do begin
    if a.e[ia]<low then begin
        inc(i1);
    end;
    if a.e[ia]<high then begin
        inc(i2)
    end else break;
    ia := ia + a.dc
  end;
  lgn := lgn(a,i1,i2)
end;


{**
@abstract(@html(
<pre><span style="font-weight: bold; color: rgb(0, 0, 0);">procedure</span><span
 style="color: rgb(0, 0, 0);"> lgn(</span><span
 style="font-weight: bold; color: rgb(0, 0, 0);">var</span><span
 style="color: rgb(0, 0, 0);"> a:mtx; ri, rx:mtx);</span>
</pre>
  TODO check this and preconditions...
  Selection and in situ deletion and shift of rows to
  select regions of interest in a histogram stored in "a".
  The first column of a must increase monotonously.<br>
  "ri": include rows if a(i,1) in any ri bracket
  "rx": exclude rows if a(i,1) in any rx bracket
  a should be plain (Warning: no sanity check done !)<br>
  replaced by rix() based on roi() and rex() cleaner/faster hopefully safer !
  <a href="ex/lgn_rix.pas.html">Example</a>
))}
procedure lgn(var a:mtx; ri, rx:mtx); inline;
begin
  rix(a,ri,rx)
end;

{**
@abstract(@html(
<pre><span style="font-weight: bold; color: rgb(0, 0, 0);">procedure</span><span
 style="color: rgb(0, 0, 0);"> lgn(</span><span
 style="font-weight: bold; color: rgb(0, 0, 0);">var</span><span
 style="color: rgb(0, 0, 0);"> a:mtx; ri:mtx);</span>
</pre>
Select rows of interest depending on values in the 1st column of "a" compared
with min max brackets given in each row of "ri".
Note: "a" size may decrease by removing all lines not in the roi brackets.
Assert that the first col "a" is sorted and monotonously increasing.<br>
Obsolete use roi() instead
<a href="ex/roi0.pas.html">Example</a><br>
))}
procedure lgn(var a:mtx; ri:mtx); inline;
begin
  roi(a,ri)
end;

{**
@abstract(@html(
<pre><span style="font-weight: bold; color: rgb(0, 0, 0);">function</span><span
 style="color: rgb(0, 0, 0);"> sli(a:mtx; l1, l2, c1, c2: </span><span
 style="color: rgb(128, 0, 0);">longint</span><span
 style="color: rgb(0, 0, 0);">): mtx;</span></pre>
 Return a slice of a; same as Scilab a(l1:l2,c1:c2)<br>
 <a href="ex/sli_1.pas.html">Example</a>
))}
function sli(a:mtx; l1, l2, c1, c2: longint): mtx; inline;
// TODO ??? why NOT use directly smx ??? and inline ???
// I avoid the name blk misleading with the GSL terminology...
// but we want a block of contiguous rows and columns.
// WARNING NOTE TODO no range check on lgn or col... ? it
// may or may not spawn an obvious error ! TODO
// TODO range check somewhere ... in col or and lgn ?
begin
  sli := col( lgn(a,l1,l2), c1, c2 )
end;

{**
@abstract(@html(
<pre><span style="font-weight: bold; color: rgb(0, 0, 0);">function</span><span
 style="color: rgb(0, 0, 0);"> trp(a:mtx):mtx;</span>
</pre>
Return a slice: transposed access to "a"<br>
))}
function trp(a:mtx):mtx;
// TODO implement an in situ transpose using ***GSL***
// procedure trp(var r:mtx; a:mtx);
// prototype:  trs(var a:mtx); ONLY for plain mtx of course !
begin
  trp := a;
  trp.t := not a.t; // TODO redundant ?
  trp.lgn := a.col;
  trp.col := a.lgn;
  trp.dl :=  a.dc;
  trp.dc :=  a.dl;
  trp.smx := true
end;

// TODO check that smx field is ok for all slice functions
// a slice function returns an mtx record...
// another special slice function
// TODO returns a row slice should it be col
// TODO check sane behaviour in various places
// including trp(dia(a))...

{**
@abstract(@html(
<pre><span style="font-weight: bold; color: rgb(0, 0, 0);">function</span><span
 style="color: rgb(0, 0, 0);"> dia(a:mtx):mtx;</span>
</pre>
Return a slice: diagonal elements of "a"<br>
"a" Must be a plain matrix.<br>
<a href="ex/dia0.pas.html">Example</a>

))}
function dia(a:mtx):mtx;
begin
  assert( not a.smx );
  dia := a;
  // dia.col := a.col;
  if a.col>a.lgn then dia.col:=a.lgn else dia.col:=a.col;
  dia.lgn := 1;
  dia.dc := a.dc+a.dl;
  dia.dl := a.dc+a.dl;
  dia.smx := true
end;

{**
@abstract(@html(
<pre><span style="font-weight: bold; color: rgb(0, 0, 0);">procedure</span><span
 style="color: rgb(0, 0, 0);"> zer(</span><span
 style="font-weight: bold; color: rgb(0, 0, 0);">var</span><span
 style="color: rgb(0, 0, 0);"> r:mtx; lgn, col: </span><span
 style="color: rgb(128, 0, 0);">longint</span><span
 style="color: rgb(0, 0, 0);">); </span><span
 style="font-weight: bold; color: rgb(0, 0, 128);">inline</span><span
 style="color: rgb(0, 0, 0);">;</span>
</pre>
Create "r(lgn,col)" filled with 0;<br>
Scilab: <tt>r = zeros(lgn,col);</tt><br>
<a href="ex/zer0.pas.html">Example</a>
))}
procedure zer(var r:mtx; lgn, col: longint); inline;
begin
  dim(r,lgn,col,0.0);
  wtrace(r,'zer')
end;

{**
@abstract(@html(
<pre><span style="font-weight: bold; color: rgb(0, 0, 0);">procedure</span><span
 style="color: rgb(0, 0, 0);"> zer(</span><span
 style="font-weight: bold; color: rgb(0, 0, 0);">var</span><span
 style="color: rgb(0, 0, 0);"> r:mtx; a:mtx); </span><span
 style="font-weight: bold; color: rgb(0, 0, 128);">inline</span><span
 style="color: rgb(0, 0, 0);">;</span>
</pre>
Create "r" same dimensions as "a" filled with 0; <br>
Scilab: r=zeros(a)<br>
<a href="ex/zer1.pas.html">Example</a>
))}
procedure zer(var r:mtx; a:mtx); inline;
begin
  dim(r,a.lgn,a.col,0.0);
  wtrace(r,a,'zer')
end;

{**
@abstract(@html(
<pre><span style="font-weight: bold; color: rgb(0, 0, 0);">procedure</span><span
 style="color: rgb(0, 0, 0);"> one(</span><span
 style="font-weight: bold; color: rgb(0, 0, 0);">var</span><span
 style="color: rgb(0, 0, 0);"> r:mtx; lgn, col: </span><span
 style="color: rgb(128, 0, 0);">longint</span><span
 style="color: rgb(0, 0, 0);">); </span><span
 style="font-weight: bold; color: rgb(0, 0, 128);">inline</span><span
 style="color: rgb(0, 0, 0);">;</span>
</pre>
Create "r(lgn,col)" filled with 1.0;<br>
Scilab: <tt>r = ones(lgn,col);</tt><br>
<a href="ex/one0.pas.html">Example</a>
))}
procedure one(var r:mtx; lgn, col: longint); inline;
begin
  dim(r,lgn,col,1.0);
  wtrace(r,'one')
end;

{**
@abstract(@html(
<pre><span style="font-weight: bold; color: rgb(0, 0, 0);">procedure</span><span
 style="color: rgb(0, 0, 0);"> one(</span><span
 style="font-weight: bold; color: rgb(0, 0, 0);">var</span><span
 style="color: rgb(0, 0, 0);"> r:mtx; a:mtx); </span><span
 style="font-weight: bold; color: rgb(0, 0, 128);">inline</span><span
 style="color: rgb(0, 0, 0);">;</span>
</pre>
Create "r" of the same dim. as "a" filled with 1.0;<br>
Scilab: <tt>r=ones(a)</tt><br>
<a href="ex/one0.pas.html">Example</a>
))}
procedure one(var r:mtx; a:mtx); inline;
begin
  dim(r,a.lgn,a.col,1.0);
  wtrace(r,a,'one')
end;

{**
@abstract(@html(
<pre><span style="font-weight: bold; color: rgb(0, 0, 0);">procedure</span><span
 style="color: rgb(0, 0, 0);"> eye(</span><span
 style="font-weight: bold; color: rgb(0, 0, 0);">var</span><span
 style="color: rgb(0, 0, 0);"> r:mtx; lgn, col: </span><span
 style="color: rgb(128, 0, 0);">longint</span><span
 style="color: rgb(0, 0, 0);">);</span>
</pre>
Equivalent to Scilab <tt>r = eye(lgn,col)</tt><br>
<a href="ex/eye0.pas.html">Example</a>
))}
procedure eye(var r:mtx; lgn, col: longint);
var
  n, k: longint;
begin
  dim(r,lgn,col,0.0); // GSL compat in dim()
  k := 1;
  n := col;
  if n > lgn then n := lgn;
  while n>0 do begin
    r.e[k] := 1.0;
    dec(n);
    k := k + col + 1
  end;
  wtrace(r,'eye')
end;

{**
@abstract(@html(
<pre><span style="font-weight: bold; color: rgb(0, 0, 0);">procedure</span><span
 style="color: rgb(0, 0, 0);"> eye(</span><span
 style="font-weight: bold; color: rgb(0, 0, 0);">var</span><span
 style="color: rgb(0, 0, 0);"> r:mtx; a:mtx); </span><span
 style="font-weight: bold; color: rgb(0, 0, 128);">inline</span><span
 style="color: rgb(0, 0, 0);">;</span>
</pre>
Equivalent to Scilab <tt>r = eye(a)</tt><br>
<a href="ex/eye0.pas.html">Example</a>
))}
procedure eye(var r:mtx; a:mtx); inline;
begin
  eye(r,a.lgn,a.col)
end;

{**
@abstract(@html(
<pre><span style="font-weight: bold; color: rgb(0, 0, 0);">procedure</span><span
 style="color: rgb(0, 0, 0);"> lsp(</span><span
 style="font-weight: bold; color: rgb(0, 0, 0);">var</span><span
 style="color: rgb(0, 0, 0);"> r:mtx; x0, x1: </span><span
 style="color: rgb(128, 0, 0);">double</span><span
 style="color: rgb(0, 0, 0);">; n:</span><span
 style="color: rgb(128, 0, 0);">longint</span><span
 style="color: rgb(0, 0, 0);">);</span>
</pre>
Returns "r", a column vector of "n"+1 elements evenly spaced
from "x0" to "x1". Same as Scilab's linspace() function.<br>
<a href="ex/lsp0.pas.html">Example</a><br>

))}
procedure lsp(var r:mtx; x0, x1: double; n:longint);
// Scilab linspace -> r[1,n+1]
var
  k: longint;
begin
  dim(r,n+1,1,0.0);
  for k := 0 to n do begin
    r.e[k+1] := (x0*(n-k)+x1*k)/n;
  end;
  r.nam := '';
  wtrace(r,'lsp');
end;

{**
@abstract(@html(
<pre><span style="font-weight: bold; color: rgb(0, 0, 0);">procedure</span><span
 style="color: rgb(0, 0, 0);"> lsp(a:mtx; x0, x1:</span><span
 style="color: rgb(128, 0, 0);">double</span><span
 style="color: rgb(0, 0, 0);">);</span>
</pre>
Enable to fill a column or row
of a table with linearly spaced values.
Intended use: lsp(col(table,1),0.0,100.0);
to fill the first col. of a table with lin. spaced values of the
independant variable.<br>
<a href="ex/lsp1.pas.html">Example</a><br>
))}
procedure lsp(a:mtx; x0, x1:double);
var
  i,j,k,m,n,p: longint;
begin
  m := lgn(a);
  n := col(a);
  p := m*n-1;
  k := 0;
  for i := 1 to m do begin
    for j := 1 to n do begin
      elm(a,i,j,(x0*(p-k)+x1*k)/p);
      k := succ(k)
    end
  end
end;

{**
@abstract(@html(
<pre><span style="font-weight: bold; color: rgb(0, 0, 0);">procedure</span><span
 style="color: rgb(0, 0, 0);"> equ(</span><span
 style="font-weight: bold; color: rgb(0, 0, 0);">var</span><span
 style="color: rgb(0, 0, 0);"> r:mtx; x:</span><span
 style="color: rgb(128, 0, 0);">double</span><span
 style="color: rgb(0, 0, 0);">; m, n:</span><span
 style="color: rgb(128, 0, 0);">longint</span><span
 style="color: rgb(0, 0, 0);">); </span><span
 style="font-weight: bold; color: rgb(0, 0, 128);">inline</span><span
 style="color: rgb(0, 0, 0);">;</span>
</pre>
Create matrix "r" with "m" rows and "n" columns, with all
elements equal to "x".<br>
<a href="ex/equrxmn.pas.html">Example</a>
))}
procedure equ(var r:mtx; x:double; m, n:longint); inline;
begin
  dim(r,m,n,x);
  wtrace(r,'equ');
end;

{**
@abstract(@html(
<pre><span style="font-weight: bold; color: rgb(0, 0, 0);">procedure</span><span
 style="color: rgb(0, 0, 0);"> equ(</span><span
 style="font-weight: bold; color: rgb(0, 0, 0);">var</span><span
 style="color: rgb(0, 0, 0);"> r:mtx; x: </span><span
 style="color: rgb(128, 0, 0);">double</span><span
 style="color: rgb(0, 0, 0);">; a:mtx); </span><span
 style="font-weight: bold; color: rgb(0, 0, 128);">inline</span><span
 style="color: rgb(0, 0, 0);">;</span>
</pre>
Create matrix "r" of the same size as "a" all elements of "r" equal to "x".<br>
<a href="ex/equrxa.pas.html">Example</a>
))}
procedure equ(var r:mtx; x: double; a:mtx); inline;
begin
  dim(r,a.lgn,a.col,x);
  wtrace(r,a,'equ');
end;

{**
@abstract(@html(
<pre><span style="font-weight: bold; color: rgb(0, 0, 0);">procedure</span><span
 style="color: rgb(0, 0, 0);"> equ(</span><span
 style="font-weight: bold; color: rgb(0, 0, 0);">var</span><span
 style="color: rgb(0, 0, 0);"> r:mtx; a:mtx); </span><span
 style="font-weight: bold; color: rgb(0, 0, 128);">inline</span><span
 style="color: rgb(0, 0, 0);">;</span>
</pre>
 Create a copy "r" of matrix "a"; "r" will be a plain matrix whereas "a"
 can be a slice.<br>
 <a href="ex/equra.pas.html">Example</a>
 ))}
procedure equ(var r:mtx; a:mtx); inline;
begin
  cpy(r,a);
  wtrace(r,a,'equ')
end;

{**
@abstract(@html(
<pre><span style="font-weight: bold; color: rgb(0, 0, 0);">procedure</span><span
 style="color: rgb(0, 0, 0);"> equ(</span><span
 style="font-weight: bold; color: rgb(0, 0, 0);">var</span><span
 style="color: rgb(0, 0, 0);"> r:mtx; s:</span><span
 style="color: rgb(128, 0, 0);">string</span><span
 style="color: rgb(0, 0, 0);">);</span>
</pre>
Convert the string representation "s" of a matrix into "r";
Within a row elements are separated by spaces or commas and rows
are separated by semi-colons.<br>

Alternatively s can be filename[,mtxname] to read from. This is
redundant with mrd() but convenient in certain cases e.g. to
communicate a table of value smaal (big) in a command line option limited to
255 char. in free Pascal;

<a href="ex/equrs.pas.html">Example</a>
))}
procedure equ(var r:mtx; s:string);
const
  delim = [','];
var
  m: longint = 0;
  n: longint = 0;
  s0, s1: string;
  i, j: longint;
  x: double;
  code, nw: integer;
  fnam, mnam: string;
begin
  if pos('/',s)>0 then begin
    // it is a filename[,mtxname]
    mnam := '*';
    fnam := ExtractWord(1,s,delim);
    nw := WordCount(s,delim);
    if nw>1 then mnam := ExtractWord(2,s,delim);
    mrd(r,fnam,mnam);
    s0 := 'equ:'+s;
    wtrace(r,s0)
  end else begin
    m := WordCount(s,RowDelim);
    if m=0 then begin // empty matrix
       dim(r,0,0);
       wtrace(r,'equ');
       exit
    end;
    for i := 1 to m do begin
      s0 := ExtractWord(i,s,RowDelim);
      if n=0 then begin
        // mem. alloc. and dim.
        n := WordCount(s0,ElmDelim);
        assert( n>0 ); // TODO should this be handled in dim ? TODO
        dim(r,m,n,0.0);
      end else begin
        // sanity check: TODO may be relax and fill with 0 missing col. ?
        assert(WordCount(s0,ElmDelim)=n)
      end;
      for j := 1 to n do begin
        s1 := ExtractWord(j,s0,ElmDelim);
        val(s1,x,code);
        assert(code=0); // TODO ? or NaN or 0 ?
        elm(r,i,j,x)
      end
    end;
    wtrace(r,'equ')
  end;
end;

{**
@abstract(@html(
<pre><span style="font-weight: bold; color: rgb(0, 0, 0);">procedure</span><span
 style="color: rgb(0, 0, 0);"> cat(</span><span
 style="font-weight: bold; color: rgb(0, 0, 0);">var</span><span
 style="color: rgb(0, 0, 0);"> r:mtx; a,b:mtx);</span>
</pre>
Matrix right concatenation. <br>
Scilab: <tt>r = [a,b]</tt><br>
<a href="ex/cat1.pas.html">Example</a>
))}
procedure cat(var r:mtx; a,b:mtx);
var
  ca, cb, la: longint;
begin
  la := lgn(a);
  assert((lgn(b)=la));
  ca := col(a);
  cb := col(b);
  dim(r,la,ca+cb);
  elm(col(r,1,ca),a);
  elm(col(r,ca+1,ca+cb),b);
  wtrace(r,a,b,'cat');
end;

{**
@abstract(@html(
<pre><span style="font-weight: bold; color: rgb(0, 0, 0);">procedure</span><span
 style="color: rgb(0, 0, 0);"> cat(</span><span
 style="font-weight: bold; color: rgb(0, 0, 0);">var</span><span
 style="color: rgb(0, 0, 0);"> r:mtx; a:mtx);</span></pre>
 Matrix right concatenation; i.e. <tt>r = [r,a]</tt> using Scilab notation.<br>
<a href="ex/stkcat.pas.html">Example</a>
))}
procedure cat(var r:mtx; a:mtx);
begin
  if mty(a) then exit
  else if mty(r) then cpy(r,a)
  else begin
    cpy(cat0,r);
    cat(r,cat0,a)
  end
end;

{**
@abstract(@html(
<pre><span style="font-weight: bold; color: rgb(0, 0, 0);">procedure</span><span
 style="color: rgb(0, 0, 0);"> stk(</span><span
 style="font-weight: bold; color: rgb(0, 0, 0);">var</span><span
 style="color: rgb(0, 0, 0);"> r:mtx; a,b:mtx);</span>
</pre>
Stack "a" and "b" and copy in "r"; "a" and "b" must have the same number
of columns.<br>
Note: empty mtx are not handled properly, use the other avatar.<br>
Scilab: <tt>r=[a;b]</tt><br>
<a href="ex/stk1.pas.html">Example</a>
))}
procedure stk(var r:mtx; a,b:mtx);
// TODO check this: in this avatar empty mtx NOT handled properly
// TODO use the other version stk(r,a) !
var
  la, lb, ca: longint;
begin
  ca := col(a);
  assert(col(b)=ca);
  la := lgn(a);
  lb := lgn(b);
  dim(r,la+lb,ca);
  elm(lgn(r,1,la),a);
  elm(lgn(r,la+1,la+lb),b);
  wtrace(r,a,b,'stk');
end;

{**
@abstract(@html(
<pre><span style="font-weight: bold; color: rgb(0, 0, 0);">procedure</span><span
 style="color: rgb(0, 0, 0);"> stk(</span><span
 style="font-weight: bold; color: rgb(0, 0, 0);">var</span><span
 style="color: rgb(0, 0, 0);"> r:mtx; a:mtx);</span>
</pre>
Append "a" below "r"; "a" and "r" must have the same number of columns or
be empty.<br>
Scilab: <tt>r = [r;a]</tt><br>
<a href="ex/stk2.pas.html">Example</a>
))}
procedure stk(var r:mtx; a:mtx);
begin
  if mty(a) then exit
  else if mty(r) then cpy(r,a)
  else begin
    cpy(cat0,r);
    stk(r,cat0,a)
  end
end;

//=============================================================================
// ancillary output for dbug... depends on smx()
//=============================================================================

{**
@abstract(@html(
output fields of <tt>a</tt>
This procedure is intended for debug only.<br>
<a href="ex/mpt_.pas.html">Example</a>
))}
procedure mpt(a:mtx);
begin
  writeln('# ====================================');
  writeln('# smx, trp, dyn: ',a.smx,' ',a.t,' ',a.dyn);
  writeln('# max, lgn, col, dl, dc: ',a.max,'   ',a.lgn,' ',a.col,'   ',a.dl,' ',a.dc);
  writeln('# ====================================');
  writeln;
{ ifdef gsl}
  if (not smx(a)) and (a.e<>NIL) then begin // ***GSL*** data records...
    writeln('a.mat^.block^.size: ',a.mat^.block^.size);
    writeln('a.mat^.block^.data: ',longint(a.mat^.block^.data));
    writeln('a.mat^.size1:       ',a.mat^.size1);
    writeln('a.mat^.size2:       ',a.mat^.size2);
    writeln('a.mat^.tda:         ',a.mat^.tda);
    writeln('a.mat^.data:        ',longint(a.mat^.data));
    writeln('a.vec^.block^.size: ',a.vec^.block^.size);
    writeln('a.vec^.block^.data: ',longint(a.vec^.block^.data));
    writeln('a.vec^.size:        ',a.vec^.size);
    writeln('a.vec^.stride:      ',a.vec^.stride);
    writeln('a.vec^.data:        ',longint(a.vec^.data));
  end;
{ endif} // gsl ==========================================
  flush(output);
end;

{**
@abstract(@html(
output fields of <tt>a</tt> along with the comment string <tt>s</tt>
This procedure is intended for debug only.<br>
<a href="ex/mpt_.pas.html">Example</a>
))}
procedure mpt(a:mtx; s:string); inline;
begin
  writeln('# ========== mpt ====',s,'====');
  flush(output);
  mpt(a);
end;

//=============================================================================

{$include opr.inc}
{$include add.inc}
{$include sub.inc}
{$include mul.inc}
{$include dvd.inc}

//=============================================================================

{**
@abstract(@html(
<pre><span style="font-weight: bold; color: rgb(0, 0, 0);">procedure</span><span
 style="color: rgb(0, 0, 0);"> prd(</span><span
 style="font-weight: bold; color: rgb(0, 0, 0);">var</span><span
 style="color: rgb(0, 0, 0);"> r: mtx; a, b: mtx);</span>
</pre>
matrix product; Scilab: <tt>r = a * b</tt><br>
The dimensions of "a" and "b" must be compatible with matrix product.<br>
<a href="ex/prd1.pas.html">Example</a>
))}
procedure prd(var r: mtx; a, b: mtx);
var
  m, n, p: longint;
  i, j, k: longint;
  ia, ib, ir: longint;
  ia1, ib1, ir1: longint;
  t0: double;
begin
  chkprd(r,a,b);
  m := a.lgn;
  n := a.col;
  p := b.col;
  ia1 := 1;
  ib1 := 1;
  ir1 := 1;
  for i := 1 to m do begin
    ir := ir1;
    for j := 1 to p do begin
      t0 := 0;
      ia := ia1;
      ib := ib1;
      for k := 1 to n do begin
        t0 := t0 + a.e[ia]*b.e[ib];
        ia := ia + a.dl;
        ib := ib + b.dc;
      end;
      //writeln(stderr,i:5,j:5,ir:5);
      r.e[ir] := t0;
      ir := ir + r.dl;
      ib1 := ib1 + b.dl;
    end;
    ia1 := ia1 + a.dc;
    ib1 := 1;
    ir1 := ir1 + r.dc;
  end;
  wtrace(r,a,b,'prd')
end;

{**
@abstract(@html(
<pre><span style="font-weight: bold; color: rgb(0, 0, 0);">function</span><span
 style="color: rgb(0, 0, 0);"> prd(a, b:mtx):</span><span
 style="color: rgb(128, 0, 0);">double</span><span
 style="color: rgb(0, 0, 0);">; </span><span
 style="font-weight: bold; color: rgb(0, 0, 128);">inline</span><span
 style="color: rgb(0, 0, 0);">;</span>
</pre>
Returns the scalar product of vectors "a" and "b"; any combination
of row or column vectors with the same number of elements;<br>
<a href="ex/prd_xab.pas.html">Example</a>
))}
function prd(a, b:mtx):double; inline;
var
  a1, b1: mtx; // map a and b to row & col vectors
begin
  if lgn(a)>1 then a1 := trp(a) else a1 := a;
  if col(b)>1 then b1 := trp(b) else b1 := b;
  prd(x11,a1,b1);
  prd := elm(x11,1,1)
end;

{**
@abstract(@html(
<pre><span style="font-weight: bold; color: rgb(0, 0, 0);">function</span><span
 style="color: rgb(0, 0, 0);"> inv(</span><span
 style="font-weight: bold; color: rgb(0, 0, 0);">var</span><span
 style="color: rgb(0, 0, 0);"> r:mtx; a:mtx; </span><span
 style="font-weight: bold; color: rgb(0, 0, 0);">var</span><span
 style="color: rgb(0, 0, 0);"> d:</span><span
 style="color: rgb(128, 0, 0);">double</span><span
 style="color: rgb(0, 0, 0);">):</span><span
 style="color: rgb(128, 0, 0);">boolean</span><span
 style="color: rgb(0, 0, 0);">;</span>
</pre>
   Matrix inversion by Gauss/Jordan with maximum pivot. The function
returns <bf>true</bf> if successful.
<ul>
   <li> <bf>r</bf>: inverse of a </li>
   <li> <bf>a</bf>: plain mtx or a slice; not modified </li>
   <li> <bf>d</bf>: determinant of a </li>
</ul>
<a href="ex/inv_rad.pas.html">Example</a>
))}
function inv(var r:mtx; a:mtx; var d:double):boolean; inline;
begin
  cpy(r,a); // private function no trace as in equ()
  inv := matinv(r,d);
end;

{**
@abstract(@html(
<pre><span style="font-weight: bold; color: rgb(0, 0, 0);">function</span><span
 style="color: rgb(0, 0, 0);"> inv(a:mtx; </span><span
 style="font-weight: bold; color: rgb(0, 0, 0);">var</span><span
 style="color: rgb(0, 0, 0);"> d:</span><span
 style="color: rgb(128, 0, 0);">double</span><span
 style="color: rgb(0, 0, 0);">):</span><span
 style="color: rgb(128, 0, 0);">boolean</span><span
 style="color: rgb(0, 0, 0);">; </span><span
 style="font-weight: bold; color: rgb(0, 0, 128);">inline</span><span
 style="color: rgb(0, 0, 0);">;</span>
</pre>
 Matrix inversion by Gauss/Jordan with maximum pivot. The function
 returns <bf>true</bf> if successful.
<ul>
   <li> <bf>a</bf>: plain mtx, inverted in situ.</li>
   <li> <bf>d</bf>: determinant of "a" </li>
</ul>
<a href="ex/inv_rd.pas.html">Example</a>
))}
function inv(a:mtx; var d:double):boolean; inline;
begin
  assert( not a.smx );
  inv := matinv(a,d);
end;

{**
@abstract(@html(
<pre><span style="font-weight: bold; color: rgb(0, 0, 0);">function</span><span
 style="color: rgb(0, 0, 0);"> spl( </span><span
 style="font-weight: bold; color: rgb(0, 0, 0);">var</span><span
 style="color: rgb(0, 0, 0);"> s:mtx; x, y: mtx):</span><span
 style="color: rgb(128, 0, 0);">boolean</span><span
 style="color: rgb(0, 0, 0);">;</span>
</pre>
  Wrapper of the <a href="ex/spline.f.html">fortran source</a>.<br>
<ul>
<li> x independant variable vector </li>
<li> y dependant variable vector </li>
<li> s the returnes spline table of coefficient</li>
</ul>
  See also sev()<br>
<a href="ex/fmm.pas.html">Example</a>
))}
function spl( var s:mtx; x, y: mtx):boolean;
var
  sx, sy, b, c, d: mtx;
  n: longint;
begin
  spl := false;
  assert( (x.lgn = y.lgn) and (x.col=y.col) );
  assert( (x.lgn=1) or (x.col=1) );
  assert( x.lgn<>x.col );
  if col(x)<lgn(x) then begin
    x := trp(x);
    y := trp(y);
  end;
  // sanity check of arguments... lgn or col vectors should be ok
  // so may be transpose if col vector..
  // sev() should also handle gracefully line/col vectors TODO
  // m := lgn(x);
  n := col(x);
  //assert( ( col(y)=n ) and (n>1)  );
  //assert( (lgn(y)=1) and (m=1) );
  // build s matrix with proper dim and copy x and y in first 2 lines
  dim(s, 5, n);
  elm(lgn(s, 1),x);
  elm(lgn(s, 2),y);
  // line 1 and 2 are copy of x and y and contiguous
  // so x and y arg. can be slices TODO check this
  // lines 3 4 5 will hold b, c and d coeff.
  sx := lgn(s,1);
  sy := lgn(s,2);
  b := lgn(s,3);
  c := lgn(s,4);
  d := lgn(s,5);
  spline_(n, ptr(sx), ptr(sy), ptr(b), ptr(c), ptr(d) );
  spl := true
end;

{**
@abstract(@html(
<pre><span style="font-weight: bold; color: rgb(0, 0, 0);">function</span><span
 style="color: rgb(0, 0, 0);"> sev(u:</span><span
 style="color: rgb(128, 0, 0);">double</span><span
 style="color: rgb(0, 0, 0);">; s:mtx):</span><span
 style="color: rgb(128, 0, 0);">double</span><span
 style="color: rgb(0, 0, 0);">;</span>
</pre>
<ul>
  <li>"s" spline interpolation table calculated by spl</li>
  <li>return: interpolated value for independant variable = "u"</li>
</ul>
<a href="ex/fmm.pas.html">Example</a>
))}
function sev(u:double; s:mtx):double;
var
  x, y, b, c, d: pointer;
  n: longint;
begin
  assert( (lgn(s)=5) and (col(s)>1 ) );
  x := ptr( lgn(s,1) );
  y := ptr( lgn(s,2) );
  b := ptr( lgn(s,3) );
  c := ptr( lgn(s,4) );
  d := ptr( lgn(s,5) );
  n := col(s);
  sev := seval_(n,u,x,y,b,c,d);
end;

{**
@abstract(@html(
<pre><span style="font-weight: bold; color: rgb(0, 0, 0);">function</span><span
 style="color: rgb(0, 0, 0);"> sev(</span><span
 style="font-weight: bold; color: rgb(0, 0, 0);">var</span><span
 style="color: rgb(0, 0, 0);"> r:mtx; u:mtx; s:mtx):</span><span
 style="color: rgb(128, 0, 0);">boolean</span><span
 style="color: rgb(0, 0, 0);">;</span>
</pre>
<ul>
  <li>"r": interpolated table for the vector u and spline table s</li>
  <li>"u": a vector sorted by increasing values of the independant variable</li>
           "u" must be plain or a compact slice (1 row but NOT 1 column).</li>
  <li>"s": a cubic spline interpolation table calculated by <tt>spl()</tt></li>
  <li>return: always true (useless);
              note todo exception handling harmonisation ?</li>
</ul>
  <a href="ex/fmm.pas.html">Example</a>
))}
function sev(var r:mtx; u:mtx; s:mtx):boolean;
var
  x, y, b, c, d: pointer;
  n: longint;
  k: longint;
  uk: double;
begin
  sev := false;
  assert( not u.smx ); // interim u is plain mtx
  assert( (lgn(s)=5) and (col(s)>1 ) );
  x := ptr( lgn(s,1) );
  y := ptr( lgn(s,2) );
  b := ptr( lgn(s,3) );
  c := ptr( lgn(s,4) );
  d := ptr( lgn(s,5) );
  n := col(s);
  dim(r,lgn(u),col(u),0.0);
  for k := 1 to col(u)*lgn(u) do begin
    uk := u.e[k];
    r.e[k] := seval_(n,uk,x,y,b,c,d);
  end;
  sev := true
end;

{**
@abstract(@html(
map a mtx into a GSL vector pointer.<br>
<ul>
<li>"a" a plain mtx (or compact slice TODO check this);</li>
<li>return a GSL's vector pointer.</li>
</ul>
<a href="http://www.gnu.org/software/gsl/">The Gnu Scientific Library</a>
))}
function vec(a:mtx):Pvector;
begin
  assert( not a.smx );
  vec := a.vec
end;

{**
@abstract(@html(
<pre><span style="font-weight: bold; color: rgb(0, 0, 0);">function</span><span
 style="color: rgb(0, 0, 0);"> mat(a:mtx):Pmatrix;</span>
</pre>
map mtx "a" into a GSL matrix pointer (returned by the function).<br>
Note: this is just a hook to use a very small part of the GSL library and
mostly to test the concept.
<a href="http://www.gnu.org/software/gsl/">The Gnu Scientific Library</a>
))}
function mat(a:mtx):Pmatrix;
begin
  assert( not a.smx );
  mat := a.mat
end;


//==============================================================================

// selection of contiguous blocks of rows and in situ shift
// of excluded zones... function roi ( regions of interest and
// rex() excluded regions ... usage for peak bkg processing...
// try to fix a bug in previous lgn() avatar

{**
@abstract(@html(
<pre><span style="font-weight: bold; color: rgb(0, 0, 0);">procedure</span><span
 style="color: rgb(0, 0, 0);"> rex(</span><span
 style="font-weight: bold; color: rgb(0, 0, 0);">var</span><span
 style="color: rgb(0, 0, 0);"> m:mtx; lo, hi: </span><span
 style="color: rgb(128, 0, 0);">double</span><span
 style="color: rgb(0, 0, 0);">);</span>
  </pre>
  <a href="ex/trex.pas.html">Example</a>
))}
procedure rex(var m:mtx; lo, hi: double);
var
  r, c: longint; // current # of rows and columns of m
  il: longint;   // index of first row in the excluded region [1..n+1]
  ih: longint;   // index of last row in excl. reg. [1..n+1]
  nx: longint;   // # of rows in the excluded region
  na: longint;   // # of rows above the excluded region
  no: longint;   // # of bytes in block above region
  i0: longint;   // index of first elem. of excluded region
  i1: longint;   // index of first elem. of block above region
begin
  assert( not m.smx );
  assert( not m.t );
  r := m.lgn;
  c := m.col;
  il := fnd(m,lo);   // TODO should fnd return 0 instead of n+1 ?
  if il>r then exit; // empty region do nothing TODO see above
  ih := pred(fnd(m,hi));
  nx := ih - il + 1;
  if nx<=0 then exit; // region empty
  na := r - nx - il + 1;
  if na > 0 then begin
    // shift a contiguous non empty block of elements
    // down over the excluded region...
    no := na * c * sizeof(double);
    i0 := idx(m,il,1);
    i1 := idx(m,ih+1,1);
    move(m.e[i1],m.e[i0],no);
  end;
  // ...update the number of rows of mtx along with gsl hooks
  dim(m,lgn(m)-nx,col(m))
end;

{**
@abstract(@html(
<pre><span style="font-weight: bold; color: rgb(0, 0, 0);">procedure</span><span
 style="color: rgb(0, 0, 0);"> rex(</span><span
 style="font-weight: bold; color: rgb(0, 0, 0);">var</span><span
 style="color: rgb(0, 0, 0);"> r:mtx; rx:mtx); </span><span
 style="font-weight: bold; color: rgb(0, 0, 128);">inline</span><span
 style="color: rgb(0, 0, 0);">;</span></pre>
  <a href="ex/rex0.pas.html">Example</a>
))}
procedure rex(var r:mtx; rx:mtx); inline;
var i:longint;
begin
  for i := 1 to lgn(rx) do begin
    rex(r,elm(rx,i,1),elm(rx,i,2))
  end
end;

var
  roi_: mtx; // roi workspace;

{**
@abstract(@html(
<pre>
  <span style="font-weight: bold; color: rgb(0, 0, 0);">procedure</span><span
  style="color: rgb(0, 0, 0);"> roi(</span><span
  style="font-weight: bold; color: rgb(0, 0, 0);">var</span><span
  style="color: rgb(0, 0, 0);"> m:mtx; ri:mtx);</span>
</pre>
<br>
Trim a table of data in one or several regions of interest.
This function operates on tables (stored in matrix m) whose first column must be
a monotonously increasing independant variable. Each row of ri (a 2 column
matrix) defines a low and high bracket. Rows of m, whose independant variable
(1st col.) is in none of the given brackets are discarded from m. The size
of m will therefore eventually decrease.<br> Regions of interest may overlap,
however the first column of ri must increase monotonously.

<a href="ex/roi0.pas.html">Example</a>
))}
procedure roi(var m:mtx; ri:mtx);
var
  j: longint;  // region loop index
  l: longint;  // # count rows processed so far
  c: longint;  // # of columns of m
  r: longint;  // # of rows (i.e. rois) of ri
  lj: longint; // row index of first row of jeme roi
  hj: longint; // row index of last row of jeme roi
  bj: longint; // m.e[bj] is m(lj,1) move starting point and...
  dj: longint; // destination in m.e[]
  nl: longint; // # of rows in roi
  nj: longint; // # of elements to shift
  r1, r2: double; // low high of current roi in loop
begin
  if mty(ri) then exit;
  assert(col(ri)=2);
  c := col(m);
  r := lgn(ri);
  zer(roi_,r,3);
  l := 0;
  r2 := -Infinity;
  for j := 1 to r do begin
    r1 := elm(ri,j,1);
    if r1<r2 then r1:= r2; // fix region overlap
    r2 := elm(ri,j,2);
    if r1<r2 then begin
      lj := fnd(m,r1);
      hj := fnd(m,r2);
      dec(hj);
      assert(hj<=lgn(m));
      nl := (hj-lj+1);
      nj := nl*c;
      if nl>0 then begin
        bj := idx(m,lj,1);
        dj := idx(m,l+1,1);
        //writeln('j':5,'l':5,'lj':5,'hj':5,'bj':5,'nl':5,'nj':5,'dj':5);
        //writeln(j:5,l:5,lj:5,hj:5,bj:5,nl:5,nj:5,dj:5);
        l := l + nl // increment AFTER calc. of dj
      end
    end else begin // empty region
      nj := 0;
      bj := 0;
      dj := 0
    end;
    elm(roi_,j,1,bj);
    elm(roi_,j,2,dj);
    elm(roi_,j,3,nj)
  end;
  for j := 1 to r do begin
    bj := round(elm(roi_,j,1));
    dj := round(elm(roi_,j,2));
    nj := round(elm(roi_,j,3));
    if nj>0 then move(m.e[bj],m.e[dj],nj*sizeof(double))
  end;
  dim(m,l,col(m))
end;

{**
@abstract(@html(
<pre>
  <span style="font-weight: bold; color: rgb(0, 0, 0);">procedure</span><span
  style="color: rgb(0, 0, 0);"> rix(</span><span
  style="font-weight: bold; color: rgb(0, 0, 0);">var</span><span
  style="color: rgb(0, 0, 0);"> r:mtx; ri, rx:mtx); </span><span
  style="font-weight: bold; color: rgb(0, 0, 128);">inline</span><span
  style="color: rgb(0, 0, 0);">;</span>
</pre>
<a href="ex/rix0.pas.html">Example</a>
))}
procedure rix(var r:mtx; ri, rx:mtx); inline;
begin
  roi(r,ri);
  rex(r,rx)
end;

end.
