path: root/src/algebra/fortpak.spad.pamphlet

\documentclass{article}
\usepackage{axiom}
\begin{document}
\title{\$SPAD/src/algebra fortpak.spad}
\author{Grant Keady, Godfrey Nolan, Mike Dewar, Themos Tsikas}
\maketitle
\begin{abstract}
\end{abstract}
\eject
\tableofcontents
\eject
\section{package FCPAK1 FortranCodePackage1}
<<package FCPAK1 FortranCodePackage1>>=
)abbrev package FCPAK1 FortranCodePackage1
++ Author: Grant Keady and Godfrey Nolan
++ Date Created: April 1993
++ Date Last Updated: 
++ Basic Operations:
++ Related Constructors:
++ Also See:
++ AMS Classifications:
++ Keywords:
++ References:
++ Description:
++  \spadtype{FortranCodePackage1} provides some utilities for
++  producing useful objects in FortranCode domain.
++  The Package may be used with the FortranCode domain and its
++  \spad{printCode} or possibly via an outputAsFortran.
++  (The package provides items of use in connection with ASPs
++  in the AXIOM-NAG link and, where appropriate, naming accords
++  with that in IRENA.)
++  The easy-to-use functions use Fortran loop variables I1, I2,
++  and it is users' responsibility to check that this is sensible.
++  The advanced functions use SegmentBinding to allow users control
++  over Fortran loop variable names.
-- Later might add functions to build
-- diagonalMatrix from List, i.e. the FC version of the corresponding
-- AXIOM function from MatrixCategory;
-- bandedMatrix, i.e. the full-matrix-FC version of the corresponding
-- AXIOM function in BandedMatrix Domain
-- bandedSymmetricMatrix, i.e. the full-matrix-FC version of the corresponding
-- AXIOM function in BandedSymmetricMatrix Domain

