# Block matrices package for Maple 6
# Copyright (c) 2001 Francis J. Wright

BlockMatrices := module()
   local init, BlockMatrixInverse;
   export BlockMatrix, `+`, `-`, `.`, `/`, `*`, evalb;
   option package, load=init, `Copyright (c) 2001 Francis J. Wright`;
   description "Prototype support for block matrices";

   init := proc()
      # Types MUST be global!
      global `type/BlockMatrix`, `type/BlockMatrixSetTypes`;
      `type/BlockMatrix` := '`module`(get, getRowp, getColp, set)';
      `type/BlockMatrixSetTypes` :=
         '{BlockMatrix,matrix,list(list),symbol,0,1}'
   end proc;

   init();  # run also when module defined

   # Constructor function:
   BlockMatrix := proc(Elements::BlockMatrixSetTypes,
                       Rowp::list(nonnegint),
                       Colp::list(nonnegint))
      local m;
      m := module()
         local ranges, elements, rowp, colp, unblockIndex;
         export get, getRowp, getColp, set,
                getElements, getElement, setElements, setElement, clone,
                isDiagonal, isUpperTriangular, isLowerTriangular,
                matadd, negate, inverse, scalarmul, multiply,
                equal, iszero, isunit, transpose, exponential;
         description "Block matrix";

         ranges := proc(P)  # P::list(posint)
            # Convert a partition (list of submatrix dimensions)
            # to a list of matrix sub-ranges.
            local beg, Rs, r, nxt;
            beg := 1;  Rs := NULL;
            for r in P do
               nxt := :-`+`(beg,r);
               Rs := Rs, beg..nxt-1;
               beg := nxt
            end do;
            [Rs] # [ ] in case of single range!
         end proc;

         get := proc()
            # Return a matrix of matrices
            local rows, cols;
             if not assigned(elements) then
               error "Block matrix not initialized!"
            end if;
            rows := ranges(rowp);
            cols := ranges(colp);
            matrix([seq([seq(linalg[submatrix](elements,i,j),
               j=cols)],i=rows)])
         end proc;

         getRowp := ()->
            # Return the row partition list
            if not assigned(rowp) then
               error "Block matrix not initialized!"
            else rowp end if;

         getColp := ()->
            # Return the column partition list
            if not assigned(colp) then
               error "Block matrix not initialized!"
            else colp end if;

         set := proc(Elements::BlockMatrixSetTypes,
                       Rowp::list(nonnegint),
                       Colp::list(nonnegint))
            # Set the block matrix data flexibly

            local unblock;

            unblock := proc(LLM)
               # LLM is a list of lists of matrices
               # Return: elements, rowp, colp
               linalg[blockmatrix](                 # elements
                  nops(LLM), nops(LLM[1]), # rows, cols
                  map(op, LLM)),  # flat list of matrices
               map(linalg[rowdim], map2(op,1,LLM)), # rowp
               map(linalg[coldim], LLM[1])          # colp
            end proc;

            unassign('rowp', 'colp');

            if Elements::'BlockMatrix' then
               elements, rowp, colp :=
                  Elements:-getElements(),
                     Elements:-getRowp(), Elements:-getColp()
            elif Elements::'matrix(matrix)' then
               elements, rowp, colp :=
                  unblock(convert(Elements,listlist))
            elif Elements::'list(list(matrix))' then
               elements, rowp, colp :=
                  unblock(Elements)
            elif Elements::'matrix' then
               elements := copy(Elements)
            elif Elements::'list(list)' then
               elements := matrix(Elements)
            elif Elements = 1 then  # unit matrix
               if nargs < 2 then
                  error "Must specify at least one dimension for unit block matrix!"
               elif nargs = 3 and Rowp <> Colp then
                  error "Unit block matrix must be square!"
               end if;
               elements := array(identity,
                  1..:-`+`(op(Rowp)), 1..:-`+`(op(Rowp)));
               rowp := Rowp;  colp := Rowp;
               return  # nothing
            else  # Elements::{symbol, 0}
               if nargs < 3 then
                  error "Must specify dimensions for symbolic or zero block matrix!"
               end if;
               if Elements::symbol then
                  elements := matrix(:-`+`(op(Rowp)),
                     :-`+`(op(Colp)), (i,j)->Elements[i,j])
               else  # Elements = 0
                  elements := array(sparse,
                     1..:-`+`(op(Rowp)), 1..:-`+`(op(Colp)))
               end if;
               rowp := Rowp;  colp := Colp;
               return  # nothing
            end if;

            if nargs > 1 then
               if assigned(rowp) then
                  if rowp <> Rowp then
                     error "Incompatible initialization data!"
                  end if
               elif linalg[rowdim](Elements) <> :-`+`(op(Rowp)) then
                  error "Incompatible initialization data!"
               else
                  rowp := Rowp
               end if
            end if;

            if nargs > 2 then
               if assigned(colp) then
                  if colp <> Colp then
                     error "Incompatible initialization data!"
                  end if
               elif linalg[coldim](Elements) <> :-`+`(op(Colp)) then
                  error "Incompatible initialization data!"
               else
                  colp := Colp
               end if
            end if;

            if not(assigned(elements) and
                   assigned(rowp) and assigned(colp)) then
               error "Must initialize fully or not at all!"
            end if;

         end proc;

         unblockIndex := (index,partition)->
            if index::list then  # unblock
               :-`+`(op(1..index[1]-1,partition),index[2])
            else index end if;

         getElements := ()->
            # Return all the elements as a (flat) matrix
            if not assigned(elements) then
               error "Block matrix not initialized!"
            else elements end if;

         getElement := proc(a::[{posint,[posint,posint]},
                                {posint,[posint,posint]}])
            # Return one element
            if not assigned(elements) then
               error "Block matrix not initialized!"
            end if;
            elements[unblockIndex(a[1],rowp),
                     unblockIndex(a[2],colp)]
         end proc;

         setElements := proc()
            # Set one or more elements
            local arg;
            for arg in args do setElement(arg) end do;
            return  # nothing
         end proc;

         setElement := proc(a::[{posint,[posint,posint]},
                                {posint,[posint,posint]}]=algebraic)
            # Set one element
            elements[unblockIndex(lhs(a)[1],rowp),
                     unblockIndex(lhs(a)[2],colp)] := rhs(a)
         end proc;

         clone := ()->
            # Return copy (clone) of block matrix
            BlockMatrix(copy(elements), rowp, colp);

         isDiagonal := proc()
            # Return true if the matrix is block diagonal
            local rngs, i, j, ii, jj;
            if rowp <> colp then return false end if;
            rngs := ranges(rowp);
            for i to nops(rngs) do
               for j to nops(rngs) do
                  if i <> j then
                     # Check submatrix is zero:
                     for ii in rngs[i] do
                        for jj in rngs[j] do
                           if elements[ii,jj] <> 0 then return false end if
                        end do
                     end do
                  end if
               end do
            end do;
            true
         end proc;

         isUpperTriangular := proc()
            # Return true if the matrix is block upper triangular
            local rngs, i, j, ii, jj;
            if rowp <> colp then return false end if;
            rngs := ranges(rowp);
            for i from 2 to nops(rngs) do
               for j to i-1 do
                  # Check submatrix is zero:
                  for ii in rngs[i] do
                     for jj in rngs[j] do
                        if elements[ii,jj] <> 0 then return false end if
                     end do
                  end do
               end do
            end do;
            true
         end proc;

         isLowerTriangular := proc()
            # Return true if the matrix is block lower triangular
            local rngs, i, j, ii, jj;
            if rowp <> colp then return false end if;
            rngs := ranges(rowp);
            for j from 2 to nops(rngs) do
               for i to j-1 do
                  # Check submatrix is zero:
                  for ii in rngs[i] do
                     for jj in rngs[j] do
                        if elements[ii,jj] <> 0 then return false end if
                     end do
                  end do
               end do
            end do;
            true
         end proc;

         matadd := (M::BlockMatrix)->
            # Add block matrix M
            if rowp = M:-getRowp() and colp = M:-getColp() then
               BlockMatrix(linalg[:-matadd](elements, M:-getElements()),
                  rowp, colp)
            else
               error "Block matrices not conformable for addition!"
            end if;

         negate := ()->
            # Return negative of block matrix
            BlockMatrix(map(:-`-`, elements), rowp, colp);

         inverse := ()->
            # Return inverse of block matrix
            if rowp = colp then
               BlockMatrix(linalg[:-inverse](elements), rowp, colp)
            else
               error "Only square block matrices are invertible!"
            end if;

         scalarmul := (k::algebraic)->
            # Multiply by a scalar k
            BlockMatrix(linalg[:-scalarmul](elements, k), rowp, colp);

         multiply := (M::BlockMatrix)->
            # Multiply by a block matrix M
            if colp = M:-getRowp() then
               linalg[:-multiply](elements, M:-getElements());
               BlockMatrix(simplify(%), rowp, M:-getColp())
            else
               error "Block matrices not conformable for multiplication!"
            end if;

         equal := (M::BlockMatrix)->
            # Return true if equal to block matrix M
            rowp = M:-getRowp() and colp = M:-getColp() and
            linalg[:-equal](elements, M:-getElements());

         iszero := ()->
            # Return true if zero (block) matrix
            linalg[:-iszero](elements);

         isunit := ()->
            # Return true if unit block matrix
            rowp = colp and
            linalg[:-iszero](evalm(elements-&*()));

         transpose := ()->
            # Return transpose of block matrix
            BlockMatrix(linalg[:-transpose](elements), colp, rowp);

         exponential := (t::algebraic)->
            # Return exponential of block matrix
            if rowp = colp then
               BlockMatrix(linalg[:-exponential](elements,args), rowp, colp)
            else
               error "Only square block matrices are exponentiable!"
            end if;

      end module;

      if nargs > 0 then m:-set(args) end if;
      eval(m)
   end proc;

   # Overload arithmetic operators:

   `+` := proc(A,B)
      if A::'BlockMatrix' then
         if B::'BlockMatrix' then
            return A:-matadd(B)
         end if
      elif not B::'BlockMatrix' then
         return :-`+`(A,B)
      end if;
      error "Can add block matrices only to each other!"
   end proc;

   `-` := proc(A)
      if A::'BlockMatrix' then
         A:-negate()
      else
         :-`-`(A)
      end if
   end proc;

   `.` := proc(A,B)
      if A::'BlockMatrix' then
         if B::'BlockMatrix' then
            return A:-multiply(B)
         elif B::'algebraic' then
            return A:-scalarmul(B)
         end if
      elif B::'BlockMatrix' then
         if A::'algebraic' then
            return B:-scalarmul(A)
         end if
      else
         return :-`.`(A,B)
      end if;
      error "Can multiply block matrices only with each other and scalars!"
   end proc;

   `/` := proc(A)
      if A::'BlockMatrix' then
         BlockMatrixInverse(A)  # inert!
      else
         :-`/`(1,A)  # binary!
      end if
   end proc;

   `*` := proc(A,B)
      # Disallow BM*BM.
      # Allow only BM*k, BM/BM => BM*(/BM), k*BM, k/BM => k*(/BM).
      if A::'BlockMatrix' then
         if B::'BlockMatrix' then
            error "Must use . to multiply block matrices with each other!"
         elif B::specfunc(anything,BlockMatrixInverse) then
            A:-multiply(op(B):-inverse())
         elif B::'algebraic' then
            A:-scalarmul(B)
         else
            error "Can only multiply/divide block matrices and scalars by each other!"
         end if
      elif B::'BlockMatrix' then
         if A::'algebraic' then
            B:-scalarmul(A)
         else
            error "Can use * only to multiply block matrices by scalars!"
         end if
      elif B::specfunc(anything,BlockMatrixInverse) then
         if A::'algebraic' then
            op(B):-inverse():-scalarmul(A)
         else
            error "Can only divide block matrices and scalars by each other!"
         end if
      else
         :-`*`(A,B)
      end if
   end proc;

   evalb := proc()
      local A, B;
      if type(args, 'equation') then
         A := lhs(args);  B := rhs(args);
         if not A::'BlockMatrix' then
            A, B := B, A
         end if;
         if A::'BlockMatrix' then
            if B::'BlockMatrix' then
               return A:-equal(B)
            elif B = 0 then
               return A:-iszero()
            elif B = 1 then
               return A:-isunit()
            else
               return false
            end if
         end if
      end if;
      :-evalb(args)
      end proc;

end module;

# For Emacs:
# Local Variables:
# mode: text
# eval: (auto-fill-mode nil)
# indent-tabs-mode: nil
# End: