# Some common routines for the analysis of eigenstructure of various schemes # Denis Zorin, 1997-98 # The following packages have to be loaded before this file: # with(share): # readshare(intpak, numerics): # readlib(optimize): # readlib(C): `convert/interval` := proc (e) local ope, mope, fn; option system; if type(e,'interval') or type(e,'interval_comp') then e elif type(e,`+`) then Interval_add(`convert/interval`(op(1,e)),`convert/interval`(e-op(1,e))) elif type(e,`*`) then Interval_times(`convert/interval`(op(1,e)),`convert/interval`(e/op(1,e))) elif type(e,`^`) then if type(op(2,e),posint) then Interval_Integerpower(`convert/interval`(op(1,e)),op(2,e)) elif type(op(2,e),integer) then Interval_reciprocal(Interval_Integerpower(`convert/interval`(op(1,e)),-op(2,e))) else Interval_power(`convert/interval`(op(1,e)),`convert/interval`(op(2,e))) fi elif type(e,`fraction`) then Interval_divide(op(1,e),op(2,e)) elif type(e,function) then ope := [op(e)]; mope := op(map(`convert/interval`,ope)); fn := subs(Interval_fnlist,op(0,e)); fn(mope) else e fi end: # This is a hack to make Maple's C-conversion routine to print to strings, which are added to a global list S, rather than to the standard output; then we can do postprocessing `C/lprintf` := proc (s) local Sadd; global S; if nargs = 0 then Sadd := sprintf(`\n`); S := [op(S), Sadd]; RETURN(NULL) fi; if 1 < nargs then ERROR(`not implemented`) fi; Sadd := sprintf(`%s\n`,s); S := [op(S),Sadd]; NULL; end: # Replace all occurences of x with y in a string replaceall := proc(s, x, y) local n, ns; n := SearchText(x,s); if n <> 0 then cat( substring(s,1..n-1),y,replaceall(substring(s,n+length(x)..-1),x,y) ); else s; fi; end: # Add explicit conversions required for the generic float class, encapsulating interval arithmetics in C++; # we need 64-bit integers, for which there is no standard; for VC++, the postfix for 64-bit integer constants is I64, # for SGI C++ it is LL. We assume that a macro CONST64() is defined, which appends the proper postfix. int64convert := `CONST64`: ConvertToFloat := proc(x) global int64convert; if type(x,fraction) and type(op(1,x), integer) and type(op(2,x), integer) then 'FR'(cat( int64convert, `(`, op(1,x),`)`) , cat( int64convert, `(`, op(2,x), `)`) ); elif type(x,integer) then 'Float'(cat(int64convert, `(`, x, `)`)); elif type(x, `^`) and type(op(2,x), posint) then ConvertToFloat(op(1,x))^op(2,x); elif type(x, `^`) and type(op(2,x), integer) then 'Fdiv'( Float(1), ConvertToFloat(op(1,x))^(-op(2,x))); elif type(x, `^`) and type(op(2,x), fraction) and (op(2,x) = 1/2) then 'sqrt'(ConvertToFloat(op(1,x))); elif type(x,numeric) then print(whattype(x)); Float(x); elif type(x,indexed) then x; elif x = Pi then Fpi(); elif op(x) <> x then map(ConvertToFloat, x); else x; fi; end: # Generrate code from a complex vector. Two functions are generated, for the real and complex parts. GenerateEigenvectorCode := proc( Eigenvector::array, OutputFile) local FixedHdr, EigenvectRe, EigenvectIm, EigenvectorReal, EigenvectorImaginary,i, S1; global S; EigenvectRe := map( evalc, map( Re, Eigenvector)); EigenvectIm := map( simplify, evalm( (1/s) * map( evalc, map( Im, Eigenvector)))): EigenvectorReal := makevoid(`optimize/makeproc`(subs( EigenvectRe = EVre, [`EVre` = 'array'(1..vectdim(EigenvectRe)), op(ConvertToFloat([optimize(EigenvectRe)]))]), parameters=['c',lambda,EVre])); S := []; C(EigenvectorReal,ansi); FixedHdr := `virtual void EigenvectorReal( Float c, Float lambda, Float* EVre )`; # cat(`virtual void `,substring(op(1,S), length(`double`)+1..-1)); # assume the first line is include S := subsop( (2)=FixedHdr, S); for i from 1 to vectdim(S) do if OutputFile = `` then # use lprint rather than printf because for some reson printf to stdout barely work lprint( replaceall(op(i, S),`double`, `Float`) ); else fprintf(OutputFile, replaceall(op(i, S),`double`, `Float`) ); fi; od; EigenvectorImaginary := makevoid(`optimize/makeproc`(subs(EigenvectIm = EVim, [`EVim` = 'array'(1..vectdim(EigenvectIm)), op(ConvertToFloat([optimize(EigenvectIm)]))]), parameters=['c',lambda,EVim])); S := []; C(EigenvectorImaginary,ansi); FixedHdr := `virtual void EigenvectorImaginary( Float c, Float lambda, Float* EVim )`; # cat(`void `,substring(op(1,S), length(`double`)+1..-1)); # assume the first line is include S := subsop( 2=FixedHdr, S); for i from 1 to vectdim(S) do if OutputFile = `` then # use lprint rather than printf because for some reson prontf to stdout barely works on PC lprint(replaceall(op(i, S),`double`, `Float`) ); else S1 := replaceall(op(i, S),`double`, `Float`); fprintf(OutputFile, replaceall(op(i, S),`double`, `Float`) ); fi; od; NULL; end: MakeClassHeader := proc( OutputFile, ClassName, Inner::integer, Outer::integer, radius::integer, RegScheme) fprintf( OutputFile, `class %s: public KRegScheme {\n`, ClassName); fprintf( OutputFile, ` public: %s(void): KRegScheme(%d,%d,%d, "%s", %s()) {}\n`, ClassName, Inner, Outer, radius, ClassName, RegScheme); end: