axiom-source 20170501-3 / usr / share / axiom-20170501 / src / algebra / ASP31.spad

)abbrev domain ASP31 Asp31
++ Author: Mike Dewar, Grant Keady and Godfrey Nolan
++ Date Created: Mar 1993
++ Date Last Updated: 6 October 1994
++ References:
++ Hawk95 Two more links to NAG numerics involving CA systems
++ Kead93 Production of Argument SubPrograms in the AXIOM -- NAG link
++ Description:
++\spadtype{Asp31} produces Fortran for Type 31 ASPs, needed for NAG routine 
++d02ejf, for example:
++
++\tab{5}SUBROUTINE PEDERV(X,Y,PW)\br
++\tab{5}DOUBLE PRECISION X,Y(*)\br
++\tab{5}DOUBLE PRECISION PW(3,3)\br
++\tab{5}PW(1,1)=-0.03999999999999999D0\br
++\tab{5}PW(1,2)=10000.0D0*Y(3)\br
++\tab{5}PW(1,3)=10000.0D0*Y(2)\br
++\tab{5}PW(2,1)=0.03999999999999999D0\br
++\tab{5}PW(2,2)=(-10000.0D0*Y(3))+(-60000000.0D0*Y(2))\br
++\tab{5}PW(2,3)=-10000.0D0*Y(2)\br
++\tab{5}PW(3,1)=0.0D0\br
++\tab{5}PW(3,2)=60000000.0D0*Y(2)\br
++\tab{5}PW(3,3)=0.0D0\br
++\tab{5}RETURN\br
++\tab{5}END