FortranCodePackage1: Exports  == Implementation where

  NNI    ==> NonNegativeInteger
  PI     ==> PositiveInteger
  PIN    ==> Polynomial(Integer)
  SBINT  ==> SegmentBinding(Integer)
  SEGINT ==> Segment(Integer)
  LSBINT ==> List(SegmentBinding(Integer))
  SBPIN  ==> SegmentBinding(Polynomial(Integer))
  SEGPIN ==> Segment(Polynomial(Integer))
  LSBPIN ==> List(SegmentBinding(Polynomial(Integer)))
  FC     ==> FortranCode
  EXPRESSION  ==> Union(Expression Integer,Expression Float,Expression Complex Integer,Expression Complex Float)

  Exports == with

    zeroVector: (Symbol,PIN) -> FC
      ++ zeroVector(s,p) \undocumented{}

    zeroMatrix: (Symbol,PIN,PIN) -> FC
      ++ zeroMatrix(s,p,q) uses loop variables in the Fortran, I1 and I2

    zeroMatrix: (Symbol,SBPIN,SBPIN) -> FC
      ++ zeroMatrix(s,b,d) in this version gives the user control 
      ++ over names of Fortran variables used in loops.

    zeroSquareMatrix: (Symbol,PIN) -> FC
      ++ zeroSquareMatrix(s,p) \undocumented{}

    identitySquareMatrix: (Symbol,PIN) -> FC
      ++ identitySquareMatrix(s,p) \undocumented{}

  Implementation ==> add
    import FC

    zeroVector(fname:Symbol,n:PIN):FC ==
      ue:Expression(Integer) := 0
      i1:Symbol := "I1"::Symbol
      lp1:PIN := 1::PIN
      hp1:PIN := n
      segp1:SEGPIN:= segment(lp1,hp1)$SEGPIN
      segbp1:SBPIN := equation(i1,segp1)$SBPIN
      ip1:PIN := i1::PIN
      indices:List(PIN) := [ip1]
      fa:FC := forLoop(segbp1,assign(fname,indices,ue)$FC)$FC
      fa

    zeroMatrix(fname:Symbol,m:PIN,n:PIN):FC ==
      ue:Expression(Integer) := 0
      i1:Symbol := "I1"::Symbol
      lp1:PIN := 1::PIN
      hp1:PIN := m
      segp1:SEGPIN:= segment(lp1,hp1)$SEGPIN
      segbp1:SBPIN := equation(i1,segp1)$SBPIN
      i2:Symbol := "I2"::Symbol
      hp2:PIN := n
      segp2:SEGPIN:= segment(lp1,hp2)$SEGPIN
      segbp2:SBPIN := equation(i2,segp2)$SBPIN
      ip1:PIN := i1::PIN
      ip2:PIN := i2::PIN
      indices:List(PIN) := [ip1,ip2]
      fa:FC :=forLoop(segbp1,forLoop(segbp2,assign(fname,indices,ue)$FC)$FC)$FC
      fa

    zeroMatrix(fname:Symbol,segbp1:SBPIN,segbp2:SBPIN):FC ==
      ue:Expression(Integer) := 0
      i1:Symbol := variable(segbp1)$SBPIN
      i2:Symbol := variable(segbp2)$SBPIN
      ip1:PIN := i1::PIN
      ip2:PIN := i2::PIN
      indices:List(PIN) := [ip1,ip2]
      fa:FC :=forLoop(segbp1,forLoop(segbp2,assign(fname,indices,ue)$FC)$FC)$FC
      fa

    zeroSquareMatrix(fname:Symbol,n:PIN):FC ==
      ue:Expression(Integer) := 0
      i1:Symbol := "I1"::Symbol
      lp1:PIN := 1::PIN
      hp1:PIN := n
      segp1:SEGPIN:= segment(lp1,hp1)$SEGPIN
      segbp1:SBPIN := equation(i1,segp1)$SBPIN
      i2:Symbol := "I2"::Symbol
      segbp2:SBPIN := equation(i2,segp1)$SBPIN
      ip1:PIN := i1::PIN
      ip2:PIN := i2::PIN
      indices:List(PIN) := [ip1,ip2]
      fa:FC :=forLoop(segbp1,forLoop(segbp2,assign(fname,indices,ue)$FC)$FC)$FC
      fa

    identitySquareMatrix(fname:Symbol,n:PIN):FC ==
      ue:Expression(Integer) := 0
      u1:Expression(Integer) := 1
      i1:Symbol := "I1"::Symbol
      lp1:PIN := 1::PIN
      hp1:PIN := n
      segp1:SEGPIN:= segment(lp1,hp1)$SEGPIN
      segbp1:SBPIN := equation(i1,segp1)$SBPIN
      i2:Symbol := "I2"::Symbol
      segbp2:SBPIN := equation(i2,segp1)$SBPIN
      ip1:PIN := i1::PIN
      ip2:PIN := i2::PIN
      indice1:List(PIN) := [ip1,ip1]
      indices:List(PIN) := [ip1,ip2]
      fc:FC := forLoop(segbp2,assign(fname,indices,ue)$FC)$FC
      f1:FC := assign(fname,indice1,u1)$FC
      fl:List(FC) := [fc,f1]
      fa:FC := forLoop(segbp1,block(fl)$FC)$FC
      fa

