(*^ ::[ frontEndVersion = "NeXT Mathematica Notebook Front End Version 2.1"; next21StandardFontEncoding; paletteColors = 128; currentKernel; fontset = title, nohscroll, center, bold, L3, 24, "Times"; ; fontset = subtitle, nohscroll, center, bold, 18, "Times"; ; fontset = subsubtitle, nohscroll, center, bold, 14, "Times"; ; fontset = section, nohscroll, grayBox, bold, 14, "Times"; ; fontset = subsection, nohscroll, blackBox, bold, 12, "Times"; ; fontset = subsubsection, nohscroll, whiteBox, bold, 10, "Times"; ; fontset = text, nohscroll, 12; fontset = smalltext, nohscroll, 10, "Times"; ; fontset = input, nowordwrap, bold, 12, "Courier"; ; fontset = output, nowordwrap, 12, "Courier"; ; fontset = message, nowordwrap, R65535, 12, "Courier"; ; fontset = print, nowordwrap, 12, "Courier"; ; fontset = info, nowordwrap, 12, "Courier"; ; fontset = postscript, nowordwrap, 12, "Courier"; ; fontset = name, nowordwrap, nohscroll, italic, B65535, 10, "Times"; ; fontset = header, 10, "Times"; ; fontset = Left Header, inactive, nohscroll, noKeepOnOnePage, preserveAspect, M7, italic, L1, 12, "Times"; ; fontset = footer, center, 12; fontset = Left Footer, inactive, nohscroll, noKeepOnOnePage, preserveAspect, center, M7, italic, L1, 12, "Times"; ; fontset = help, nohscroll, 10, "Times"; ; fontset = clipboard, 12; fontset = completions, nowordwrap, 12; fontset = special1, nowordwrap, 12; fontset = special2, nowordwrap, center, 12; fontset = special3, nowordwrap, right, 12; fontset = special4, nowordwrap, 12; fontset = special5, nowordwrap, 12; ] :[font = subsubtitle; inactive; ] FortranDef by P.J. 14.06.90. :[font = section; inactive; startGroup; ] Usage :[font = text; inactive; ] This file is PD software by Pekka Janhunen Finnish Meteorological Institute Geophysics Dept. pjanhune@finsun.csc.fi Correct working at all conditions is not guaranteed! :[font = text; inactive; ] This Mathematica code defines a function FortranDef which takes a Mathematica expression, writes an corresponding optimized Fortran program, compiles it and defines a Mathematica function which calls that program. Requires UNIX. Porting to other operating systems (but not to Macintosh, of course) should be easy. ;[s] 7:0,0;5,1;16,2;67,3;78,4;168,5;179,6;315,-1; 7:1,11,8,Times,0,12,0,0,0;1,10,8,Times,2,12,0,0,0;1,11,8,Times,0,12,0,0,0;1,10,8,Times,2,12,0,0,0;1,11,8,Times,0,12,0,0,0;1,10,8,Times,2,12,0,0,0;1,11,8,Times,0,12,0,0,0; :[font = input; initialization; startGroup; ] *) BeginPackage["FortranDef`"] (* :[font = output; endGroup; ] "FortranDef`" ;[o] FortranDef` :[font = input; initialization; startGroup; ] *) FortranDef::usage = "FortranDef[fn,expr,parlist,{opts}] takes expr, creates a corresponding Fortran program, compiles it and defines fn[parlist] so that it calls that program. Use VariableType -> \"complex\" if your function contains complex numbers."; FortranDef::notes = "System dependence: Modify the function CompileFortran if necessary." MultiEval::usage = "MultiEval[f,{arglist}] evaluates f at many points. f must have been defined with FortranDef earlier."; FortranIntegrate2DDef::usage = "FortranIntegrate2DDef[fn,f,{intvars},{othervars}] defines a two dimensional numerical integral routine in Fortran.\n Example: After the call\n\n \tFortranIntegrate2DDef[fn,f,{x,y},{u,v,..}]\n\n where f depends on x,y,u,v and w, the call\n\n \tFortranIntegrate2D[fn,{-1,1},{-1,1},{u0,v0,..}]\n\n gives the integral of f over the unit square and at fixed u,v,.."; FortranIntegrate2D::usage = "Use FortranIntegrat2DDef prior to FortranIntegrate2D.\n Do ?FortranIntegrate2DDef."; VariableType::usage = "VariableType is an option for FortranDef and FortranIntegrate2DDef.\n Use e.g. VariableType -> \"complex\".\n It determines the type of all variables used in internal computations."; ParameterType::usage = "ParameterType is an option for FortranDef and FortranIntegrate2DDef.\n Use e.g. ParameterType -> "real*8".\n It determines the type of external parameters."; IntvarType::usage = "IntvarType is an option for FortranDef and FortranIntegrate2DDef.\n Use e.g. IntvarType -> "real*4".\n It determines the type of integration variables.\n Making them complex makes no sense."; Nx::usage = "Nx is an option for FortranIntegrate2DDef.\n Use e.g. Nx -> 100.\n It determines the number of grid points in x-direction."; Ny::usage = "Ny is an option for FortranIntegrate2DDef.\n Use e.g. Ny -> 100.\n It determines the number of grid points in y-direction."; DontCompile::usage = "DontCompile is an option for FortranDef and FortranIntegrate2DDef.\n Set DontCompile -> True to prevent unwanted recompilation\n and so to save time."; (* :[font = message; ] MessageName::messg: Message 4 <<2>> real is not a string. :[font = message; ] MessageName::messg: Message 8 <<2>> real is not a string. :[font = message; ] General::spell1: Possible spelling error: new symbol name "real" is similar to existing symbol "Real". :[font = output; endGroup; ] "System dependence: Modify the function CompileFortran if\ necessary." ;[o] System dependence: Modify the function CompileFortran if\ necessary. :[font = input; initialization; startGroup; ] *) Begin["`Private`"] (* :[font = output; endGroup; endGroup; ] "FortranDef`Private`" ;[o] FortranDef`Private` :[font = section; inactive; startGroup; ] Implementation :[font = text; inactive; whiteBox; startGroup; ] Define function FastSequence to give an "optimized" computation sequence utilizing auxiliary variables... :[font = input; initialization; startGroup; ] *) fast[a_?AtomQ]:= a; FastOper[var_,expr_]:= If[ MemberQ[$rules,mySet[_,expr]], $rules[[ Position[$rules,mySet[_,expr]][[1,1]], 1 ]], AppendTo[$rules,mySet[var,expr]]; var ]; fast[u_ + v_]:= FastOper[Global`azxq[$Nrules+=1],fast[u]+fast[v]]; fast[u_ * v_]:= FastOper[Global`azxq[$Nrules+=1],fast[u]*fast[v]]; fast[u_ ^ v_]:= FastOper[Global`azxq[$Nrules+=1],fast[u]^fast[v]]; fast[f_[u_]]:= FastOper[Global`azxq[$Nrules+=1],f[fast[u]]]; (*general rule; use with caution !*) FastSequence[f_]:= ($rules = {}; $Nrules = 0; fast[f]; {$rules,$Nrules}); (* :[font = message; endGroup; endGroup; ] General::spell1: Possible spelling error: new symbol name "$Nrules" is similar to existing symbol "$rules". :[font = text; inactive; whiteBox; startGroup; ] Modify FortranForm to work with list arguments. Enhance FortranForm by giving rules for sqrt,power,asin etc. :[font = input; initialization; endGroup; ] *) Unprotect[FortranForm]; FortranForm[x_List]:= OutputForm[SequenceForm @@ Flatten[{{FortranForm[#],","}& /@ Drop[x,-1], FortranForm[Last[x]]}]]; FortranForm[E^x_]:= OutputForm[SequenceForm["exp(",FortranForm[x],")"]]; FortranForm[x_^(1/2)]:= OutputForm[SequenceForm["sqrt(",FortranForm[x],")"]]; FortranForm[a_^b_]:= OutputForm[SequenceForm["(",FortranForm[a],")**(",FortranForm[b],")"]]; FortranForm[ArcSin[x_]]:= OutputForm[SequenceForm["asin(",FortranForm[x],")"]]; FortranForm[ArcCos[x_]]:= OutputForm[SequenceForm["acos(",FortranForm[x],")"]]; FortranForm[ArcTan[x_]]:= OutputForm[SequenceForm["atan(",FortranForm[x],")"]]; Protect[FortranForm]; (* :[font = text; inactive; whiteBox; startGroup; ] Define CompileFortran which calls the compiler. Define "take" which is the same as Take but does not give error message if range is exceeded. Define also FortranCall1, FortranCall and RatToReal that are needed later. :[font = input; initialization; startGroup; ] *) CompileFortran[file_String]:= (Get[StringJoin["!f77 -o ",file,".prg ",file,".f"]]; Null); take[list_List,n_Integer?Positive]:= If[n < Length[list], Take[list,n],list]; FortranCall[ProgramFileName_String,arglist_List,OutputType_String]:= Block[ {tmpfile,result,i}, tmpfile = OpenTemporary[]; Write[tmpfile,Length[arglist]]; Write[tmpfile,FortranForm[#]]& /@ arglist; Close[tmpfile]; result = ReadList[ StringJoin["!",ProgramFileName," <",tmpfile, " | tr '(' ' ' | tr ')' ' ' | tr ',' ' ' | tr 'D' 'E' "], Number ]; Get[StringJoin["!rm ",tmpfile]]; If[StringMatchQ[StringJoin @@ (ToString /@ take[Characters[OutputType],7]), "complex"], result = Table[result[[i]]+I result[[i+1]], {i,1,Length[result],2}]]; result]/; Count[Flatten[arglist],_?NumberQ] == Length[Flatten[arglist]]; FortranCall1[ProgramFileName_String,args_List,OutputType_String]:= (FortranCall[ProgramFileName,{args},OutputType] [[1]]) /; Count[args,_?NumberQ] == Length[args] ; RatToReal[f_]:= (f//. Rational[n_,m_] :> N[n]/N[m]); (* :[font = message; endGroup; endGroup; ] General::spell1: Possible spelling error: new symbol name "take" is similar to existing symbol "Take". :[font = text; inactive; whiteBox; startGroup; ] Define FortranDef which writes the Fortran program and compiles it. After the call FortranDef[fn,f,{x,y}] fn[x,y] gives the value of f at particular x,y. Likewise, MultiEval[fn,{{x1,y1},{x2,y2},...}] yields a list where f is evaluated at many points {xi,yi}. :[font = input; initialization; startGroup; ] *) Options[FortranDef] = {VariableType -> "Real*8", ParameterType -> "Real*8", DontCompile -> False}; FortranDef[fn_,f_,vars_List,opts___]:= Block[ {rules,Nrules,file,fileF,parlist, ProgramFileName,InputParams,OutputType,VariableType,ParameterType}, VariableType = VariableType/.{opts}/.Options[FortranDef]; ParameterType = ParameterType/.{opts}/.Options[FortranDef]; parlist = Pattern[#,Blank[]]& /@ vars; file = ToString[fn]; fileF = StringJoin[file,".f"]; Clear[fn]; OutputType = VariableType; If[!TrueQ[DontCompile/.{opts}/.Options[FortranDef]], {rules,Nrules} = FastSequence[f]; OpenWrite[fileF]; Write[fileF,OutputForm["\tprogram "], fn]; Write[fileF,OutputForm["C Creator: Mathematica/FortranDef."]]; Write[fileF,OutputForm["C Computes values of "],fn, OutputForm[" in many points."] ]; Write[fileF,OutputForm["\tinteger ij22f4"]]; Write[fileF,OutputForm["\tinteger kk893g"]]; Write[fileF,OutputForm[StringJoin["\t",VariableType," "]], FortranForm[Global`azxq] ]; Write[fileF,OutputForm[StringJoin["\t",ParameterType," "]], FortranForm[vars] ]; Write[fileF,OutputForm["\tdimension azxq("],Nrules+5,OutputForm[")"] ]; Write[fileF,OutputForm["\tread(*,*) kk893g"]]; Write[fileF,OutputForm["\tdo 10 ij22f4=1,kk893g"]]; Write[fileF,OutputForm["\tread(*,*) "], FortranForm[vars] ]; Write[fileF,OutputForm["\t"], FortranForm[#[[1]]], OutputForm["="], FortranForm[RatToReal[#[[2]]]]]& /@ rules; Write[fileF,OutputForm["\twrite(*,*) "], FortranForm[Global`azxq[1]]]; Write[fileF,OutputForm["10\tcontinue"]]; Write[fileF,OutputForm["\tend"]]; Close[fileF]; CompileFortran[file] ]; ProgramFileName = StringJoin["./",file,".prg"]; InputParams = vars; arglist = .; MultiEval[fn,arglist_List] = FortranCall[ProgramFileName,arglist,OutputType]; Release[fn @@ parlist] = FortranCall1[ProgramFileName, InputParams,OutputType];]; (* ;[s] 3:0,0;775,1;786,2;2203,-1; 3:1,10,8,Courier,1,12,0,0,0;1,10,8,Courier,3,12,0,0,0;1,10,8,Courier,1,12,0,0,0; :[font = message; ] General::spell1: Possible spelling error: new symbol name "fileF" is similar to existing symbol "file". :[font = message; ] General::spell: Possible spelling error: new symbol name "Nrules" is similar to existing symbols {$Nrules, rules}. :[font = message; endGroup; endGroup; ] General::spell1: Possible spelling error: new symbol name "rules" is similar to existing symbol "$rules". :[font = text; inactive; whiteBox; ] Define FortranIntegrate2DDef. After the call FortranIntegrate2DDef[fn,f,{x,y},{u,v,w}] FortranIntegrate2D[fn,{x1,x2},{y1,y2},{u0,v0,w0}] gives the definite integral of f when x runs from x1 to x2, y runs from y1 to y2 and the variables u,v,w have the particular numerical values u0,v0,w0. FortranIntegrate2DDef writes a Fortran program that reads in 4+N+2 numbers (N=3 here) x1,x2; y1,y2; u0,v0,w0; Nx,Ny in this order, performs the integral and writes the resulting single number to output. Nx and Ny are the numbers of grid points in x and y directions, respectively. Semicolons above indicate newlines. :[font = input; initialization; startGroup; ] *) Options[FortranIntegrate2DDef] = {VariableType -> "real*8", ParameterType -> "real*8", IntvarType -> "real*8", DontCompile -> False}; Options[FortranIntegrate2D] = {Nx -> 30, Ny -> 30}; FortranIntegrate2DDef\ [fn_,integrand_,intvars_List,othervars_List,opts___]:= Block[ {rules,Nrules,file,fileF,vars,opts2, ProgramFileName,InputParams,OutputType, variableType,parameterType,intvarType, x1,x2,y1,y2,others}, variableType = VariableType/.{opts}/.Options[FortranIntegrate2DDef]; parameterType = ParameterType/.{opts}/.Options[FortranIntegrate2DDef]; intvarType = IntvarType/.{opts}/.Options[FortranIntegrate2DDef]; parlist = If[TrueQ[Head[#]==Pattern],#,Pattern[#,Blank[]]]& /@ parlist; vars = #[[1]]& /@ parlist; file = StringJoin[ToString[fn],"Integrate2D"]; fileF = StringJoin[file,".f"]; OutputType = VariableType/.{opts}/.Options[FortranDef]; If[!TrueQ[DontCompile/.{opts}/.Options[FortranDef]], {rules,Nrules} = FastSequence[integrand]; OpenWrite[fileF]; Write[fileF,OutputForm["\tprogram "], fn,OutputForm["INT2D"]]; Write[fileF,OutputForm["C Creator: Mathematica/FortranIntegrate2DDef."]]; Write[fileF,OutputForm["C Computes the 2D definite integral of "],fn]; Write[fileF,OutputForm[StringJoin["\t",variableType," "]], FortranForm[Global`azxq]]; Write[fileF,OutputForm[StringJoin["\t",parameterType," "]], FortranForm[othervars]]; Write[fileF,OutputForm[StringJoin["\t",intvarType," "]], FortranForm[intvars]]; Write[fileF, OutputForm[StringJoin["\t",intvarType, " x1azxq,x2azxq,y1azxq,y2azxq,dxazxq,dyazxq"]]]; Write[fileF,OutputForm[StringJoin["\t",variableType," sazxq"]]]; Write[fileF,OutputForm["\tinteger Nxazxq,Nyazxq,Iazxq,Jazxq"]]; Write[fileF,OutputForm["\tdimension "], FortranForm[Global`azxq], OutputForm["("], Nrules+5, OutputForm[")"] ]; Write[fileF,OutputForm["C Read in integration limits, external parameters"]]; Write[fileF,OutputForm["C and grid dimensions"]]; Write[fileF,OutputForm["\tread(*,*) x1azxq,x2azxq"]]; Write[fileF,OutputForm["\tread(*,*) y1azxq,y2azxq"]]; Write[fileF,OutputForm["\tread(*,*) "], FortranForm[othervars]]; Write[fileF,OutputForm["\tread(*,*) Nxazxq,Nyazxq"]]; Write[fileF,OutputForm["\tsazxq=0."]]; Write[fileF,OutputForm["\tdxazxq=(x2azxq-x1azxq)/Nxazxq"]]; Write[fileF,OutputForm["\tdyazxq=(y2azxq-y1azxq)/Nyazxq"]]; Write[fileF,OutputForm["C Start of integration loop"]]; Write[fileF,OutputForm["\tdo 10 Iazxq=1,Nxazxq"]]; Write[fileF,OutputForm["\tdo 20 Jazxq=1,Nyazxq"]]; Write[fileF,OutputForm["\t"], OutputForm[intvars[[1]]], OutputForm["=x1azxq+dxazxq*(Iazxq-1)"]]; Write[fileF,OutputForm["\t"], OutputForm[intvars[[2]]], OutputForm["=y1azxq+dyazxq*(Jazxq-1)"]]; Write[fileF,OutputForm["\t"], FortranForm[#[[1]]], OutputForm["="], FortranForm[RatToReal[#[[2]]]]]& /@ rules; Write[fileF,OutputForm["\tsazxq=sazxq+"], FortranForm[Global`azxq[1]]]; Write[fileF,OutputForm["20\tcontinue"]]; Write[fileF,OutputForm["10\tcontinue"]]; Write[fileF,OutputForm["\twrite(*,*) sazxq*dxazxq*dyazxq"]]; Write[fileF,OutputForm["\tend"]]; Close[fileF]; CompileFortran[file] ]; ProgramFileName = StringJoin["./",file,".prg"]; FortranIntegrate2D [fn,{x1_Real,x2_Real},{y1_Real,y2_Real}, others_List?(Length[#]==Length[othervars]&),opts2___] = IntegrateCall[ProgramFileName,x1,x2,y1,y2,others, OutputType,opts2];]; IntegrateCall\ [ProgramFileName_String, x1_Real,x2_Real,y1_Real,y2_Real,others_List, OutputType_String, opts___]:= Block[ {tmpfile,result,i,nx,ny}, nx = Nx/.{opts}/.Options[FortranIntegrate2D]; ny = Ny/.{opts}/.Options[FortranIntegrate2D]; tmpfile = OpenTemporary[]; Write[tmpfile,FortranForm[{x1,x2}]]; Write[tmpfile,FortranForm[{y1,y2}]]; Write[tmpfile,FortranForm[others]]; Write[tmpfile,FortranForm[{nx,ny}]]; Close[tmpfile]; result = ReadList[ StringJoin["!",ProgramFileName," <",tmpfile, " | tr '(' ' ' | tr ')' ' ' | tr ',' ' ' | tr 'D' 'E' "], Number ]; Get[StringJoin["!rm ",tmpfile]]; If[StringMatchQ[StringJoin @@ (ToString /@ take[Characters[OutputType],7]), "complex"], result = Table[result[[i]]+I result[[i+1]], {i,1,Length[result],2}]]; result[[1]]]/; Count[Flatten[others],_?NumberQ] == Length[Flatten[others]]; (* ;[s] 4:0,0;1178,1;1189,2;4164,3;5030,-1; 4:1,10,8,Courier,1,12,0,0,0;1,10,8,Courier,3,12,0,0,0;1,10,8,Courier,1,12,0,0,0;1,10,8,Courier,1,12,0,0,0; :[font = message; ] General::stop: Further output of General::spell1 will be suppressed during this calculation. :[font = message; ] General::spell1: Possible spelling error: new symbol name "intvarType" is similar to existing symbol "IntvarType". :[font = message; ] General::spell1: Possible spelling error: new symbol name "parameterType" is similar to existing symbol "ParameterType". :[font = message; endGroup; ] General::spell1: Possible spelling error: new symbol name "variableType" is similar to existing symbol "VariableType". :[font = text; inactive; whiteBox; startGroup; ] End package :[font = input; initialization; startGroup; ] *) Print[{FortranDef,MultiEval,FortranIntegrate2DDef,FortranIntegrate2D}] (* :[font = print; endGroup; ] {FortranDef, MultiEval, FortranIntegrate2DDef, FortranIntegrate2D} :[font = input; initialization; startGroup; ] *) End[] (* :[font = output; endGroup; ] "FortranDef`Private`" ;[o] FortranDef`Private` :[font = input; initialization; endGroup; endGroup; ] *) EndPackage[] (* :[font = section; inactive; Cclosed; startGroup; ] Example :[font = input; startGroup; ] expr = D[ Exp[Sqrt[x^2+y^4]] Sin[x^2 + y^4], {x,0},{y,0}]; Clear[FortranExpr]; Expr[x_,y_] = expr; FortranDef[FortranExpr, expr, {x_,y_}, VariableType -> "real*8" ] :[font = print; inactive; endGroup; endGroup; ] :[font = section; inactive; Cclosed; startGroup; ] Testing :[font = input; ] inp = Table[{x},{x,-1.,1.,.1}]; :[font = input; startGroup; ] Timing[res = (Expr /@ inp);] :[font = output; output; inactive; endGroup; ] {6.700000000000001*Second, Null} ;[o] {6.7 Second, Null} :[font = input; startGroup; ] Timing[fres = MultiEval[FortranExpr,inp];] :[font = output; output; inactive; endGroup; ] {0.43*Second, Null} ;[o] {0.43 Second, Null} :[font = input; ] Table[(fres[[i]]-res[[i]])/res[[i]],{i,1,Length[res]}] :[font = input; ] res :[font = input; endGroup; ] fres ^*)