--* From postmaster%watson.vnet.ibm.com@yktvmv.watson.ibm.com  Tue Sep  6 06:03:06 1994
--* Received: from yktvmv-ob.watson.ibm.com by asharp.watson.ibm.com (AIX 3.2/UCB 5.64/930311)
--*           id AA21412; Tue, 6 Sep 1994 06:03:06 -0400
--* Received: from watson.vnet.ibm.com by yktvmv.watson.ibm.com (IBM VM SMTP V2R3)
--*    with BSMTP id 9879; Tue, 06 Sep 94 06:03:10 EDT
--* Received: from YKTVMV by watson.vnet.ibm.com with "VAGENT.V1.0"
--*           id <A.BILL.NOTE.YKTVMV.2507.Sep.06.06:03:09.-0400>
--*           for asbugs@watson; Tue, 06 Sep 94 06:03:10 -0400
--* Received: from ben.britain.eu.net by watson.ibm.com (IBM VM SMTP V2R3) with TCP;
--*    Tue, 06 Sep 94 06:03:01 EDT
--* Received: from iec.co.uk by ben.britain.eu.net via JANET with NIFTP (PP)
--*           id <sg.25481-0@ben.britain.eu.net>; Tue, 6 Sep 1994 11:00:31 +0100
--* From: "W.A.Naylor" <bill@num-alg-grp.co.uk>
--* Date: Tue, 6 Sep 94 10:59:47 BST
--* Message-Id: <7336.9409060959@nags2.nag.co.uk>
--* Received: from co.uk (nags8) by nags2.nag.co.uk (4.1/UK-2.1) id AA07336;
--*           Tue, 6 Sep 94 10:59:47 BST
--* To: asbugs@watson.ibm.com
--* Subject: [2] compilation Bug: libSymeSyme:

--@ Fixed  by:  SSD   Fri Nov 11 10:13:22 EST 1994 
--@ Tested by:  none 
--@ Summary:    Error message improved to more clearly indicate that the error results because eigenCode is out-of-date wrt axextend. 

-- Command line: asharp -O -Fasy -Faso -Flsp -laxiom -Mno-AS_W_WillObsolete -DAxiom -M no-emax aspB.as
-- Version: A# version 0.36.5 for SPARC
-- Original bug file name: aspB.as

--+ terminates compilation with error message
--+
--+ Bug: libSymeSyme:  cannot find 'Symbol' in 'axextend' while reading 'eigenCode'
--+
--+ N.B. eigenCode contains my own asharp code, which does compile.
#include "axiom.as"
#library OtherPackages "roeslib2.asl"

import from OtherPackages;

+++ Author: Bill Naylor
+++ Date Created: Mon Aug  8th
+++ Date Last Updated: Mon Aug  8th
+++ Related Constructors:
+++ \spadtype{AspB} produces Fortran
+++ needed for NAG routine D03PLF

  EIGCODE ==> EigenCode;
  ASPB ==>AspB;
  ROEG ==> RoesSchemeGeneral;

  B    ==> Boolean;
  S    ==> Symbol;
  INT  ==> Integer;
  FLT  ==> Float;
  L    ==> List;
  SW   ==> Switch;
  V2   ==> VectorFunctions2;
  E2   ==> ExpressionFunctions2;
  LS   ==> List Symbol;
  VS   ==> Vector Symbol;
  FT   ==> FortranType;
  FC   ==> FortranCode;
  PI   ==> Polynomial Integer;
  EI   ==> Expression Integer;
  EF   ==> Expression Float;
  EQ   ==> Equation;
  OF   ==> OutputForm;
  SEG  ==> Segment;
  SWT  ==> Union(I: Expression Integer,F: Expression Float,CF: Expression Complex Float,switch: Switch);
  STR  ==> String;
  VEI  ==> Vector Expression Integer;
  VEF  ==> Vector Expression Float;
  MEI  ==> Matrix Expression Integer;
  MEF  ==> Matrix Expression Float;
  NNI  ==> NonNegativeInteger;
  FST  ==> FortranScalarType;
  MF2  ==> MatrixCategoryFunctions2;
  FPI  ==> Fraction Polynomial Integer;
  LPI  ==> List Polynomial Integer;
  LOF  ==> List OutputForm;
  LLOF ==> List List OutputForm;
  LLEI ==> List List Expression Integer;
  VFPI ==> Vector Fraction Polynomial Integer;
  MFPI ==> Matrix Fraction Polynomial Integer;
  SYMTAB  ==> SymbolTable;

  GEVT ==> Record(genVec:Vector Fraction Polynomial Integer,symb:Symbol);
  ETYPEUN ==> Union(S:List List Expression Integer,N:'requireNumeric');
  RETYPE  ==> Record(retEval:ETYPEUN,retEvec:GEVT);

  RRT ==> Record(aBarMatrixField : Matrix Fraction Polynomial Integer,
                 eigenField : RETYPE,
                 alphasField : Vector Expression Integer,
                 numericalFluxField : Vector Expression Integer);