@
\section{package NAGSP NAGLinkSupportPackage}
<<package NAGSP NAGLinkSupportPackage>>=
)abbrev package NAGSP NAGLinkSupportPackage
++ Author: Mike Dewar and Godfrey Nolan
++ Date Created:  March 1993
++ Date Last Updated: March 4 1994
++                    October 6 1994
++ Basic Operations:
++ Related Domains:
++ Also See:
++ AMS Classifications:
++ Keywords:
++ Examples:
++ References:
++ Description: Support functions for the NAG Library Link functions
NAGLinkSupportPackage() : exports == implementation where

  exports ==> with
    fortranCompilerName : () -> String
      ++ fortranCompilerName() returns the name of the currently selected
      ++  Fortran compiler
    fortranLinkerArgs   : () -> String
      ++ fortranLinkerArgs() returns the current linker arguments
    aspFilename         : String -> String
      ++ aspFilename("f") returns a String consisting of "f" suffixed with
      ++  an extension identifying the current AXIOM session.
    dimensionsOf	: (Symbol, Matrix DoubleFloat) -> SExpression
	++ dimensionsOf(s,m) \undocumented{}
    dimensionsOf	: (Symbol, Matrix Integer) -> SExpression
	++ dimensionsOf(s,m) \undocumented{}
    checkPrecision      : () -> Boolean
	++ checkPrecision() \undocumented{}
    restorePrecision    : () -> Void
	++ restorePrecision() \undocumented{}

  implementation ==> add
    makeAs:                   (Symbol,Symbol) -> Symbol
    changeVariables:          (Expression Integer,Symbol) -> Expression Integer
    changeVariablesF:         (Expression Float,Symbol) -> Expression Float

    import String
    import Symbol

    checkPrecision():Boolean ==
      (_$fortranPrecision$Lisp = "single"::Symbol) and (_$nagEnforceDouble$Lisp)  =>
        systemCommand("set fortran precision double")$MoreSystemCommands
        if _$nagMessages$Lisp  then 
          print("*** Warning: Resetting fortran precision to double")$PrintPackage
        true
      false

    restorePrecision():Void ==
      systemCommand("set fortran precision single")$MoreSystemCommands
      if _$nagMessages$Lisp  then 
        print("** Warning: Restoring fortran precision to single")$PrintPackage
      void()$Void

    uniqueId : String := ""
    counter : Integer := 0
    getUniqueId():String ==
      if uniqueId = "" then
        uniqueId := concat(getEnv("HOST")$Lisp,getEnv("SPADNUM")$Lisp)
      concat(uniqueId,string (counter:=counter+1))

    fortranCompilerName() == string _$fortranCompilerName$Lisp
    fortranLinkerArgs() == string _$fortranLibraries$Lisp

    aspFilename(f:String):String == concat ["/tmp/",f,getUniqueId(),".f"]

    dimensionsOf(u:Symbol,m:Matrix DoubleFloat):SExpression ==
      [u,nrows m,ncols m]$Lisp
    dimensionsOf(u:Symbol,m:Matrix Integer):SExpression ==
      [u,nrows m,ncols m]$Lisp

@
\section{package FORT FortranPackage}
<<package FORT FortranPackage>>=
)abbrev package FORT FortranPackage
-- Because of a bug in the compiler:
)bo $noSubsumption:=true 

++ Author: Mike Dewar
++ Date Created: October 6 1991
++ Date Last Updated: 13 July 1994
++ Basic Operations: linkToFortran
++ Related Constructors:
++ Also See:
++ AMS Classifications:
++ Keywords: 
++ References:
++ Description: provides an interface to the boot code for calling Fortran
FortranPackage(): Exports == Implementation where
 FST ==> FortranScalarType
 SEX ==> SExpression
 L   ==> List
 S   ==> Symbol
 FOP ==> FortranOutputStackPackage
 U   ==> Union(array:L S,scalar:S)

 Exports ==> with
  linkToFortran: (S, L U, L L U, L S) -> SEX
	++ linkToFortran(s,l,ll,lv) \undocumented{}
  linkToFortran: (S, L U, L L U, L S, S) -> SEX
	++ linkToFortran(s,l,ll,lv,t) \undocumented{}
  linkToFortran: (S,L S,TheSymbolTable,L S) -> SEX
	++ linkToFortran(s,l,t,lv) \undocumented{}
  outputAsFortran: FileName -> Void
	++ outputAsFortran(fn) \undocumented{}
  setLegalFortranSourceExtensions: List String -> List String
	++ setLegalFortranSourceExtensions(l) \undocumented{}

 Implementation ==> add

  legalFortranSourceExtensions : List String := ["f"]

  setLegalFortranSourceExtensions(l:List String):List String ==
    legalFortranSourceExtensions := l
    
  checkExtension(fn : FileName) : String ==
    -- Does it end in a legal extension ?
    stringFn := fn::String
    not member?(extension fn,legalFortranSourceExtensions) =>
      error [stringFn,"is not a legal Fortran Source File."]
    stringFn

  outputAsFortran(fn:FileName):Void ==
--    source : String := checkExtension fn
    source : String := fn::String
    not readable? fn => 
      popFortranOutputStack()$FOP
      error([source,"is not readable"]@List(String))
    target : String := topFortranOutputStack()$FOP
    command : String := 
      concat(["sys rm -f ",target," ; cp ",source," ",target])$String
    systemCommand(command)$MoreSystemCommands
    void()$Void

  linkToFortran(name:S,args:L U, decls:L L U, res:L(S)):SEX == 
    makeFort(name,args,decls,res,NIL$Lisp,NIL$Lisp)$Lisp

  linkToFortran(name:S,args:L U, decls:L L U, res:L(S),returnType:S):SEX == 
    makeFort(name,args,decls,res,returnType,NIL$Lisp)$Lisp

  dimensions(type:FortranType):SEX ==
    convert([convert(convert(u)@InputForm)@SEX _
      for u in dimensionsOf(type)])@SEX

  ftype(name:S,type:FortranType):SEX ==
    [name,scalarTypeOf(type),dimensions(type),external? type]$Lisp

  makeAspList(asp:S,syms:TheSymbolTable):SExpression==
    symtab : SymbolTable := symbolTableOf(asp,syms)
    [asp,returnTypeOf(asp,syms),argumentListOf(asp,syms), _
     [ftype(u,fortranTypeOf(u,symtab)) for u in parametersOf symtab]]$Lisp

  linkToFortran(name:S,aArgs:L S,syms:TheSymbolTable,res:L S):SEX ==
    arguments : L S := argumentListOf(name,syms)$TheSymbolTable
    dummies : L S := setDifference(arguments,aArgs)
    symbolTable:SymbolTable := symbolTableOf(name,syms)
    symbolList := newTypeLists(symbolTable)
    rt:Union(fst: FST,void: "void") := returnTypeOf(name,syms)$TheSymbolTable

    -- Look for arguments which are subprograms
    asps :=[makeAspList(u,syms) for u in externalList(symbolTable)$SymbolTable]
    rt case fst =>
      makeFort1(name,arguments,aArgs,dummies,symbolList,res,(rt.fst)::S,asps)$Lisp
    makeFort1(name,arguments,aArgs,dummies,symbolList,res,NIL$Lisp,asps)$Lisp

@
\section{package FOP FortranOutputStackPackage}
<<package FOP FortranOutputStackPackage>>=
)abbrev package FOP FortranOutputStackPackage
-- Because of a bug in the compiler:
)bo $noSubsumption:=false

++ Author: Mike Dewar
++ Date Created:  October 1992
++ Date Last Updated: 
++ Basic Operations:
++ Related Domains:
++ Also See:
++ AMS Classifications:
++ Keywords:
++ Examples:
++ References:
++ Description: Code to manipulate Fortran Output Stack
FortranOutputStackPackage() : specification == implementation where

  specification == with

    clearFortranOutputStack : () -> Stack String
      ++ clearFortranOutputStack() clears the Fortran output stack
    showFortranOutputStack : () -> Stack String
      ++ showFortranOutputStack() returns the Fortran output stack
    popFortranOutputStack : () -> Void
      ++ popFortranOutputStack() pops the Fortran output stack
    pushFortranOutputStack : FileName -> Void
      ++ pushFortranOutputStack(f) pushes f onto the Fortran output stack
    pushFortranOutputStack : String -> Void
      ++ pushFortranOutputStack(f) pushes f onto the Fortran output stack
    topFortranOutputStack : () -> String
      ++ topFortranOutputStack() returns the top element of the Fortran
      ++ output stack

  implementation == add

    import MoreSystemCommands

    -- A stack of filenames for Fortran output.  We are sharing this with
    -- the standard Fortran output code, so want to be a bit careful about
    -- how we interact with what the user does independently.  We get round
    -- potential problems by always examining the top element of the stack 
    -- before we push.  If the user has redirected output then we alter our
    -- top value accordingly.
    fortranOutputStack : Stack String := empty()@(Stack String)

    topFortranOutputStack():String == string(_$fortranOutputFile$Lisp)

    pushFortranOutputStack(fn:FileName):Void ==
      if empty? fortranOutputStack then
        push!(string(_$fortranOutputFile$Lisp),fortranOutputStack)
      else if not(top(fortranOutputStack)=string(_$fortranOutputFile$Lisp)) then
        pop! fortranOutputStack
        push!(string(_$fortranOutputFile$Lisp),fortranOutputStack)
      push!( fn::String,fortranOutputStack)
      systemCommand concat(["set output fortran quiet ", fn::String])$String
      void()

    pushFortranOutputStack(fn:String):Void ==
      if empty? fortranOutputStack then
        push!(string(_$fortranOutputFile$Lisp),fortranOutputStack)
      else if not(top(fortranOutputStack)=string(_$fortranOutputFile$Lisp)) then
        pop! fortranOutputStack
        push!(string(_$fortranOutputFile$Lisp),fortranOutputStack)
      push!( fn,fortranOutputStack)
      systemCommand concat(["set output fortran quiet ", fn])$String
      void()

    popFortranOutputStack():Void ==
      if not empty? fortranOutputStack then pop! fortranOutputStack
      if empty? fortranOutputStack then push!("CONSOLE",fortranOutputStack)
      systemCommand concat(["set output fortran quiet append ",_
                           top fortranOutputStack])$String
      void()

    clearFortranOutputStack():Stack String ==
      fortranOutputStack := empty()@(Stack String)

    showFortranOutputStack():Stack String ==
      fortranOutputStack

@
\section{package TEMUTL TemplateUtilities}
<<package TEMUTL TemplateUtilities>>=
)abbrev package TEMUTL TemplateUtilities
++ Author: Mike Dewar
++ Date Created:  October 1992
++ Date Last Updated: 
++ Basic Operations:
++ Related Domains:
++ Also See:
++ AMS Classifications:
++ Keywords:
++ Examples:
++ References:
++ Description: This package provides functions for template manipulation
TemplateUtilities(): Exports == Implementation where

  Exports == with
    interpretString : String -> Any
      ++ interpretString(s) treats a string as a piece of AXIOM input, by
      ++ parsing and interpreting it.
    stripCommentsAndBlanks : String -> String
      ++ stripCommentsAndBlanks(s) treats s as a piece of AXIOM input, and
      ++ removes comments, and leading and trailing blanks.

  Implementation == add

    import InputForm

    stripC(s:String,u:String):String ==
      i : Integer := position(u,s,1)
      i = 0 => s
      delete(s,i..)

    stripCommentsAndBlanks(s:String):String ==
      trim(stripC(stripC(s,"++"),"--"),char " ")

    parse(s:String):InputForm ==
      ncParseFromString(s)$Lisp::InputForm 

    interpretString(s:String):Any ==
      interpret parse s

