From: Peter Broadbery <peterb>
Date: Mon, 17 Jun 96 16:07:21 BST
Received: from co.uk (nags8) by nags2.nag.co.uk (4.1/UK-2.1)
	id AA03227; Mon, 17 Jun 96 16:07:24 BST
To: ax-bugs

Subject: fixbug
By: PAB
Fixed: bug878.as
--* From postmaster%watson.vnet.ibm.com@yktvmv.watson.ibm.com  Thu Oct  6 16:19:14 1994
--* Received: from yktvmv-ob.watson.ibm.com by watson.ibm.com (AIX 3.2/UCB 5.64/930311)
--*           id AA16974; Thu, 6 Oct 1994 16:19:14 -0400
--* Received: from watson.vnet.ibm.com by yktvmv.watson.ibm.com (IBM VM SMTP V2R3)
--*    with BSMTP id 5887; Thu, 06 Oct 94 16:19:19 EDT
--* Received: from YKTVMV by watson.vnet.ibm.com with "VAGENT.V1.0"
--*           id <A.PETEB.NOTE.VAGENT2.9615.Oct.06.16:19:18.-0400>
--*           for asbugs@watson; Thu, 06 Oct 94 16:19:19 -0400
--* Received: from YKTVMV by watson.vnet.ibm.com with "VAGENT.V1.0"
--*           id 9611; Thu, 6 Oct 1994 16:19:18 EDT
--* Received: from watson.ibm.com by yktvmv.watson.ibm.com (IBM VM SMTP V2R3)
--*    with TCP; Thu, 06 Oct 94 16:19:17 EDT
--* Received: by watson.ibm.com (AIX 3.2/UCB 5.64/930311)
--*           id AA19271; Thu, 6 Oct 1994 16:19:08 -0400
--* Date: Thu, 6 Oct 1994 16:19:08 -0400
--* From: peteb@watson.ibm.com (Peter A. Broadbery)
--* Message-Id: <9410062019.AA19271@watson.ibm.com>
--* To: asbugs@watson.ibm.com
--* Subject: [2] bug: tfGetExportList

--@ Bug Number:  bug878.as 
--@ Fixed  by:  PAB   
--@ Tested by:  none.as 
--@ Summary:    Fixed by SSD sometime back. 

-- Command line: -laxiom
-- Version: current
-- Original bug file name: eigvec.as


#include "axiom.as"

S ==> Symbol;
V ==>Vector;
I ==> Integer;
M ==> Matrix;
L ==> List;
OP ==> OutputPackage;
OF ==> OutputForm;
M2 ==> MatrixCategoryFunctions2;
UP ==> UnivariatePolynomial;
EP ==> EigenPackage;
NNI ==> NonNegativeInteger;
FPI ==> Fraction Polynomial Integer;
SAE ==> SimpleAlgebraicExtension;
SEG ==> Segment;

+++ Author: Bill Naylor
+++ Date Created: Fri Sep 23rd
+++ Date Last Updated: Mon Oct 3rd
+++
+++ Related Constructors: RoesLinearisedJacobian,
+++                       MyEigenFunctions,
+++                       RoesSchemeGeneral,
+++                       AspB. etc.
+++
+++ Description:
+++ MyEigVec(M,x,cP) enables the production of a 'deconstructed' form of the eigen
+++ vectors of a matrix, This is done by traversing a Gaussian elimination procedure
+++ which is performed on a matrix in the SimpleAlgebraicExtension(Fraction Polynomial
+++ Integer,UP,cp) which is equivalent to M. At every step of the elimination, an
+++ expression is removed, and replaced by a generic SAE, this avoids complexity
+++ building up in the coefficients of the SAEs. N.B. for small matrices (< rank 3) or
+++ sparse matrices this is not an efficient procedure.