AspB(name:Symbol):with {

    construct : (VFPI,VS,B) -> %;
    coerce : % -> OutputForm;
    outputAsFortran : % -> Void;

}

== add {

    import from EigenCode;
    import from RoesSchemeGeneral;

    import from FT;
    import from SYMTAB;
    import from FortranCode;
    import from Union(fst:FST,void:'void');
    import from S;

    default len:NNI := 1@NNI;
    default lenu:INT := len pretend INT;

    syms:SYMTAB == empty();
    {import {coerce:STR->FST;} from FST;
    default real:FST := ("real"@STR)::FST;}

    declare!(["NPDE"::S,"NCODE"::S]$LS,fortranInteger(),syms);

    vType : FT := construct([real]$Union(fst:FST,void:'void'),["NCODE"::S]$LS,false$B)$FT;

    Rep ==> FortranProgram(name,[void]$Union(fst:FST,void:Enumeration(void)),
["NPDE"::S,"T"::S,"X"::S,"NCODE"::S,"V"::S,"ULEFT"::S,"URIGHT"::S,"FLUX"::S]$LS,syms);

--    import from Rep;

    leftScript(u : S) : S == {
      (concat(string(u)$S,"L")$String)::S;
    }

    rightScript(u : S) : S == {
      (concat(string(u)$S,"R")$String)::S;
    }

    -- my coerce to convert a Fortran Code object into an Expression Float
--    myCoerce(c:FC):EF == {
--      import from S;
--      code(c).callBranch::S::EF;}

    -- mySubst is a substitute function for FPI
    mySubst(f:FPI,le:(L EQ PI)):FPI == {
        (eval(numer(f)$FPI,le)$PI /$FPI eval(denom(f)$FPI,le)$PI);
    }

    -- N.B. This function will construct NAG library calls, which evaluate
    -- eigen values and vectors.
--(possible routines :-
-- F08NHF - balances a real general matrix in order to improve the accuracy
--          of computed eigenvalues /  vectors  ---Perhaps we don't want this
-- F08NEF - reduces a real general matrix to Hessenberg form
-- F08NFF - generates the real orthogonal matrix Q, generated by F08NEF, above
-- F08PEF - computes all the eigenvalues, (optionally the Schur factorization)
--          of a real Hessenberg matrix or a real general matrix which has been
--            reduced to Hessenberg form.
-- F08QKF - computes selected left and/or right eigen vectors of a real upper
--            quasi-triangular matrix
-- F08NJF - transforms eigenvectors of a balanced matrix to those of the
--          original real nonsymmetric matrix)    -- Perhaps we don't want this
-- Alternatively call F02EBF, which calls all of the above

    constrEcall(lwork:INT):L(FC) == {
      import from EI;
      import from S;
      strlen:STR := string(lenu)$STR;
      strlwork:STR := string(lwork)$STR;
      theString:STR := concat(["F02EBF,N,",strlen,",ABAR,",strlen,
",EVALS,EVALSI,EVECS,",strlen,",EVECSI,",strlen,",WORK,",strlwork,
"IFAIL"]$L(STR))$STR;
      theCall:FC := call(theString)$FC;

      [assign("IFAIL"::S,0@EI)$FC,theCall]$L(FC);
    }

    fluxCodeFn(flux:VEF):L(FC) == {

      import from S;
      import from SEG(INT);
      import from OF;
      import from SW;
      import from SWT;

      local fluxNum:EF;
      local fluxDen:EF;
      local fluxCodetmp:L(FC) := nil()$L(FC);
      local fluxQuo:EF := ("FLUXN"::S::EF) /$EF ("FLUXD"::S::EF);
      local elseBlock:FC;
      local thenBlock:FC;
      for i in 1@INT..lenu repeat {
        fluxDen := denominator(flux.i)$EF;
        fluxNum := numerator(flux.i)$EF;
        fluxCodetmp := cons(assign("FLUXN"::S,fluxNum)$FC,fluxCodetmp)$L(FC);
        fluxCodetmp := cons(assign("FLUXD"::S,fluxDen)$FC,fluxCodetmp)$L(FC);
        -- shall only test to see if the denominator is zero
        elseBlock := assign(script("FLUX"::S,[[i::OF]$LOF]$LLOF),fluxQuo)$FC;
        thenBlock := assign(script("FLUX"::S,[[i::OF]$LOF]$LLOF),0@EF)$FC;
        ifStmnt := cond(LT(["FLUXD"::S::EF]$SWT,["LTLNUM"::S::EF]$SWT)$SW_
                          ,thenBlock,elseBlock)$FC;
        fluxCodetmp := cons(ifStmnt,fluxCodetmp)$L(FC);
      }
      fluxCode:L(FC) := reverse(fluxCodetmp)$L(FC);
    }

    aBarCodeFn(aBar:MEF):L(FC) == {

      import from S;
      import from SEG(INT);
      import from OF;
      import from SW;
      import from SWT;

      local eltNum : EF;
      local eltDen : EF;
      local aBarCodetmp:L(FC) := nil()$L(FC);
      local eltQuo:EF := ("ELTN"::S::EF) /$EF ("ELTD"::S::EF);
      local elseBlock:FC;
      local thenBlock:FC;
      for i in 1@INT..lenu repeat {
        for j in 1@INT..lenu repeat {
          eltNum := numerator(apply(aBar,i,j))$EF;
          eltDen := denominator(apply(aBar,i,j))$EF;
          aBarCodetmp := cons(assign("ELTN"::S,eltNum)$FC,aBarCodetmp)$L(FC);
          aBarCodetmp := cons(assign("ELTD"::S,eltDen)$FC,aBarCodetmp)$L(FC);
          -- shall only test to see if the denominator is zero
          elseBlock := assign(script("ABAR"::S,[[i::OF,j::OF]$LOF]$LLOF),eltQuo)$FC;
          thenBlock := assign(script("ABAR"::S,[[i::OF,j::OF]$LOF]$LLOF),0@EF)$FC;
          ifStmnt := cond(LT(["ELTD"::S::EF]$SWT,["LTLNUM"::S::EF]$SWT)$SW_
                          ,thenBlock,elseBlock)$FC;
          aBarCodetmp := cons(ifStmnt,aBarCodetmp)$L(FC);
        }
      }
      eltCode:L(FC) := reverse aBarCodetmp;
    }

    alphaCodeFn(alphaField:VEF):L(FC) == {

      import from EF;
      import from S;
      import from SEG(INT);
      import from OF;
      import from SW;
      import from SWT;

      local alphaNum : EF;
      local alphaDenom : EF;
      local alphaCodetmp:L(FC) := nil()$L(FC);
      local alphaQuo:EF := ("ALPHAN"::S::EF) /$EF ("ALPHAD"::S::EF);
      local elseBlock:FC;
      local innerElseBlock:FC;
      local innerThenBlock:FC;
      local thenBlock:FC;
      local ifStmnt:FC;

      -- for each alpha_i
      for i in 1@INT..lenu repeat {
        -- split alpha_i into numerator / denominator
        alphaNum := numerator(alphaField.i)$EF;
        alphaDenom := denominator(alphaField.i)$EF;
        -- do Fortran evaluation of both
        alphaCodetmp := cons(assign("ALPHAN"::S,alphaNum)$FC,alphaCodetmp)$L(FC);
        alphaCodetmp := cons(assign("ALPHAD"::S,alphaDenom)$FC,alphaCodetmp)$L(FC);
        -- do Fortran if denominator = 0
        -- then       if numerator = 0 then 0 else error
        -- else numerator/denominator
        elseBlock := assign(script("ALPHAS"::S,[[i::OF]$LOF]$LLOF),alphaQuo)$FC;
        innerElseBlock := call("STOP")$FC;
-- ** this is not to cool, as 0/0 not= 0 in general the value will be determined by
--    the limiting value of the alphas around the point
--    an approach might be something like this,
--    calculate all the denominator zero's,
--    then calculate expressions for the limits of the alphas at this point.
--    then in the Fortran, we must determine, which of these points the 0/0 point is
--    and evaluate the respective expression, at that point. (in some cases, AXIOM
--    may be able to give the value, not just an expression) **
        innerThenBlock :=assign(script("ALPHAS"::S,[[i::OF]$LOF]$LLOF),0@EF)$FC;
        thenBlock := cond(LT(["ALPHAN"::S::EF]$SWT,["LTLNUM"::S::EF]$SWT)$SW
                          ,innerThenBlock,innerElseBlock)$FC;
        ifStmnt := cond(LT(["ALPHAD"::S::EF]$SWT,["LTLNUM"::S::EF]$SWT)$SW
                          ,thenBlock,elseBlock)$FC;
        alphaCodetmp := cons(ifStmnt,alphaCodetmp)$L(FC);
      }
      alphaCode:L(FC):=reverse alphaCodetmp;
    }

-- create the symbol tables
    declareSymTab():SYMTAB == {-- note that syms is global, so only return
                              --  locals
      uType : FT := construct([real]$Union(fst:FST,void:'void'),[lenu::PI]$LPI,false$B);
      tmpMatrixType : FT := construct([real]$Union(fst:FST,void:'void'),[lenu::PI,lenu::PI],false$B);

      declare!("T"::S,fortranReal(),syms);
      declare!("X"::S,fortranReal(),syms);
      declare!("V"::S,vType,syms);
      declare!(["ULEFT"::S,"URIGHT"::S,"FLUX"::S]$LS,uType,syms);

      locals : SYMTAB := empty()$SYMTAB;

-- declare the temporary variables
      declare!(["ABAR"::S,"EVECS"::S]$LS,tmpMatrixType,locals);
      declare!(["EVALS"::S,"ALPHAS"::S]$LS,uType,locals);
      declare!(["ABARN"::S,"ABARD"::S,"ALPHAN"::S,"ALPHAD"::S,"FLUXN"::S,"FLUXD"::S]$LS,fortranReal(),locals);
      locals;
    }

    declareSymTab2(locals:SYMTAB):SYMTAB == {

      import from S;

      -- declare IFAIL, EVALSI, EVECSI and WORK
      uType : FT := construct([real]$Union(fst:FST,void:'void'),[lenu::PI]$LPI,false$B);
      declare!(["IFAIL"::S]$LS,fortranReal(),locals);
      declare!(["EVALSI"::S,"EVECSI"::S,"WORK"::S]$LS,uType,locals);
      locals;
    }

-- declare the variables involved in evaluating  eigen values and vectors,
-- in the case that we may do this symbolicaly.
    declareSymTab3(noTempVars:NNI,noCoeffVars:NNI,locals:SYMTAB):SYMTAB == {

      import {coerce:INT -> PI;} from PI;
      import {coerce:STR -> S;} from S;

      noTempVars2:INT := noTempVars pretend INT;
      noCoeffVars2:INT := noCoeffVars pretend INT;

      -- now we know enough to do the local declarations for the temporary
      --  variables
      -- Declare Temporary E. value vars

      tType : FT := construct([real]$Union(fst:FST,void:'void'),
                              [noTempVars2::PI]$LPI,false$B);
      declare!(["T"::S]$LS,tType,locals);

      -- Declare Temporary Coefficient vars
      tcType : FT := construct([real]$Union(fst:FST,void:'void'),
                              [noCoeffVars2::PI]$LPI,false$B);
      declare!(["TC"::S]$LS,tcType,locals);

      -- Declare the Denominator vars
      declare!(["D"::S]$LS,fortranReal(),locals);
      locals;
    }
    map2(x:EI):EF == map(coerce$FLT,x)$E2(INT,FLT);
    map3(x:FPI):EI == x::EI;
    map4(x:EI):EF == map((coerce$FLT)@(INT->FLT),x)$E2(INT,FLT);

-- makeCode is the function which makes the code, that fits the specification
--   for the subroutine NUMFLX
    makeCode(f:VFPI,u:VS,rqrNumEval:B):
               Record(localSymbols:SymbolTable,code:(List FortranCode)) == {

      import from SEG(INT);
      import from GEVT;
      import from RETYPE;
      import from ETYPEUN;
      import from RRT;
      import from PI;
      import from S;

      len := #u;
      lenu := len pretend INT;

      local retVal : RRT := [new(len,len,0$FPI)$MFPI
                      ,[[requireNumeric]$ETYPEUN
                        ,[new(0$NNI,1$FPI)$VFPI,new()$S]$GEVT]$RETYPE
                      ,new(len,0$EI)$VEI
                      ,new(len,0$EI)$VEI]$RRT;

-- create the symbol tables
      local locals : SYMTAB := declareSymTab();

-- form a list of unique variables. (we can't have u+uL, u & uL unrelated
-- variables
      local uniqueVars : VS := vector([new()$S for i in 1@INT..lenu]$LS)$VS;

      local eigcode : L(FC);
-- replace the symbols in the input equations by the new unique ones
      local newExpr : VFPI := new(len,0$FPI);
      for i in 1@INT..lenu repeat {
        newExpr.i := (eval(numer(f.i)$FPI,[(u.j) for j in 1@INT..lenu]$LS,
                     [(uniqueVars.j)::PI for j in 1@INT..lenu]$L(PI))$PI
                   /$FPI eval(denom(f.i)$FPI,[(u.j) for j in 1@INT..lenu]$LS,
                     [(uniqueVars.j)::PI for j in 1@INT..lenu]$L(PI))$PI);
      }

      -- form linearised jacobian
      local roeResults:RRT := roeG(uniqueVars,newExpr,rqrNumEval)$ROEG;
      local abarmat:MFPI := roeResults.aBarMatrixField;

      -- replace the variables uL by corresponding ULEFT_i and URIGHT_i,
      -- in aBarMatrixField, alphasField, and numericalFluxField
      -- N.B. it is not required to replace other vars, as the expressions,
      -- are created in the correct variables.

      local uLeftEqnList : (L EQ PI) := [equation(leftScript(uniqueVars.i)::PI,
                            script("ULEFT"::S,[[i::OF]$LOF]$LLOF)::PI)$(EQ PI)
                            for i in 1@INT..lenu]$L(EQ(PI));
      local uRightEqnList : (L EQ PI) := [equation(rightScript(uniqueVars.i
                            )::PI,script("URIGHT"::S,[[i::OF]$LOF]$LLOF)
                            ::PI)$(EQ PI) for i in 1@INT..lenu]$L(EQ(PI));
      local uLeftEqnListE : (L EQ EI) := [equation(leftScript(uniqueVars.i)::EI
                            ,script("ULEFT"::S,[[i::OF]$LOF]$LLOF)::EI
                            )$(EQ EI) for i in 1@INT..lenu]$L(EQ(EI));

      local uRightEqnListE : (L EQ EI) := [equation(rightScript(uniqueVars.i
                            )::EI,script("URIGHT"::S,[[i::OF]$LOF]$LLOF)
                            ::EI)$(EQ EI) for i in 1@INT..lenu]$L(EQ(EI));

      abarmat := map((x:FPI):FPI +-> mySubst(x,uLeftEqnList),abarmat
                     )$MF2(FPI,VFPI,VFPI,MFPI,FPI,VFPI,VFPI,MFPI);
      abarmat := map((x:FPI):FPI +-> mySubst(x,uRightEqnList),abarmat
                     )$MF2(FPI,VFPI,VFPI,MFPI,FPI,VFPI,VFPI,MFPI);
      local eigcode:L(FC);
      -- do eigen value / vector stuff,
      if retVal.eigenField.retEval case 'requireNumeric'
      then  {-- construct nag routine call to evaluate eigen values / vectors
        local lwork:INT := lenu*(64@INT);
        locals:=declareSymTab2(locals);
        eigcode := constrEcall(lwork);
      }
      else {
        locals := declareSymTab3(#((roeResults.eigenField).retEval.S),
                                 #(roeResults.eigenField.retEvec.genVec),
                                 locals);
        eigcode := makeBlock(roeResults.eigenField)$EIGCODE;

      }

      local alphaFieldEI:VEI := map((x:EI):EI +-> subst(x,uLeftEqnListE)$EI,roeResults.alphasField
                     )$V2(EI,EI);
      alphaFieldEI := map((x:EI):EI +-> subst(x,uRightEqnListE)$EI,alphaFieldEI
                     )$V2(EI,EI);
      local alphaField:VEF := map(map2,alphaFieldEI)$V2(EI,EF);

      local fluxFieldEI:VEI := map((x:EI):EI +-> subst(x,uLeftEqnListE)$EI,
                        roeResults.numericalFluxField)$V2(EI,EI);
      fluxFieldEI := map((x:EI):EI +-> subst(x,uRightEqnListE)$EI,fluxFieldEI)$V2(EI,EI);
      local fluxField:VEF := map(map2,fluxFieldEI)$V2(EI,EF);

      local abarMEI:MEI := map(map3,abarmat)$MF2(FPI,VFPI,VFPI,MFPI,EI,VEI,VEI,MEI);
      abarMEF:MEF := map(map4,abarMEI
                        )$MF2(EI,VEI,VEI,MEI,EF,VEF,VEF,MEF);
      local littleNum:L(FC) := [assign("LTLNUM"::S,"X02AKF()"::S::EF)$FC];
      local aBarCode:L(FC) := aBarCodeFn abarMEF;
      local alphaCode:L(FC) := alphaCodeFn alphaField;
      local fluxCode:L(FC) := fluxCodeFn fluxField;

      local codeList:L(FC) := append(littleNum,append(aBarCode,append(alphaCode,eigcode)));
      [locals,codeList]$Record(localSymbols:SymbolTable,code:L(FC));
    }

    construct(f:VFPI,u:VS,rqrNumeric:B):% == {makeCode(f,u,rqrNumeric)::Rep;}

    coerce(u:%):OF == {coerce(u)$Rep;}

    outputAsFortran(u:%):Void == {outputAsFortran(u)$Rep;}
}