# The linalg package for Linear Algebra

# initiated September 04, 2008

# Change History

# v1.0.0 as of August 01, 2010:
#   - initial release
# v1.0.1 as of December 23, 2010:
#   - internal migration of the code to the new substring operator.
# v1.0.2 as of February 13, 2010:
#   - improved error handling of `linalg.hilbert` and `linalg.matrix`.
# v1.1.0 as of September 30, 2012
#   - The new function `linalg.norm` returns the norm of a matrix or vector.
#   - Improved error handling of `linalg.checkmatrix` and `linalg.matrix`.
# v1.1.1 as of February 21, 2014
#   - See change log for Agena 2.1.3
# v1.2.0 as of March 01, 2014
#   - See change log for Agena 2.1.4

package.checkclib('linalg');  # has C dynamic link library been found ?

linalg.initstring := 'linalg v1.2.4 as of July 16, 2014\n';

create table linalg.aux;

# define metamethods for vectors (do not change the name of this metatable for vectors)
linalg.vmt := [
   '__add' ~ proc(a, b) is  # the Agena version is 7 % slower
      return linalg.add(a, b)
   end,
   '__sub' ~ proc(a, b) is  # the Agena version is 7 % slower
      return linalg.sub(a, b)
   end,
   '__mul' ~ proc(n, A) is
      return linalg.scalarmul(n, A)
   end,
   '__unm' ~ << a -> (-1)*a >>,
   # the `sqrt(qsadd(a))` and `sqrt(linalg.dotprod(v, v))` versions are much slower
   '__abs' ~ proc(v) is  # this Agena version is 28 % faster than the C version
      local s := 0;
      for i to v.dim do
         inc s, v[i]**2
      od;
      return sqrt(s)
   end,
   '__qsadd' ~ proc(v) is
      local s := 0;
      for i to v.dim do
         inc s, v[i]**2
      od;
      return s
   end,
   '__aeq' ~ proc(A, B) is  # added 2.1.3
      return linalg.vaeq(A, B)
   end,
   '__eeq' ~ proc(A, B) is  # added 2.1.4
      return linalg.veeq(A, B)
   end,
   '__tostring' ~ proc(v) is
      if type(v.dim) <> number then
         error('Error in `print`: internal vector dimension of incorrect type, must be a number.')
      fi;
      local str := '[ ';
      for i to v.dim do  # concat works on arrays only, but vectors are dictionaries
         str := str & v[i] & ', '
      od;
      return str[1 to -3] & ' ]'
   end,
   '__writeindex' ~ proc(v, key, value) is linalg.setvelem(v, key, value) end,
   '__index' ~ proc(v, n) is
      if n < 1 or n > v.dim then
         error('Error via `linalg.vmt` metatable: index out of range.')
      elif rawget(v, n) = null then
         return 0
      else
        return rawget(v, n)
      fi
   end
];


# define metamethods for matrices (do not change the name of this metatable for matrices)
linalg.mmt := [
   '__add' ~ proc(A, B) is
      local rowdim, coldim, R;
      linalg.checkmatrix(A, B);  # both matrices must be checked
      if A.dim <> B.dim then
         error('Error via `linalg.mmt` metatable: matrix dimensions incompatible.')
      fi;
      rowdim, coldim := left(A.dim), right(A.dim);
      R := linalg.matrix(rowdim, coldim);
      for i to rowdim do  # traverse row vectors
         for j to coldim do  # traverse components, faster than calling linalg.scalarmul
            R[i, j] := A[i, j] + B[i, j]
         od
      od;
      return R
   end,
   '__sub' ~ proc(A, B) is
      return A + (-1)*B
   end,
   '__mul' ~ proc(a, b) is  # extended 2.1.3
      if b :: number then
         a, b := b, a;
      fi;
      if a :: number and b :: 'matrix' then
         local rowdim, coldim, R;
         rowdim, coldim := left(b.dim), right(b.dim);
         R := linalg.matrix(rowdim, coldim);
         for i to rowdim do  # traverse row vectors
            for j to coldim do  # traverse components, faster than calling linalg.scalarmul
               R[i, j] := b[i, j] * a
            od
         od;
         return R
      elif a :: 'matrix' and b :: 'matrix' then
         return linalg.mmul(a, b)
      else
         error('Error in `*`: wrong kind of arguments.')
      fi
   end,
   '__aeq' ~ proc(A, B) is  # added 2.1.4
      return linalg.maeq(A, B)
   end,
   '__eeq' ~ proc(A, B) is  # added 2.1.3
      return linalg.meeq(A, B)
   end,
   '__index' ~ proc(A, n) is
      if n = null or n < 1 or n > left(A.dim) then
         error('Error via `linalg.mmt` metatable: index out of range or null.')
      elif rawget(A, n) = null then
         return 0
      else
         return rawget(A, n)
      fi
   end,
   # a __writeindex metamethod is not needed since the vector __index metamethod prevents out-of-range assignments
   '__tostring' ~ proc(m) is
      local maxlen, str, rowdim, coldim;
      rowdim := linalg.rowdim(m);
      coldim := linalg.coldim(m);
      # vector which holds the length of the longest number in each row
      maxlen := linalg.zero(coldim);
      # first determine the `longest number` column by column
      for i to coldim do
         for j in tostring @ linalg.column(m, i) do
            if size j > maxlen[i] then maxlen[i] := size j fi
         od
      od;
      # now build matrix output
      str := '';
      for i to rowdim do
         str := str & '[ ';
         for j to coldim do  # traverse rows
            str := str & strings.format('%' & maxlen[j] & 's', m[i, j]) & ', '
         od;
         str := str[1 to -3] & ' ]\n'
      od;
      return str[1 to -2]  # return string without trailing newline
   end
];


# matrix constructor, 28.07.2008; changed 06.09.2008; improved & extended 15.12.2008 .. 19.12.2008; improved 30.09.2012
# accepts sequences, tables, or vectors (all indicating row vectors, even empty ones !)
linalg.matrix := proc(?) is
   local A, r, rowdim, coldim;
   if nargs >= 2 and varargs[1] :: number and varargs[2] :: number then
      rowdim := varargs[1];
      coldim := varargs[2];
      create table r(rowdim);
      if nargs = 2 or varargs[3] = [] then  # create a sparse matrix
         to rowdim do
            insert linalg.vector(coldim, []) into r
         od
      elif varargs[3] :: table then  # process Maple-like syntax
         local rowv, tbl, c;
         tbl := varargs[3];
         c := 0;
         for i to rowdim do
            rowv := linalg.vector(coldim, []);
            for j to coldim do
               inc c;
               if isstructure(tbl[c]) then
                  error('Error in `linalg.matrix`: structures not allowed.')
               fi;
               rowv[j] := tbl[c]
            od;
            r[i] := rowv
         od
      else
         error('Error in `linalg.matrix`: table expected as third argument, got ' & type(varargs[3]) & '.')
      fi;
   else  # create a `filled` matrix
      rowdim := nargs;
      if rowdim = 0 then  # Agena 1.8.0 u1
         error('Error in `linalg.matrix`: expected at least one argument.')
      fi;
      create table r(rowdim);
      for i in varargs do
         A := i;  # do not change control variable !
         if typeof(A) in {table, sequence} then  # convert to sequence
            A := linalg.vector(unpack(A))
         else
            linalg.checkvector(A)
         fi;
         insert A into r
      od;
      coldim := r[1].dim;
      # check whether all row vectors have the same size (dimension)
      for i to size(r)-1 do
         if r[i].dim <> r[i+1].dim then
            error('Error in `linalg.matrix`: row vectors of different dimensions.')
         fi
      od
   fi;
   settype(r, 'matrix');
   setmetatable(r, linalg.mmt);
   r.dim := rowdim:coldim;
   return r
end;


# returns the colum dimension of a matrix; 27.07.2008; extended 06.09.2008; simplified 18.12.2008
linalg.coldim := proc(A) is
   return right(linalg.checkmatrix(A, true))
end;


# returns the row dimension of the matrix A; 27.07.2008; extended 06.09.2008; simplified 18.12.2008
linalg.rowdim := proc(A) is
   return left(linalg.checkmatrix(A, true))
end;


# determines the dimension of a matrix or a vector. If A is a matrix, the result is a pair
# with the left-hand side the number of rows and the right-hand side the number of columns.
# If A is a vector, the size of the vector is determined. 18.12.2008, 30.09.2012
linalg.dim := proc(A) is
   if not(linalg.ismatrix(A) or linalg.isvector(A)) then
      error('Error in `linalg.dim`: matrix or vector expected, got ' & type(A) & '.')
   fi;
   return A.dim
end;


# checks whether the matrix A is a square matrix, 27.07.2008; simplified 18.12.2008
linalg.issquare := proc(A) is
   local a := linalg.checkmatrix(A, true);
   return left(a) = right(a)
end;


# checks whether a matrix is a square matrix and issues an error if not, 27.07.2008;
# patched & simplified 18.12.2008
linalg.checksquare := proc(A, option) is
   local dims := linalg.checkmatrix(A, true);
   if left(dims) <> right(dims) then
      error('Error in `linalg.checksquare`: square matrix expected.')
   elif option = true then
      return dims
   fi
end;


# return the n-th column of a matrix or row vector o as a new vector; 18.12.2008
linalg.column := proc(o, n :: number) is
   local R;
   if float(n) then
      error('Error in `linalg.column`: only positive integers accepted, got ' & n & '.')
   fi;
   case typeof(o)
      of 'matrix' then
         local rowdim;
         if n < 1 or n > right(o.dim) then  # not(n in 1:right(o.dim)) is slower
            error('Error in `linalg.column`: column ' & n & ' does not exist.')
         fi;
         rowdim := left(o.dim);
         R := linalg.vector(rowdim, []);
         for i to rowdim do
            R[i] := o[i, n]
         od
      of 'vector' then
         if n < 1 or n > o.dim then
            error('Error in `linalg.column`: column ' & n & ' does not exist.')
         fi;
         R := linalg.vector(o[n])
      else
         error('Error in `linalg.column`: matrix or vector expected, got ' & typeof(o) & '.')
   esac;
   return R
end;


# Create a generalized n x n Hilbert matrix, it has 1/(i+j-x) as its (i,j)-th entry.
# If x is not specified, then x is 1. Written 18.12.2008, improved 1.3.3, January 31, 2011.
linalg.hilbert := proc(n :: number, x) is
   local R;
   if isnonposint(n) then
      error('Error in `linalg.hilbert`: posint expected as first argument, got ' & tostring(x) & '.')
   fi;
   if x = null then
      x := 1
   elif x :- number then
      error('Error in `linalg.hilbert`: number expected as second argument, got ' & type(x) & '.')
   fi;
   R := linalg.matrix(n, n);
   for i to n do
      for j to n do
         R[i, j] := 1/(i+j-x)
      od
   od;
   return R
end;


# Compute the cross-product of two vectors of dimension 3. 20.12.2008
linalg.crossprod := proc(u, v) is
   if linalg.checkvector(u, v) <> 3 then
      error('Error in `linalg.crossprod`: vectors of dimension 3 expected.')
   fi;
   return linalg.vector(u[2]*v[3]-u[3]*v[2], u[3]*v[1]-u[1]*v[3], u[1]*v[2]-u[2]*v[1])
end;


# Specialized version of map for matrices. 19.12.2008
linalg.mmap := proc(f :: procedure, A :: matrix, ?) is
   local R, dims, rowdim;
   dims := linalg.checkmatrix(A, true);
   rowdim := left(dims);
   R := linalg.matrix(rowdim, right(dims));
   for i to rowdim do
      R[i] := linalg.vmap(f, A[i], unpack(varargs))
   od;
   return R
end;