MyEigVec(m : Matrix Fraction Polynomial Integer,x : Symbol
    ,cp : UnivariatePolynomial(x,Fraction Polynomial Integer)) : with {

  myGauss : () -> L(L(L(L(FPI))));
  destructSAE : (UP(x,FPI),NNI) -> L(FPI);
  constructSAE : NNI -> SAE(FPI,UP(x,FPI),cp);
  makeTempCoeff : I -> S;

}
== add {

  import from I;
  import from OP;
  import from OF;
  import from SEG(I);

  local len:I := nrows(m)$M(FPI) pretend I;
--  x:S == "x"::S;
  local tmpVarCnt:I := 0$I;
--  local cp:FPI ==> (characteristicPolynomial(m,x)$EP(I))::FPI;

  import from UP(x,FPI);

--  local cpUP ==> cp :: UP(x,FPI);

  import from SAE(FPI,UP(x,FPI),cp);

  local m2:M(SAE(FPI,UP(x,FPI),cp));
--  SAE2 := SAE(FPI,UP(x,FPI),cpUP)
  id:M(SAE(FPI,UP(x,FPI),cp)) := matrix([[0$SAE(FPI,UP(x,FPI),cp) for i in (1$I..len)$SEG(I)]$L(SAE(FPI,UP(x,FPI),cp)) for j in (1$I..len)$SEG(I)]$L(L(SAE(FPI,UP(x,FPI),cp))))$M(SAE(FPI,UP(x,FPI),cp));
  for i in (1$I..len)$SEG(I) repeat {
    set!(id,i,i,1$SAE(FPI,UP(x,FPI),cp));
  }
  local gen:SAE(FPI,UP(x,FPI),cp):=generator()$SAE(FPI,UP(x,FPI),cp);
  local genUP:UP(x,FPI) := lift(gen)$SAE(FPI,UP(x,FPI),cp);
  local m2:M(SAE(FPI,UP(x,FPI),cp)) := map((e:FPI):SAE(FPI,UP(x,FPI),cp) +->
coerce(e)$SAE(FPI,UP(x,FPI),
cp),m)$M2(FPI,V(FPI),V(FPI),M(FPI),SAE(FPI,UP(x,FPI),cp),V(SAE(FPI,UP(x,FPI),
cp)),V(SAE(FPI,UP(x,FPI),cp)),M(SAE(FPI,UP(x,FPI),cp))) - gen * id;

  output(m2::OF);

  makeTempCoeff(deg:I):S == {
    free tmpVarCnt:I;
    local retVal:S := script("TEV"::S,[[(coerce(tmpVarCnt)$I)@OF,(coerce(deg)$I)@OF]$L(OF)]$L(L(OF)));
    tmpVarCnt := tmpVarCnt +$I 1$I;
    retVal
  }

  destructSAE(up:UP(x,FPI),deg:NNI):L(FPI) == {
    import from SEG(NNI);
    coeffList:L(FPI) := []$L(FPI);
    for i in (0$NNI..deg)$SEG(NNI) repeat {
      -- remove coefficient
      coeffList := cons(coefficient(up,i)$UP(x,FPI),coeffList)$L(FPI);
    }
    coeffList;
  }

  --function to generate a generic polynomial from SimpleAlgebraicExtension
  constructSAE(d:NNI):SAE(FPI,UP(x,FPI),cp) == {
    output("start construct");
    import from SEG(I);
    import from FPI;
    output("after imports");
    free genUP:UP(x,FPI);
    output("after free");
    local genericUp:UP(x,FPI) := 0$UP(x,FPI);
    output("after local");
    for i in (0$I..(d pretend I))$SEG(I) repeat {
      output("in for");
      genericUp := genericUp +$(UP(x,FPI)) (coerce(makeTempCoeff(i))$FPI)@FPI *$(UP(x,FPI)) (genUP **$(UP(x,FPI)) (i pretend NNI));
    }
    output("after for");
    reduce(genericUp)$SAE(FPI,UP(x,FPI),cp);
  }

  -- function to perform piecewise gaussian elimination on a matrix, removing
  -- subexpressions at each stage of the elimination procedure.
  myGauss():L(L(L(L(FPI)))) == {

    free len;
    import from SEG(I);				
    import from I;

    output("start myGauss");

    local retVal:L(L(L(L(FPI)))) := []$L(L(L(L(FPI)))); --the final value

    output("before for");
    --loop once for each pass of the elimination
    for i in (1$I..len)$SEG(I) repeat{

      local pass:L(L(L(FPI))); --variable to store each pass
      output("local1");
      -- the first non-zero element of the row is replaced by a symbol
      local value:SAE(FPI,UP(x,FPI),cp);
      output("local2");
      local up1:UP(x,FPI) := lift(apply(m2,i,i)$M(SAE(FPI,UP(x,FPI),cp)))$SAE(FPI,UP(x,FPI),cp);
      output("local3");
      local d1:NNI := degree(up1)$UP(x,FPI);
      output("local4");
      set!(m2,i,i,constructSAE(d1))$M(SAE(FPI,UP(x,FPI),cp));
      output("i",(coerce(i)$I)@OF);
      firstRow:L(L(FPI)):=[]$L(L(FPI));

      --divide the rest of the first row by the first elt
      for j in (len..(i +$I (1$I)))$SEG(I) by (-$I (1$I)) repeat {

        value := apply(m2,i,j)/apply(m2,i,i);
        up := lift(value)$SAE(FPI,UP(x,FPI),cp);
        d:=degree(up)$UP(x,FPI);
        firstRow:=cons(destructSAE(up,d),firstRow);

        output("j1",(coerce(j)$I)@OF);
-- this next line is not necessary
        set!(m2,i,j,constructSAE(d));
        output("m2",m2::OF);
      }

      --now put the first element on the return list
      --N.B. cons is used for efficiency considerations
      firstRow:=cons(destructSAE(up1,d1),firstRow);
      --store in pass variable
      pass := [firstRow]$L(L(L(FPI)));
--This next line is not necessary
      set!(m2,i,i,1$FPI :: SAE(FPI,UP(x,FPI),cp));

      --subtract stuff from the rest of the rows
      for k in ((i +$I (1$I))..len)$SEG(I) repeat {

        --the first non-zero element of each row is replaced by a symbol
        up := lift(apply(m2,k,i))$SAE(FPI,UP(x,FPI),cp);
        d  := degree(up)$UP(x,FPI);
        set!(m2,k,i,constructSAE(d))$M(SAE(FPI,UP(x,FPI),cp));
        local nextRow:L(L(FPI)):=[destructSAE(up,d)]$L(L(FPI));

        --now do the rest of the row
        for j in ((i +$I (1$I))..len)$SEG(I) repeat {

          value := apply(m2,k,j)-apply(m2,k,i)*apply(m2,i,j);
          up := lift(value)$SAE(FPI,UP(x,FPI),cp);
          d:=degree(up)$UP(x,FPI);
          nextRow:=cons(destructSAE(up,d),nextRow);
          output("j2",(coerce(j)$I)@OF);
          set!(m2,k,j,constructSAE(d));
          output("elt(m2,k,j)",apply(m2,k,j)::OF);
        }

--        output("k",k::OF);
        nextRow:=reverse(nextRow)$(L(L(FPI)));
        pass := cons(nextRow,pass)$L(L(L(FPI)));
--This next line is not necessary
        set!(m2,k,i,0$FPI :: SAE(FPI,UP(x,FPI),cp));
      }
      pass := reverse(pass)$L(L(L(FPI)));
    }

    retVal:=cons(pass,retVal)$L(L(L(L(FPI))));
--    output(m2::OF);
--    matrix([[1$FPI]$L(FPI)]$L(L(FPI)));
    reverse(retVal)$L(L(L(L(FPI))));
  }
}