@
\section{package MCALCFN MultiVariableCalculusFunctions}
<<package MCALCFN MultiVariableCalculusFunctions>>=
)abbrev package MCALCFN MultiVariableCalculusFunctions
++ Author: Themos Tsikas, Grant Keady
++ Date Created: December 1992
++ Date Last Updated: June 1993
++ Basic Operations:
++ Related Constructors:
++ Also See:
++ AMS Classifications:
++ Keywords:
++ References:
++ Description:
++  \spadtype{MultiVariableCalculusFunctions} Package provides several
++  functions for multivariable calculus.
++ These include gradient, hessian and jacobian,
++ divergence and laplacian.
++ Various forms for banded and sparse storage of matrices are
++ included.
MultiVariableCalculusFunctions(S,F,FLAF,FLAS) : Exports == Implementation where
  PI ==> PositiveInteger
  NNI ==> NonNegativeInteger

  S: SetCategory
  F: PartialDifferentialRing(S)
  FLAS: FiniteLinearAggregate(S)
    with finiteAggregate
  FLAF: FiniteLinearAggregate(F)

  Exports ==> with
    gradient: (F,FLAS) -> Vector F
     ++ \spad{gradient(v,xlist)}
     ++ computes the gradient, the vector of first partial derivatives,
     ++ of the scalar field v,
     ++ v a function of the variables listed in xlist.
    divergence: (FLAF,FLAS) ->  F
     ++ \spad{divergence(vf,xlist)}
     ++ computes the divergence of the vector field vf,
     ++ vf a vector function of the variables listed in xlist.
    laplacian: (F,FLAS) -> F
     ++ \spad{laplacian(v,xlist)}
     ++ computes the laplacian of the scalar field v,
     ++ v a function of the variables listed in xlist.
    hessian: (F,FLAS) -> Matrix F
     ++ \spad{hessian(v,xlist)}
     ++ computes the hessian, the matrix of second partial derivatives,
     ++ of the scalar field v,
     ++ v a function of the variables listed in xlist.
    bandedHessian: (F,FLAS,NNI) -> Matrix F
     ++ \spad{bandedHessian(v,xlist,k)}
     ++ computes the hessian, the matrix of second partial derivatives,
     ++ of the scalar field v,
     ++ v a function of the variables listed in xlist,
     ++ k is the semi-bandwidth, the number of nonzero subdiagonals,
     ++ 2*k+1 being actual bandwidth.
     ++ Stores the nonzero band in lower triangle in a matrix, 
     ++ dimensions k+1 by #xlist,
     ++ whose rows are the vectors formed by diagonal, subdiagonal, etc.
     ++ of the real, full-matrix, hessian.
     ++ (The notation conforms to LAPACK/NAG-F07 conventions.)
    -- At one stage it seemed a good idea to help the ASP<n> domains
    -- with the types of their input arguments and this led to the
    -- standard Gradient|Hessian|Jacobian functions.
    --standardJacobian: (Vector(F),List(S)) -> Matrix F
    -- ++ \spad{jacobian(vf,xlist)}
    -- ++ computes the jacobian, the matrix of first partial derivatives,
    -- ++ of the vector field vf,
    -- ++ vf a vector function of the variables listed in xlist.
    jacobian: (FLAF,FLAS) -> Matrix F
     ++ \spad{jacobian(vf,xlist)}
     ++ computes the jacobian, the matrix of first partial derivatives,
     ++ of the vector field vf,
     ++ vf a vector function of the variables listed in xlist.
    bandedJacobian: (FLAF,FLAS,NNI,NNI) -> Matrix F
     ++ \spad{bandedJacobian(vf,xlist,kl,ku)}
     ++ computes the jacobian, the matrix of first partial derivatives,
     ++ of the vector field vf,
     ++ vf a vector function of the variables listed in xlist,
     ++ kl is the number of nonzero subdiagonals,
     ++ ku is the number of nonzero superdiagonals,
     ++ kl+ku+1 being actual bandwidth.
     ++ Stores the nonzero band in a matrix, 
     ++ dimensions kl+ku+1 by #xlist.
     ++ The upper triangle is in the top ku rows,
     ++ the diagonal is in row ku+1,
     ++ the lower triangle in the last kl rows.
     ++ Entries in a column in the band store correspond to entries
     ++ in same column of full store.
     ++ (The notation conforms to LAPACK/NAG-F07 conventions.)

  Implementation ==> add
    localGradient(v:F,xlist:List(S)):Vector(F) ==
       vector([D(v,x) for x in xlist])
    gradient(v,xflas) ==
       --xlist:List(S) := [xflas(i) for i in 1 .. maxIndex(xflas)]
       xlist:List(S) := parts(xflas)
       localGradient(v,xlist)
    localDivergence(vf:Vector(F),xlist:List(S)):F ==
       i: PI
       n: NNI
       ans: F
       -- Perhaps should report error if two args of min different
       n:= min(#(xlist),((maxIndex(vf))::NNI))$NNI
       ans:= 0
       for i in 1 .. n repeat ans := ans + D(vf(i),xlist(i)) 
       ans
    divergence(vf,xflas) ==
       xlist:List(S) := parts(xflas)
       i: PI
       n: NNI
       ans: F
       -- Perhaps should report error if two args of min different
       n:= min(#(xlist),((maxIndex(vf))::NNI))$NNI
       ans:= 0
       for i in 1 .. n repeat ans := ans + D(vf(i),xlist(i)) 
       ans
    laplacian(v,xflas) ==
       xlist:List(S) := parts(xflas)
       gv:Vector(F) := localGradient(v,xlist)
       localDivergence(gv,xlist)
    hessian(v,xflas) ==
       xlist:List(S) := parts(xflas)
       matrix([[D(v,[x,y]) for x in xlist] for y in xlist])
    --standardJacobian(vf,xlist) ==
    --   i: PI
    --   matrix([[D(vf(i),x) for x in xlist] for i in 1 .. maxIndex(vf)])
    jacobian(vf,xflas) ==
       xlist:List(S) := parts(xflas)
       i: PI
       matrix([[D(vf(i),x) for x in xlist] for i in 1 .. maxIndex(vf)])
    bandedHessian(v,xflas,k) ==
       xlist:List(S) := parts(xflas)
       j,iw: PI
       n: NNI
       bandM: Matrix F
       n:= #(xlist)
       bandM:= new(k+1,n,0)
       for j in 1 .. n repeat setelt(bandM,1,j,D(v,xlist(j),2))
       for iw in 2 .. (k+1) repeat (_
         for j in 1 .. (n-iw+1) repeat (_
           setelt(bandM,iw,j,D(v,[xlist(j),xlist(j+iw-1)])) ) )
       bandM
    bandedJacobian(vf,xflas,kl,ku) ==
       xlist:List(S) := parts(xflas)
       j,iw: PI
       n: NNI
       bandM: Matrix F
       n:= #(xlist)
       bandM:= new(kl+ku+1,n,0)
       for j in 1 .. n repeat setelt(bandM,ku+1,j,D(vf(j),xlist(j)))
       for iw in (ku+2) .. (ku+kl+1) repeat (_
         for j in 1 .. (n-iw+ku+1) repeat (_
           setelt(bandM,iw,j,D(vf(j+iw-1-ku),xlist(j))) ) )
       for iw in 1 .. ku repeat (_
         for j in (ku+2-iw) .. n repeat (_
           setelt(bandM,iw,j,D(vf(j+iw-1-ku),xlist(j))) ) )
       bandM

@
\section{License}
<<license>>=
--Copyright (c) 1991-2002, The Numerical ALgorithms Group Ltd.
--All rights reserved.
--
--Redistribution and use in source and binary forms, with or without
--modification, are permitted provided that the following conditions are
--met:
--
--    - Redistributions of source code must retain the above copyright
--      notice, this list of conditions and the following disclaimer.
--
--    - Redistributions in binary form must reproduce the above copyright
--      notice, this list of conditions and the following disclaimer in
--      the documentation and/or other materials provided with the
--      distribution.
--
--    - Neither the name of The Numerical ALgorithms Group Ltd. nor the
--      names of its contributors may be used to endorse or promote products
--      derived from this software without specific prior written permission.
--
--THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS
--IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED
--TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A
--PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER
--OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
--EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
--PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
--PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
--LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
--NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
--SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
@
<<*>>=
<<license>>

<<package FCPAK1 FortranCodePackage1>>
<<package NAGSP NAGLinkSupportPackage>>
<<package FORT FortranPackage>>
<<package FOP FortranOutputStackPackage>>
<<package TEMUTL TemplateUtilities>>
<<package MCALCFN MultiVariableCalculusFunctions>>
@
\eject
\begin{thebibliography}{99}
\bibitem{1} nothing
\end{thebibliography}
\end{document}