# Specialized version of zip for matrices. 10.01.2010, tweaked 1.3.3, January 31, 2011
linalg.mzip := proc(f :: procedure, A :: matrix, B :: matrix, ?) is
   local R, rowdim;
   if A.dim <> B.dim then
      error('Error in `linalg.mzip`: matrices must have equal dimensions.')
   fi;
   rowdim := left(A.dim);
   R := linalg.matrix(rowdim, right(A.dim));
   for i to rowdim do
      R[i] := linalg.vzip(f, A[i], B[i], unpack(varargs))
   od;
   return R
end;


# checks whether the matrix A is a diagonal matrix. If so, it returns true and false otherwise.
# Maple version March 17, 1998, Agena version December 21, 2008.
linalg.isdiagonal := proc(A) is
   local ranges, rowdim, coldim;
   ranges := linalg.checkmatrix(A, true);
   rowdim := left(ranges);
   coldim := right(ranges);
   if rowdim <> coldim then
      return false
   fi;
   for i to rowdim do
      for j to coldim do
         skip when i = j;
         if A[i, j] <> 0 then
            return false
         fi
      od
   od;
   return true
end;


# checks whether the matrix A is an identity matrix. If so, it returns true and false
# otherwise. Maple version March 16, 1998, Agena version December 21, 2008.
linalg.isidentity := proc(A) is
   local rowdim, n, i, j;
   rowdim := left(linalg.checkmatrix(A, true));
   if not(linalg.isdiagonal(A)) then
      return false
   fi;
   for i to rowdim do
      if A[i, i] <> 1 then return false fi
   od;
   return true
end;


# checks whether the matrix A is a symmetric matrix. If so, it returns true and false otherwise.
# Maple version March 16, 1998, Agena version December 21, 2008
linalg.issymmetric := proc(A) is
   local ranges, i, j, mrange, nrange;
   ranges := linalg.checkmatrix(A, true);
   if not(linalg.issquare(A)) then
      return false
   fi;
   mrange := left(ranges);
   nrange := right(ranges);
   for m from 1 to mrange-1 do
      for n from m + 1 to nrange do
         skip when m = n;
         if A[m, n] <> A[n, m] then
            return false
         fi
      od
   od;
   return true
end;


# checks whether the matrix A is an antisymmetric matrix. If so, it returns true and false otherwise.
# Maple version March 20, 1998, Agena version December 21, 2008
linalg.isantisymmetric := proc(A) is
   local ranges, i, j, mrange, nrange;
   ranges := linalg.checkmatrix(A, true);
   if not(linalg.issquare(A)) then
      return false
   fi;
   mrange := left(ranges);
   nrange := right(ranges);
   for m from 1 to mrange-1 do
      for n from m + 1 to nrange do
         skip when m = n;
         if A[m, n] <> -A[n, m] then
            return false
         fi
      od
   od;
   return true
end;


# checks whether vectors or meatrices have been passed to linalg.augment
linalg.aux.augment_check := proc(x) is
   local utype := typeof(x);
   if not(utype in {'matrix', 'vector'}) then
      error('Error in `linalg.augment`: matrix or vector expected.')
   fi;
   if utype = 'vector' then
      return x.dim, 1
   else
      return left(x.dim), right(x.dim)
   fi
end;

# inserts a vector or matrix into the resulting matrix
linalg.aux.augment_insert := proc(x, rowdim, coldim, o) is
   if typeof(x) = 'vector' then
      for i to rowdim do
         insert x[i] into o[i]
      od
   else
      for i to rowdim do
         for j to coldim do
            insert x[i, j] into o[i]
         od
      od
   fi
end;

# Joins two or more matrices or vectors together horizontally. Vectors are supposed to be
# column vectors. The matrices and vectors must have the same number of rows.
# January 11, 2010, 0.30.2
linalg.augment := proc(?) is
   local r, s, rowdim, coldim, coldims, oldrowdim;
   # get dimensions of first object
   oldrowdim, coldim := linalg.aux.augment_check(varargs[1]);
   coldims := coldim;
   # create matrix (in raw form)
   create table r;
   for i to oldrowdim do
      insert [] into r;
   od;
   linalg.aux.augment_insert(varargs[1], oldrowdim, coldim, r);
   # get dimensions of all other object
   for i from 2 to nargs do
      rowdim, coldim := linalg.aux.augment_check(varargs[i]);
      if oldrowdim <> rowdim then
         error('Error in `linalg.augment`: objects must have the same number of rows.')
      else
         oldrowdim := rowdim;
         inc coldims, coldim
      fi;
      linalg.aux.augment_insert(varargs[i], oldrowdim, coldim, r)
   od;
   # set attributes to row vectors
   s := size(r[1]);
   for i in r do
      i.dim := s;
      settype(i, 'vector');
      setmetatable(i, linalg.vmt)
   od;
   # set matrix attributes
   if rowdim = null then  # only one argument given ?
      rowdim := oldrowdim
   fi
   r.dim := rowdim:coldims;
   settype(r, 'matrix');
   setmetatable(r, linalg.mmt);
   return r
end;


linalg.aux.stack_check := proc(x) is
   local utype := typeof(x);
   if not(utype in {'matrix', 'vector'}) then
      error('Error in `linalg.stack`: matrix or vector expected.')
   fi;
   if utype = 'vector' then
      return 1, x.dim
   else
      return left(x.dim), right(x.dim)
   fi
end;


linalg.aux.stack_insert := proc(x, rowdim, o) is
   if typeof(x) = 'vector' then
      insert x into o
   else
      for i to rowdim do
         insert x[i] into o
      od
   fi
end;

# Joins two or more matrices or vectors together vertically. Vectors are supposed to be
# row vectors. The matrices and vectors must have the same number of columns.
# January 11, 2010, 0.30.2
linalg.stack := proc(?) is
   local r, s, rowdim, rowdims, coldim, oldcoldim;
   # get dimensions of first object
   rowdim, oldcoldim := linalg.aux.stack_check(varargs[1]);
   rowdims := rowdim;
   # create matrix (in raw form)
   create table r;
   linalg.aux.stack_insert(varargs[1], rowdim, r);
   # get dimensions of all other object
   for i from 2 to nargs do
      rowdim, coldim := linalg.aux.stack_check(varargs[i]);
      if oldcoldim <> coldim then
         error('Error in `linalg.stack`: objects must have the same number of columns.')
      else
         oldcoldim := coldim;
         inc rowdims, rowdim
      fi;
      linalg.aux.stack_insert(varargs[i], rowdim, r)
   od;
   # set matrix attributes
   if coldim = null then  # only one argument given ?
      coldim := oldcoldim
   fi
   r.dim := rowdims:coldim;
   settype(r, 'matrix');
   setmetatable(r, linalg.mmt);
   return r
end;


# swap rows in a matrix, January 12, 2010, Agena 0.30.0
linalg.swaprow := proc(A :: matrix, p :: number, q :: number) is
   local B, rowdim;
   rowdim := left(A.dim);
   if p < 1 or p > rowdim or float(p) then
      error('Error in `linalg.swaprow`: second argument must be a posint <= row dimension.')
   fi;
   if q < 1 or q > rowdim or float(p) then
      error('Error in `linalg.swaprow`: third argument must be a posint <= row dimension.')
   fi;
   B := copy(A);
   B[q], B[p] := B[p], B[q];  # optimised 2.1.3
   return B
end;


# swap columns in a matrix, January 12, 2010, Agena 0.30.0
linalg.swapcol := proc(A :: matrix, p :: number, q :: number) is
   local B, rowdim, coldim;
   rowdim, coldim := left(A.dim), right(A.dim);
   if p < 1 or p > coldim or float(p) then
      error('Error in `linalg.swapcol`: second argument must be a posint <= column dimension.')
   fi;
   if q < 1 or q > coldim or float(p) then
      error('Error in `linalg.swapcol`: third argument must be a posint <= column dimension.')
   fi;
   B := copy(A);
   for i to rowdim do
      B[i, p], B[i, q] := B[i, q], B[i, p]  # optimised 2.1.3
   od;
   return B
