diff options
Diffstat (limited to 'src/interp/fortcall.boot')
| -rw-r--r-- | src/interp/fortcall.boot | 798 |
1 files changed, 798 insertions, 0 deletions
diff --git a/src/interp/fortcall.boot b/src/interp/fortcall.boot new file mode 100644 index 00000000..72d79948 --- /dev/null +++ b/src/interp/fortcall.boot @@ -0,0 +1,798 @@ +-- 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. + + +makeFort(name,args,decls,results,returnType,aspInfo) == + -- Create an executable Fortran file to call a given library function, + -- and a stub Axiom function to process its arguments. + -- the following is a list of objects for which values need not be + -- passed by the user. + dummies := [SECOND(u) for u in args | EQUAL(car u,0)] + args := [untangle2(u) for u in args] -- lose spad Union representation + where untangle2 u == + atom (v := rest(u)) => v + first(v) + userArgs := [u for u in args | not member(u,dummies)] -- Temporary + decls := [untangle(u) for u in decls] -- lose spad Union representation + where untangle u == + [if atom(rest(v)) then rest(v) else _ + [if atom(w) then w else rest(w) for w in rest(v)] for v in u] + makeFort1(name,args,userArgs,dummies,decls,results,returnType,aspInfo) + +makeFort1(name,args,userArgs,dummies,decls,results,returnType,aspInfo) == + asps := [first(u) for u in aspInfo] + -- Now reorder the arguments so that all the scalars come first, so + -- that when we come to deal with arrays we know all the dimensions. + scalarArgs := [u for u in args | atom getFortranType(u,decls)] + arrayArgs := [u for u in args | not member(u,scalarArgs)] + orderedArgs := [:scalarArgs,:arrayArgs] + file := if $fortranDirectory then + STRCONC($fortranDirectory,"/",STRINGIMAGE name) + else + STRINGIMAGE name + makeFortranFun(name,orderedArgs,args,dummies,decls,results,file, + $fortranDirectory,returnType,asps) + makeSpadFun(name,userArgs,orderedArgs,dummies,decls,results,returnType,asps, + aspInfo,file) + name + +makeFortranFun(name,args,fortranArgs,dummies,decls,results,file,dir, + returnType,asps) == + -- Create a C file to call the library function, and compile it. + fp := MAKE_-OUTSTREAM(STRCONC(file,".c")) + writeCFile(name,args,fortranArgs,dummies,decls,results,returnType,asps,fp) + if null dir then dir := '"." + asps => SYSTEM STRCONC("cc -c ",file,".c ; mv ",file,".o ",dir) + SYSTEM STRCONC("cc ",file,".c -o ",file,".spadexe ",$fortranLibraries) + +writeCFile(name,args,fortranArgs,dummies,decls,results,returnType,asps,fp) == + WRITE_-LINE('"#include <stdio.h>",fp) + WRITE_-LINE('"#include <sys/select.h>",fp) + WRITE_-LINE('"#include <rpc/rpc.h>",fp) + WRITE_-LINE('"#ifndef NULL",fp) + WRITE_-LINE('"#define NULL 0",fp) + WRITE_-LINE('"#endif NULL",fp) + WRITE_-LINE('"#define MAX__ARRAY(x) (x ? x : 20000)",fp) + WRITE_-LINE('"#define CHECK(x) if (!x) {fprintf(stderr,_"xdr failed_"); exit(1);}",fp) + WRITE_-LINE('"void main()",fp) + WRITE_-LINE('"{",fp) + WRITE_-LINE('" XDR xdrs;",fp) + WRITE_-LINE('" {",fp) + if $addUnderscoreToFortranNames then + routineName := STRCONC(name,STRING(95)) + else + routineName := name + -- If it is a function then give it somewhere to stick its result: + if returnType then + returnName := INTERN STRCONC(name,"__result") + wl(['" ",getCType returnType,'" ",returnName,'",",routineName,'"();"],fp) + -- print out type declarations for the Fortran parameters, and build an + -- ordered list of pairs [<parameter> , <type>] + argList := nil + for a in args repeat + argList := [[a, getCType getFortranType(a,decls)], :argList] + printDec(SECOND first argList,a,asps,fp) + argList := nreverse argList; + -- read in the data + WRITE_-LINE('" xdrstdio__create(&xdrs, stdin, XDR__DECODE);",fp) + for a in argList repeat + if LISTP SECOND a then writeMalloc(first a,first SECOND a,rest SECOND a,fp) + not MEMQ(first a,[:dummies,:asps]) => writeXDR(a,'"&xdrs",fp) + -- now call the Library routine. FORTRAN names may have an underscore + -- appended. + if returnType then + wt(['" ",returnName,'"="],fp) + else + wt(['" "],fp) + wt([routineName,'"("],fp) + if first fortranArgs then + printCName(first fortranArgs,isPointer?(first fortranArgs,decls),asps,fp) + for a in rest fortranArgs repeat + PRINC('",",fp) + printCName(a,isPointer?(a,decls),asps,fp) + writeStringLengths(fortranArgs,decls,fp) + WRITE_-LINE('");",fp) + -- now export the results. + WRITE_-LINE('" xdrstdio__create(&xdrs, stdout, XDR__ENCODE);",fp) + if returnType then + writeXDR([returnName,getCType returnType],'"&xdrs",fp) + for r in results repeat + writeXDR([r,getCType getFortranType(r,decls)],'"&xdrs",fp) + WRITE_-LINE('" exit(0);",fp) + WRITE_-LINE('" }",fp) + WRITE_-LINE('"}",fp) + +writeStringLengths(fortranArgs,decls,fp) == + for a in fortranArgs repeat + if isString?(a,decls) then wt(['",&",a,'"__length"],fp) + +isString?(u,decls) == + EQUAL(ty := getFortranType(u,decls),"character") or + LISTP(ty) and EQUAL(first ty,"character") + +isPointer?(u,decls) == + ty := getFortranType(u,decls) + LISTP(ty) or ty in ["character","complex","double complex"] + +printCName(u,ispointer,asps,fp) == + member(u,asps) => + PRINC(u,fp) + if $addUnderscoreToFortranNames then PRINC(STRING(95),fp) + if not ispointer then PRINC('"&",fp) + PRINC(u,fp) + +getFortranType(u,decls) == + -- find u in decls, return the given (Fortran) type. + result := nil + for d in decls repeat for dec in rest d repeat + atom(dec) and dec=u => + return( result := first d ) + LISTP(dec) and first(dec)=u => + return( result := [first d,:rest dec] ) + result => result + error ['"Undeclared Fortran parameter: ",u] + +getCType t == + -- Return the equivalent C type. + LISTP(t) => + --[if first(t)="character" then '"char" else getCType first t,:rest t] + first(t)="character" => ['"char",:rest t] + first(t)="complex" => ['"float",2,:rest t] + first(t)="double complex" => ['"double",2,:rest t] + [getCType first t,:rest t] + t="double" => '"double" + t="double precision" => '"double" + t="integer" => '"int" + t="real" => '"float" + t="logical" => '"int" + t="character" => ['"char",1] + t="complex" => ['"float",2] --'"Complex" -- we use our own typedef + t="double complex" => ['"double",2] --'"DComplex" -- we use our own typedef + error ['"Unrecognised Fortran type: ",t] + +XDRFun t == + LISTP(ty := SECOND t) => + if first(ty)='"char" then '"wrapstring" else '"array" + ty + +printDec(type,dec,asps,fp) == + wt(['" ",if LISTP(type) then first(type) else type,'" "],fp) + member(dec,asps) => + if $addUnderscoreToFortranNames then + wl([dec,STRING(95),'"();"],fp) + else + wl([dec,'"();"],fp) + LISTP(type) => + wl(['"*",dec,'" = NULL;"],fp) + wl(['" u__int ",dec, '"__length = 0;"],fp) + type = '"char" => + wl(['"*",dec,'" = NULL;"],fp) + wl([dec, '";"],fp) + +writeXDR(v,str,fp) == + -- Generate the calls to the filters which will read from the temp + -- file. The CHECK macro ensures that the translation worked. + underscore := STRING CHAR("__:",0) -- to avoid a compiler bug which won't + -- parse " ... __" properly. + wt(['" CHECK(xdr",underscore, XDRFun(v), '"(", str, '",&", first(v)],fp) + if (LISTP (ty :=SECOND v)) and not EQUAL(first ty,'"char") then + wt(['",&",first(v),'"__length,MAX__ARRAY(",first(v),'"__length),"],fp) + wt(['"sizeof(",first(ty),'"),xdr",underscore,first ty],fp) + wl(['"));"],fp) + +prefix2Infix(l) == + atom(l) => [l] + #l=2 => [first l,"(",:prefix2Infix SECOND l,")"] + #l=3 => ["(",:prefix2Infix SECOND l,first l,:prefix2Infix THIRD l,")"] + error '"Function in array dimensions with more than two arguments" + +writeMalloc(name,type,dims,fp) == + -- Write out a malloc for array arguments + -- Need the size as well + wl(['" ",name,'"__length=",prefix2Infix first dims,:[:["*",:prefix2Infix u] + for u in rest dims],'";"], fp) + type = '"char" => + wl(['" ",name,'"=(",type," *)malloc((1+",name, + '"__length)*sizeof(",type,'"));"],fp) + wl(['" ",name,'"=(",type," *)malloc(",name, + '"__length*sizeof(",type,'"));"],fp) + +wl (l,fp) == + for u in l repeat PRINC(u,fp) + TERPRI(fp) + +wt (l,fp) == + for u in l repeat PRINC(u,fp) + +-- spadRecordType(v,decs) == +-- -- Build a lisp representation of the declaration of a spad record. +-- -- This will be the returned type of the spad function which calls the +-- -- Fortran code. +-- ["Record",:[spadRecordType1(u,decs) for u in v]] +-- +-- spadRecordType1(u,decls) == +-- -- Create a list of the form '( |:| u <spadTypeTTT u>) +-- [":",u,spadTypeTTT getFortranType(u,decls)] + +spadTypeTTT u == + -- Return the spad domain equivalent to the given Fortran type. + -- Changed by MCD 8/4/94 to reflect correct format for domains in + -- current system. + LISTP u => + first(u)="character" => ["String"] + first(u)="logical" and #u=2 => ["List",["Boolean"]] + first(u)="logical" => ["List",["List",["Boolean"]]] + #u=2 => ["Matrix",spadTypeTTT first u] + #u=3 => ["Matrix",spadTypeTTT first u] + #u=4 => ["ThreeDimensionalMatrix",spadTypeTTT first u] + error '"Can only handle one-, two- and three-dimensional matrices" + u = "double" => ["DoubleFloat"] + u = "double precision" => ["DoubleFloat"] + u = "real" => ["DoubleFloat"] + u = "integer" => ["Integer"] + u = "logical" => ["Boolean"] + u = "character" => ["String"] + u = "complex" => ["Complex",["DoubleFloat"]] + u = "double complex" => ["Complex",["DoubleFloat"]] + error ['"Unrecognised Fortran type: ",u] + +mkQuote l == + [addQuote(u)for u in l] where + addQuote u == + atom u => ['QUOTE,u] + ["construct",:[addQuote(v) for v in u]] + +makeLispList(l) == + outputList := [] + for u in l repeat + outputList := [:outputList, _ + if atom(u) then ['QUOTE,u] else [["$elt","Lisp","construct"],_ + :makeLispList(u)]] + outputList + +makeSpadFun(name,userArgs,args,dummies,decls,results,returnType,asps,aspInfo, + file) == + -- Create an interpreter function for the user to call. + + fType := ["List", ["Record" , [":","key","Symbol"], [":","entry","Any"]]] + + -- To make sure the spad interpreter isn't confused: + if returnType then + returnName := INTERN STRCONC(name,"Result") + decls := [[returnType,returnName], :decls] + results := [returnName, :results] + argNames := [INTERN STRCONC(STRINGIMAGE(u),'"__arg") for u in userArgs] + aType := [axiomType(a,decls,asps,aspInfo) for a in userArgs] + aspTypes := [SECOND NTH(POSITION(u,userArgs),aType) for u in asps] + nilLst := MAKE_-LIST(#args+1) + decPar := [["$elt","Lisp","construct"],:makeLispList decls] + fargNames := [INTERN STRCONC(STRINGIMAGE(u),'"__arg") for u in args | + not (MEMQ(u,dummies) or MEMQ(u,asps)) ] + for u in asps repeat + fargNames := delete(INTERN STRCONC(STRINGIMAGE(u),'"__arg"),fargNames) + resPar := ["construct",["@",["construct",:fargNames],_ + ["List",["Any"]]]] + call := [["$elt","Lisp","invokeFortran"],STRCONC(file,".spadexe"),_ + [["$elt","Lisp","construct"],:mkQuote args],_ + [["$elt","Lisp","construct"],:mkQuote union(asps,dummies)], decPar,_ + [["$elt","Lisp","construct"],:mkQuote results],resPar] + if asps then + -- Make a unique(ish) id for asp files + aspId := STRCONC(getEnv('"SPADNUM"), GENTEMP('"NAG")) + body := ["SEQ",:makeAspGenerators(asps,aspTypes,aspId),_ + makeCompilation(asps,file,aspId),_ + ["pretend",call,fType] ] + else + body := ["pretend",call,fType] + interpret ["DEF",[name,:argNames],["Result",:aType],nilLst,_ + [["$elt","Result","construct"],body]] + +stripNil u == + [CAR(u), ["construct",:CADR(u)], if CADDR(u) then "true" else "false"] + +makeUnion aspType == + -- The argument is the type of the asp to be generated. We would like to + -- allow the user to be able to provide a fileName as an alternative + -- argument, so this builds the Union of aspType and FileName. + ["Union",[":","fp",aspType],[":","fn","FileName"]] + +axiomType(a,decls,asps,aspInfo) == + a in asps => + entry := first [u for u in aspInfo | first(u) = a] + ftc := ["$elt","FortranType","construct"] + rc := ["$elt", _ + ["Record",[":","key","Symbol"],[":","entry","FortranType"]], _ + "construct"] + makeUnion ["FortranProgram",_ + a,_ + CADR(entry),_ + ["construct",:mkQuote CADDR entry], _ + [ ["$elt", "SymbolTable","symbolTable"],_ + ["construct",_ + :[[rc,first(v),[ftc,:stripNil rest(v)]] for v in CADDDR entry]]_ + ] ] + spadTypeTTT(getFortranType(a,decls)) + +makeAspGenerators(asps,types,aspId) == +-- The code generated here will manipulate the Fortran output stack and write +-- the asps out as Fortran. + [:makeAspGenerators1(u,v,aspId) for u in asps for v in types] + +makeAspGenerators1(asp,type,aspId) == + [[["$elt","FOP","pushFortranOutputStack"] ,_ + ["filename",'"",STRCONC(STRINGIMAGE asp,aspId),'"f"]] , _ + makeOutputAsFortran INTERN STRCONC(STRINGIMAGE(asp),'"__arg"), _ + [["$elt","FOP","popFortranOutputStack"]] _ + ] + +makeOutputAsFortran arg == + ["IF",["case",arg,"fn"],["outputAsFortran",[arg,"fn"]],_ + ["outputAsFortran",[arg,"fp"]] ] + +makeCompilation(asps,file,aspId) == + [["$elt","Lisp","compileAndLink"],_ + ["construct",:[STRCONC(STRINGIMAGE a,aspId,'".f") for a in asps]], _ + $fortranCompilerName,_ + STRCONC(file,'".o"),_ + STRCONC(file,'".spadexe"),_ + $fortranLibraries] + + +compileAndLink(fortFileList,fortCompiler,cFile,outFile,linkerArgs) == + SYSTEM STRCONC (fortCompiler, addSpaces fortFileList,_ + cFile, " -o ",outFile," ",linkerArgs) + +addSpaces(stringList) == + l := " " + for s in stringList repeat l := STRCONC(l,s," ") + l + +complexRows z == +-- Take a list of lists of complexes (i.e. pairs of floats) and +-- make them look like a Fortran vector! + [:[:pair2list(u.i) for u in z] for i in 0..#(z.0)-1] + +pair2list u == [car u,cdr u] +vec2Lists1 u == [ELT(u,i) for i in 0..#u-1] +vec2Lists u == [vec2Lists1 ELT(u,i) for i in 0..#u-1] + +spad2lisp(u) == + -- Turn complexes into arrays of floats + first first(u)="Complex" => + makeVector([makeVector([CADR u,CDDR u],'DOUBLE_-FLOAT)],NIL) + -- Turn arrays of complexes into arrays of floats so that tarnsposing + -- them puts them in the correct fortran order + first first(u)="Matrix" and first SECOND first(u) = "Complex" => + makeVector([makeVector(complexRows vec2Lists rest u,'DOUBLE_-FLOAT)],NIL) + rest(u) + +invokeFortran(objFile,args,dummies,decls,results,actual) == + actual := [spad2lisp(u) for u in first actual] + returnedValues := spadify( _ + fortCall(objFile,prepareData(args,dummies,actual,decls),_ + prepareResults(results,args,dummies,actual,decls)),_ + results,decls,inFirstNotSecond(args,dummies),actual) + +-- -- If there are one or two elements in returnedValues we must return a +-- -- cons cell, otherwise a vector. This is to match the internal +-- -- representation of an Axiom Record. +-- #returnedValues = 1 => returnedValues +-- #returnedValues = 2 => CONS(first returnedValues,SECOND returnedValues) +-- makeVector(returnedValues,nil) + +int2Bool u == + -- Return something which looks like an axiom boolean + u=1 => "TRUE" + NIL + +makeResultRecord(name,type,value) == + -- Take an object returned by the NAG routine and make it into an AXIOM + -- object of type Record(key:Symbol,entry:Any) for use by Result. + CONS(name,CONS(spadTypeTTT type,value)) + +spadify(l,results,decls,names,actual) == + -- The elements of list l are the output forms returned from the Fortran + -- code: integers, floats and vectors. Return spad forms of these, of + -- type Record(key:Symbol,entry:Any) (for use with the Result domain). + SETQ(RESULTS,l) + spadForms := nil + for i in 0..(#l -1) repeat + fort := NTH(i,l) + name := NTH(i,results) + ty := getFortranType(name,decls) + -- Result is a string + STRINGP fort => + spadForms := [makeResultRecord(name,ty,fort), :spadForms] + -- Result is a Complex Scalar + ty in ["double complex" , "complex"] => + spadForms := [makeResultRecord(name,ty, _ + CONS(ELT(fort,0),ELT(fort,1)) ),:spadForms] + -- Result is a Complex vector or array + LISTP(ty) and first(ty) in ["double complex" , "complex"] => + dims := [getVal(u,names,actual) for u in rest ty] + els := nil + if #dims=1 then + els := [makeVector([CONS(ELT(fort,2*i),ELT(fort,2*i+1)) _ + for i in 0..(first(dims)-1)],nil)] + else if #dims=2 then + for r in 0..(first(dims) - 1) repeat + innerEls := nil + for c in 0..(SECOND(dims) - 1) repeat + offset := 2*(c*first(dims)+r) + innerEls := [CONS(ELT(fort,offset),ELT(fort,offset+1)),:innerEls] + els := [makeVector(NREVERSE innerEls,nil),:els] + else + error ['"Can't cope with complex output dimensions higher than 2"] + spadForms := [makeResultRecord(name,ty,makeVector(NREVERSE els,nil)), + :spadForms] + -- Result is a Boolean vector or array + LISTP(ty) and first(ty)="logical" and #ty=2 => + dim := getVal(first rest ty,names,actual) + spadForms := [makeResultRecord(name,ty,_ + [int2Bool ELT(fort,i) for i in 0..dim-1]), :spadForms] + LISTP(ty) and first(ty)="logical" => + dims := [getVal(u,names,actual) for u in rest ty] + els := nil + if #dims=2 then + for r in 0..(first(dims) - 1) repeat + innerEls := nil + for c in 0..(SECOND(dims) - 1) repeat + innerEls := [int2Bool ELT(fort,c*first(dims)+r),:innerEls] + els := [NREVERSE innerEls,:els] + else + error ['"Can't cope with logical output dimensions higher than 2"] + spadForms := [makeResultRecord(name,ty,NREVERSE els), :spadForms] + -- Result is a vector or array + VECTORP fort => + dims := [getVal(u,names,actual) for u in rest ty] + els := nil + -- Check to see whether we are dealing with a dummy (0-dimensional) array. + if MEMQ(0,dims) then + els := [[]] + else if #dims=1 then + els := [makeVector([ELT(fort,i) for i in 0..(first(dims)-1)],nil)] + else if #dims=2 then + for r in 0..(first(dims) - 1) repeat + innerEls := nil + for c in 0..(SECOND(dims) - 1) repeat + innerEls := [ELT(fort,c*first(dims)+r),:innerEls] + els := [makeVector(NREVERSE innerEls,nil),:els] + else if #dims=3 then + iDim := first(dims) + jDim := SECOND dims + kDim := THIRD dims + for r in 0..(iDim - 1) repeat + middleEls := nil + for c in 0..(jDim - 1) repeat + innerEls := nil + for p in 0..(kDim - 1) repeat + offset := p*jDim + c*kDim + r + innerEls := [ELT(fort,offset),:innerEls] + middleEls := [makeVector(NREVERSE innerEls,nil),:middleEls] + els := [makeVector(NREVERSE middleEls,nil),:els] + else + error ['"Can't cope with output dimensions higher than 3"] + if not MEMQ(0,dims) then els := makeVector(NREVERSE els,nil) + spadForms := [makeResultRecord(name,ty,els), :spadForms] + -- Result is a Boolean Scalar + atom fort and ty="logical" => + spadForms := [makeResultRecord(name,ty,int2Bool fort), :spadForms] + -- Result is a Scalar + atom fort => + spadForms := [makeResultRecord(name,ty,fort),:spadForms] + error ['"Unrecognised output format: ",fort] + NREVERSE spadForms + +lispType u == + -- Return the lisp type equivalent to the given Fortran type. + LISTP u => lispType first u + u = "real" => "SHORT-FLOAT" + u = "double" => "DOUBLE-FLOAT" + u = "double precision" => "DOUBLE-FLOAT" + u = "integer" => "FIXNUM" + u = "logical" => "BOOLEAN" + u = "character" => "CHARACTER" + u = "complex" => "SHORT-FLOAT" + u = "double complex" => "DOUBLE-FLOAT" + error ['"Unrecognised Fortran type: ",u] + +getVal(u,names,values) == + -- if u is the i'th element of names, return the i'th element of values, + -- otherwise if it is an arithmetic expression evaluate it. + NUMBERP(u) => u + LISTP(u) => eval [first(u), :[getVal(v,names,values) for v in rest u]] + (place := POSITION(u,names)) => NTH(place,values) + error ['"No value found for parameter: ",u] + + +prepareData(args,dummies,values,decls) == +-- TTT: we don't +-- writeData handles all the mess + [args,dummies,values,decls] + + +checkForBoolean u == + u = "BOOLEAN" => "FIXNUM" + u + +prepareResults(results,args,dummies,values,decls) == + -- Create the floating point zeros (boot doesn't like 0.0d0, 0.0D0 etc) + shortZero : fluid := COERCE(0.0,'SHORT_-FLOAT) + longZero : fluid := COERCE(0.0,'DOUBLE_-FLOAT) + data := nil + for u in results repeat + type := getFortranType(u,decls) + data := [defaultValue(type,inFirstNotSecond(args,dummies),values),:data] + where defaultValue(type,argNames,actual) == + LISTP(type) and first(type)="character" => MAKE_-STRING(1) + LISTP(type) and first(type) in ["complex","double complex"] => + makeVector( makeList( + 2*APPLY('_*,[getVal(tt,argNames,actual) for tt in rest(type)]),_ + if first(type)="complex" then shortZero else longZero),_ + if first(type)="complex" then "SHORT-FLOAT" else "DOUBLE-FLOAT" ) + LISTP type => makeVector(_ + makeList( + APPLY('_*,[getVal(tt,argNames,actual) for tt in rest(type)]),_ + defaultValue(first type,argNames,actual)),_ + checkForBoolean lispType first(type) ) + type = "integer" => 0 + type = "real" => shortZero + type = "double" => longZero + type = "double precision" => longZero + type = "logical" => 0 + type = "character" => MAKE_-STRING(1) + type = "complex" => makeVector([shortZero,shortZero],'SHORT_-FLOAT) + type = "double complex" => makeVector([longZero,longZero],'LONG_-FLOAT) + error ['"Unrecognised Fortran type: ",type] + NREVERSE data + +-- TTT this is dead code now +-- transposeVector(u,type) == +-- -- Take a vector of vectors and return a single vector which is in column +-- -- order (i.e. swap from C to Fortran order). +-- els := nil +-- rows := CAR ARRAY_-DIMENSIONS(u)-1 +-- cols := CAR ARRAY_-DIMENSIONS(ELT(u,0))-1 +-- -- Could be a 3D Matrix +-- if VECTORP ELT(ELT(u,0),0) then +-- planes := CAR ARRAY_-DIMENSIONS(ELT(ELT(u,0),0))-1 +-- for k in 0..planes repeat for j in 0..cols repeat for i in 0..rows repeat +-- els := [ELT(ELT(ELT(u,i),j),k),:els] +-- else +-- for j in 0..cols repeat for i in 0..rows repeat +-- els := [ELT(ELT(u,i),j),:els] +-- makeVector(NREVERSE els,type) + + +writeData(tmpFile,indata) == + -- Write the elements of the list data to a temporary file. Return the + -- name of that file. + -- + str := MAKE_-OUTSTREAM(tmpFile) + xstr := xdrOpen(str,true) + [args,dummies,values,decls] := indata + for v in values repeat + -- the two Boolean values + v = "T" => + xdrWrite(xstr,1) + NULL v => + xdrWrite(xstr,0) + -- characters + STRINGP v => + xdrWrite(xstr,v) + -- some array + VECTORP v => + rows := CAR ARRAY_-DIMENSIONS(v) + -- is it 2d or more (most likely) ? + VECTORP ELT(v,0) => + cols := CAR ARRAY_-DIMENSIONS(ELT(v,0)) + -- is it 3d ? + VECTORP ELT(ELT(v,0),0) => + planes := CAR ARRAY_-DIMENSIONS(ELT(ELT(v,0),0)) + -- write 3d array + xdrWrite(xstr,rows*cols*planes) + for k in 0..planes-1 repeat + for j in 0..cols-1 repeat + for i in 0..rows-1 repeat + xdrWrite(xstr,ELT(ELT(ELT(v,i),j),k)) + -- write 2d array + xdrWrite(xstr,rows*cols) + for j in 0..cols-1 repeat + for i in 0..rows-1 repeat xdrWrite(xstr,ELT(ELT(v,i),j)) + -- write 1d array + xdrWrite(xstr,rows) + for i in 0..rows-1 repeat xdrWrite(xstr,ELT(v,i)) + -- this is used for lists of booleans apparently in f01 + LISTP v => + xdrWrite(xstr,LENGTH v) + for el in v repeat + if el then xdrWrite(xstr,1) else xdrWrite(xstr,0) + -- integers + INTEGERP v => + xdrWrite(xstr,v) + -- floats + FLOATP v => + xdrWrite(xstr,v) + SHUT(str) + tmpFile + +readData(tmpFile,results) == + -- read in the results from tmpFile. The list results is a list of + -- dummy objects of the correct type which will receive the data. + str := MAKE_-INSTREAM(tmpFile) + xstr := xdrOpen(str,false) + results := [xdrRead1(xstr,r) for r in results] where + xdrRead1(x,dummy) == + VECTORP(dummy) and ZEROP(LENGTH dummy) => dummy + xdrRead(x,dummy) + SHUT(str) + results + +generateDataName()==STRCONC($fortranTmpDir,getEnv('"HOST"), + getEnv('"SPADNUM"), GENTEMP('"NAG"),'"data") +generateResultsName()==STRCONC($fortranTmpDir,getEnv('"HOST"), + getEnv('"SPADNUM"), GENTEMP('"NAG"),'"results") + + +fortCall(objFile,data,results) == + tmpFile1 := writeData(generateDataName(),data) + tmpFile2 := generateResultsName() + SYSTEM STRCONC(objFile," < ",tmpFile1," > ",tmpFile2) + results := readData(tmpFile2,results) + -- SYSTEM STRCONC("rm -f ",tmpFile1," ",tmpFile2) + PROBE_-FILE(tmpFile1) and DELETE_-FILE(tmpFile1) + PROBE_-FILE(tmpFile2) and DELETE_-FILE(tmpFile2) + results + +invokeNagman(objFiles,nfile,args,dummies,decls,results,actual) == + actual := [spad2lisp(u) for u in first actual] + result := spadify(protectedNagCall(objFiles,nfile, _ + prepareData(args,dummies,actual,decls),_ + prepareResults(results,args,dummies,actual,decls)),_ + results,decls,inFirstNotSecond(args,dummies),actual) + -- Tidy up asps + -- if objFiles then SYSTEM STRCONC("rm -f ",addSpaces objFiles) + for fn in objFiles repeat PROBE_-FILE(fn) and DELETE_-FILE(fn) + result + + +nagCall(objFiles,nfile,data,results,tmpFiled,tmpFiler) == + nagMessagesString := + $nagMessages => '"on" + '"off" + writeData(tmpFiled,data) + toSend:=STRCONC($nagHost," ",nfile," ",tmpFiler," ",tmpFiled," ",_ + STRINGIMAGE($fortPersistence)," ", nagMessagesString," ",addSpaces objFiles) + sockSendString(8,toSend) + if sockGetInt(8)=1 then + results := readData(tmpFiler,results) + else + error ['"An error was detected while reading data: ", _ + '"perhaps an incorrect array index was given ?"] + results + +protectedNagCall(objFiles,nfile,data,results) == + errors :=true + val:=NIL + td:=generateDataName() + tr:=generateResultsName() + UNWIND_-PROTECT( (val:=nagCall(objFiles,nfile,data,results,td,tr) ;errors :=NIL), + errors =>( resetStackLimits(); sendNagmanErrorSignal();cleanUpAfterNagman(td,tr,objFiles))) + val + + +cleanUpAfterNagman(f1,f2,listf)== + PROBE_-FILE(f1) and DELETE_-FILE(f1) + PROBE_-FILE(f2) and DELETE_-FILE(f2) + for fn in listf repeat PROBE_-FILE(fn) and DELETE_-FILE(fn) + +sendNagmanErrorSignal()== +-- excite nagman's signal handler! + sockSendSignal(8,15) + + +-- Globals +-- $fortranDirectory := nil +-- $fortranLibraries := '"-L/usr/local/lib/f90 -lf90 -L/usr/local/lib -lnag -lm" +-- $fortranTmpDir := '"/tmp/" +-- $addUnderscoreToFortranNames := true +-- $fortranCompilerName := '"f90" + +inFirstNotSecond(f,s)== + [i for i in f | not i in s] + +-- Code for use in the Windows version of the AXIOM/NAG interface. + +multiToUnivariate f == + -- Take an AnonymousFunction, replace the bound variables by references to + -- elements of a vector, and compile it. + (first f) ^= "+->" => error "in multiToUnivariate: not an AnonymousFunction" + if PAIRP CADR f then + vars := CDADR f -- throw away 'Tuple at start of variable list + else + vars := [CADR f] + body := COPY_-TREE CADDR f + newVariable := GENSYM() + for index in 0..#vars-1 repeat + -- Remember that AXIOM lists, vectors etc are indexed from 1 + body := NSUBST(["elt",newVariable,index+1],vars.(index),body) + -- We want a Vector DoubleFloat -> DoubleFloat + target := [["DoubleFloat"],["Vector",["DoubleFloat"]]] + rest interpret ["ADEF",[newVariable],target,[[],[]],body] + +functionAndJacobian f == + -- Take a mapping into n functions of n variables, produce code which will + -- evaluate function and jacobian values. + (first f) ^= "+->" => error "in functionAndJacobian: not an AnonymousFunction" + if PAIRP CADR f then + vars := CDADR f -- throw away 'Tuple at start of variable list + else + vars := [CADR f] + #(vars) ^= #(CDADDR f) => + error "number of variables should equal number of functions" + funBodies := COPY_-TREE CDADDR f + jacBodies := [:[DF(f,v) for v in vars] for f in funBodies] where + DF(fn,var) == + ["@",["convert",["differentiate",fn,var]],"InputForm"] + jacBodies := CDDR interpret [["$elt",["List",["InputForm"]],"construct"],:jacBodies] + newVariable := GENSYM() + for index in 0..#vars-1 repeat + -- Remember that AXIOM lists, vectors etc are indexed from 1 + funBodies := NSUBST(["elt",newVariable,index+1],vars.(index),funBodies) + jacBodies := NSUBST(["elt",newVariable,index+1],vars.(index),jacBodies) + target := [["Vector",["DoubleFloat"]],["Vector",["DoubleFloat"]],["Integer"]] + rest interpret + ["ADEF",[newVariable,"flag"],target,[[],[],[]],_ + ["IF", ["=","flag",1],_ + ["vector",["construct",:funBodies]],_ + ["vector",["construct",:jacBodies]]]] + + +vectorOfFunctions f == + -- Take a mapping into n functions of m variables, produce code which will + -- evaluate function values. + (first f) ^= "+->" => error "in vectorOfFunctions: not an AnonymousFunction" + if PAIRP CADR f then + vars := CDADR f -- throw away 'Tuple at start of variable list + else + vars := [CADR f] + funBodies := COPY_-TREE CDADDR f + newVariable := GENSYM() + for index in 0..#vars-1 repeat + -- Remember that AXIOM lists, vectors etc are indexed from 1 + funBodies := NSUBST(["elt",newVariable,index+1],vars.(index),funBodies) + target := [["Vector",["DoubleFloat"]],["Vector",["DoubleFloat"]]] + rest interpret ["ADEF",[newVariable],target,[[],[]],["vector",["construct",:funBodies]]] + + + + + |