Asp31(name) : SIG == CODE where
  name : Symbol

  O      ==> OutputForm
  FST    ==> FortranScalarType
  UFST   ==> Union(fst:FST,void:"void")
  MFLOAT ==> MachineFloat
  FEXPR  ==> FortranExpression(['X],['Y],MFLOAT)
  FT     ==> FortranType
  FC     ==> FortranCode
  SYMTAB ==> SymbolTable
  RSFC   ==> Record(localSymbols:SymbolTable,code:List(FortranCode))
  FRAC   ==> Fraction
  POLY   ==> Polynomial
  EXPR   ==> Expression
  INT    ==> Integer
  FLOAT  ==> Float
  VEC    ==> Vector
  MAT    ==> Matrix
  VF2    ==> VectorFunctions2
  MF2    ==> MatrixCategoryFunctions2(FEXPR,VEC FEXPR,VEC FEXPR,MAT FEXPR,
                   EXPR MFLOAT,VEC EXPR MFLOAT,VEC EXPR MFLOAT,MAT EXPR MFLOAT)

  SIG ==> FortranVectorFunctionCategory with

    coerce : VEC FEXPR -> $
      ++coerce(f) takes objects from the appropriate instantiation of
      ++\spadtype{FortranExpression} and turns them into an ASP.

  CODE ==> add

    real : UFST := ["real"::FST]$UFST

    syms : SYMTAB := empty()

    declare!(X,fortranReal(),syms)$SYMTAB

    yType : FT := construct(real,["*"::Symbol],false)$FT

    declare!(Y,yType,syms)$SYMTAB

    Rep := FortranProgram(name,["void"]$UFST,[X,Y,PW],syms)

    -- To help the poor old compiler!
    fexpr2expr(u:FEXPR):EXPR MFLOAT == coerce(u)$FEXPR

    localAssign(s:Symbol,j:MAT FEXPR):FC ==
      j' : MAT EXPR MFLOAT := map(fexpr2expr,j)$MF2
      assign(s,j')$FC

    makeXList(n:Integer):List(Symbol) ==
      j:Integer
      y:Symbol := Y::Symbol
      p:List(Symbol) := []
      for j in 1 .. n repeat p:= cons(subscript(y,[j::OutputForm])$Symbol,p)
      p:= reverse(p)

    coerce(u:VEC FEXPR):$ == 
      dimension := #u::Polynomial Integer
      locals : SYMTAB := empty()
      declare!(PW,[real,[dimension,dimension],false]$FT,locals)$SYMTAB
      n:Integer := maxIndex(u)$VEC(FEXPR)
      p:List(Symbol) := makeXList(n)
      jac: MAT FEXPR := jacobian(u,p)$MultiVariableCalculusFunctions(_
                                     Symbol,FEXPR ,VEC FEXPR,List(Symbol))
      code : List FC := [localAssign(PW,jac),returns()$FC]$List(FC)
      ([locals,code]$RSFC)::$

    retract(u:VEC FRAC POLY INT):$ ==
      v : VEC FEXPR := map(retract,u)$VF2(FRAC POLY INT,FEXPR)
      v::$

    retractIfCan(u:VEC FRAC POLY INT):Union($,"failed") ==
      v:Union(VEC FEXPR,"failed"):=map(retractIfCan,u)$VF2(FRAC POLY INT,FEXPR)
      v case "failed" => "failed"
      (v::VEC FEXPR)::$

    retract(u:VEC FRAC POLY FLOAT):$ ==
      v : VEC FEXPR := map(retract,u)$VF2(FRAC POLY FLOAT,FEXPR)
      v::$

    retractIfCan(u:VEC FRAC POLY FLOAT):Union($,"failed") ==
      v:Union(VEC FEXPR,"failed"):=_
        map(retractIfCan,u)$VF2(FRAC POLY FLOAT,FEXPR)
      v case "failed" => "failed"
      (v::VEC FEXPR)::$

    retract(u:VEC EXPR INT):$ ==
      v : VEC FEXPR := map(retract,u)$VF2(EXPR INT,FEXPR)
      v::$

    retractIfCan(u:VEC EXPR INT):Union($,"failed") ==
      v:Union(VEC FEXPR,"failed"):=map(retractIfCan,u)$VF2(EXPR INT,FEXPR)
      v case "failed" => "failed"
      (v::VEC FEXPR)::$

    retract(u:VEC EXPR FLOAT):$ ==
      v : VEC FEXPR := map(retract,u)$VF2(EXPR FLOAT,FEXPR)
      v::$

    retractIfCan(u:VEC EXPR FLOAT):Union($,"failed") ==
      v:Union(VEC FEXPR,"failed"):=map(retractIfCan,u)$VF2(EXPR FLOAT,FEXPR)
      v case "failed" => "failed"
      (v::VEC FEXPR)::$

    retract(u:VEC POLY INT):$ ==
      v : VEC FEXPR := map(retract,u)$VF2(POLY INT,FEXPR)
      v::$

    retractIfCan(u:VEC POLY INT):Union($,"failed") ==
      v:Union(VEC FEXPR,"failed"):=map(retractIfCan,u)$VF2(POLY INT,FEXPR)
      v case "failed" => "failed"
      (v::VEC FEXPR)::$

    retract(u:VEC POLY FLOAT):$ ==
      v : VEC FEXPR := map(retract,u)$VF2(POLY FLOAT,FEXPR)
      v::$

    retractIfCan(u:VEC POLY FLOAT):Union($,"failed") ==
      v:Union(VEC FEXPR,"failed"):=map(retractIfCan,u)$VF2(POLY FLOAT,FEXPR)
      v case "failed" => "failed"
      (v::VEC FEXPR)::$

    coerce(c:List FC):$ == coerce(c)$Rep

    coerce(r:RSFC):$ == coerce(r)$Rep

    coerce(c:FC):$ == coerce(c)$Rep

    coerce(u:$):O == coerce(u)$Rep

    outputAsFortran(u):Void ==
      p := checkPrecision()$NAGLinkSupportPackage
      outputAsFortran(u)$Rep
      p => restorePrecision()$NAGLinkSupportPackage