end;


# Returns the norm of a matrix or vector, Agena 1.8.0 u1

linalg.norm := proc(A, n) is
   if linalg.ismatrix(A) then
      create local seq r(A.dim[1]);
      for i to A.dim[1] do
         insert sadd(<< x -> abs x >> @ A[i]) into r
      od;
      return max(r)
   elif linalg.isvector(A) then
      if unassigned(n) then
         n := infinity
      elif not(isposint(n) or n = infinity) then
         error('Error in `linalg.norm`: expected a posint as second argument, got ' & type(n) & '.')
      fi;
      if n <> infinity then
         return root(sadd(linalg.vmap(<< x -> abs(x)**n >>, A)), n)
      else
         return max(linalg.vmap(<< x -> abs(x) >>, A))
      fi
   else
      error('Error in `linalg.norm`: matrix or vector expected, got ' & type(A) & '.')
   fi
end;


#/

Credit: Agena's `linalg.scale` is a port of the ALGOL 60 function REASCL, being part of the NUMAL package,
originally published by The Stichting Centrum Wiskunde & Informatica, Amsterdam, The Netherlands.

The Stichting Centrum Wiskunde & Informatica (Stichting CWI) (legal successor of Stichting Mathematisch
Centrum) at Amsterdam has granted permission to Paul McJones to attach the integral NUMAL library manual
to his software preservation project web page.

( URL: http://www.softwarepreservation.org/projects/ALGOL/applications/ )

It may be freely used. It may be copied provided that the name NUMAL and the attribution to the Stichting
CWI are retained.

Original ALGOL 60 credits to REASCL:

AUTHORS  : T.J. DEKKER, W. HOFFMANN.
CONTRIBUTORS: W. HOFFMANN, S.P.N. VAN KAMPEN.
INSTITUTE: MATHEMATICAL CENTRE.
RECEIVED: 731030.

BRIEF DESCRIPTION:

The procedure REASCL (scale) normalises the (non-null) columns of a matrix
in such a way that, in each column, an element of maximum absolute value
equals 1. The normalised vectors are delivered in the corresponding columns
of the new matrix.

RUNNING TIME: PROPORTIONAL TO N * (N2 - N1 + 1).

LANGUAGE:   ALGOL 60.

METHOD AND PERFORMANCE: SEE REF [1].

REFERENCES:
     [1].T.J. DEKKER AND W. HOFFMANN.
         ALGOL 60 PROCEDURES IN NUMERICAL ALGEBRA, PART 2.
         MC TRACT 23, 1968, MATH. CENTR., AMSTERDAM.

SOURCE TEXT(S):

 CODE 34183;
     COMMENT MCA 2413;
     PROCEDURE REASCL(A, N, N1, N2); VALUE N, N1, N2;
     INTEGER N, N1, N2; ARRAY A;
     BEGIN INTEGER I, J; REAL S;
         FOR J:= N1 STEP 1 UNTIL N2 DO
         BEGIN S:= 0;
             FOR I:= 1 STEP 1 UNTIL N DO
             IF ABS(A[I,J]) > ABS(S) THEN S:= A[I,J];
             IF S ^= 0 THEN
             FOR I:= 1 STEP 1 UNTIL N DO A[I,J]:= A[I,J] / S
         END
     END REASCL;
         EOP

/#

linalg.scale := proc(M :: matrix) is
   local A, s, n;
   A := copy(M);
   n := A.dim[1];
   for j to A.dim[2] do
      s := 0;
      for i to n do
         if abs(A[i, j]) > abs(s) then s := A[i, j] fi
      od;
      if s <> 0 then
         for i to n do div A[i, j], s od
      fi
   od;
   return A
end;


# Backward compatibility

linalg.veq := linalg.vaeq;
linalg.meq := linalg.maeq;
linalg.backsubs := linalg.gsolve;

