# Contents of /trunk/ripley/src/generate_assamblage.py

Revision 5707 - (show annotations)
Mon Jun 29 03:59:06 2015 UTC (3 years, 8 months ago) by sshaw
File MIME type: text/x-python
File size: 31700 byte(s)
```adding copyright headers to files without copyright info, moved header to top of file in some cases where it wasn't
```
 1 2 ############################################################################## 3 # 4 # Copyright (c) 2003-2015 by The University of Queensland 5 6 # 7 # Primary Business: Queensland, Australia 8 # Licensed under the Open Software License version 3.0 9 10 # 11 # Development until 2012 by Earth Systems Science Computational Center (ESSCC) 12 # Development 2012-2013 by School of Earth Sciences 13 # Development from 2014 by Centre for Geoscience Computing (GeoComp) 14 # 15 ############################################################################## 16 17 from __future__ import print_function, division 18 19 20 """ 21 this script generates the assemblage routine for the ripley rectangular grid solver 22 23 24 """ 25 from sympy import * 26 27 DIGITS=20 28 29 # 30 # quadrature form in [0,1] 31 # 32 33 q0=1/sqrt(RealNumber(3)) 34 qD1=[ (1-q0)/2, (1+q0)/2 ] 35 wD1=[RealNumber(1)/2, RealNumber(1)/2 ] 36 37 38 q0=1/sqrt(RealNumber(3)) 39 qD1=[ (1-q0)/2, (1+q0)/2 ] 40 wD1=[RealNumber(1)/2, RealNumber(1)/2 ] 41 42 qD1_r=[ RealNumber(1)/2 ] 43 wD1_r= [RealNumber(1)] 44 45 print "1D quadrature nodes =",qD1 46 print "1D quadrature weights =",wD1 47 print "1D reduced quadrature nodes =",qD1_r 48 print "1D reduced quadrature weights =",wD1_r 49 50 def generate(DIM): 51 # variables 52 h=[Symbol('h%s'%i) for i in range(DIM) ] 53 x=[Symbol('x%s'%i) for i in range(DIM) ] 54 # shape functions: labeling differently from finley!!!! 55 S=[ 1 ] 56 GridOffset=[[]] 57 GridOffsetText=[ "" ] 58 idx_full=[] 59 Q = [ [] ] 60 W = [ 1 ] 61 Q_r = [ [] ] 62 W_r = [ 1 ] 63 for i in range(DIM): 64 idx_full.append(i) 65 S = [ s * (1-x[i]/h[i]) for s in S] + [ s * x[i]/h[i] for s in S] 66 GridOffset = [ s + [0] for s in GridOffset] + [ s + [1] for s in GridOffset] 67 GridOffsetText = [ s + "0" for s in GridOffsetText] + [ s + "1" for s in GridOffsetText] 68 69 # generate quadrature points in element 70 Qnew, Wnew =[], [] 71 for j in range(len(qD1)): 72 Qnew += [ q + [qD1[j]*h[i],] for q in Q ] 73 Wnew += [ w * wD1[j]*h[i] for w in W ] 74 Q, W=Qnew, Wnew 75 76 Qnew, Wnew =[], [] 77 for j in range(len(qD1_r)): 78 Qnew += [ q + [qD1_r[j]*h[i],] for q in Q_r ] 79 Wnew += [ w * wD1_r[j]*h[i] for w in W_r ] 80 Q_r, W_r=Qnew, Wnew 81 82 # now the same thing for the faces: 83 Q_faces=[] 84 W_faces=[] 85 Q_r_faces=[] 86 W_r_faces=[] 87 idx_faces=[] 88 for face in range(DIM): 89 Q_left, W_left = [ [] ], [ 1 ] 90 Q_left_r, W_left_r = [ [] ], [ 1 ] 91 Q_right, W_right = [ [] ], [ 1 ] 92 Q_right_r, W_right_r = [ [] ], [ 1 ] 93 idx_left,idx_right=[],[] 94 for i in range(DIM): 95 # generate quadrature points in element 96 Q_left_new, W_left_new =[], [] 97 Q_right_new, W_right_new =[], [] 98 99 if face == i : 100 idx_left.append(-1) 101 idx_right.append(-2) 102 Q_left_new += [ q + [0.,] for q in Q_left ] 103 W_left_new += [ w * 1. for w in W_left ] 104 Q_right_new += [ q + [ h[i], ] for q in Q_right ] 105 W_right_new += [ w * 1. for w in W_right ] 106 else: 107 idx_left.append(i) 108 idx_right.append(i) 109 for j in range(len(qD1)): 110 Q_left_new += [ q + [qD1[j]*h[i],] for q in Q_left ] 111 W_left_new += [ w * wD1[j]*h[i] for w in W_left ] 112 Q_right_new += [ q + [qD1[j]*h[i],] for q in Q_right ] 113 W_right_new += [ w * wD1[j]*h[i] for w in W_right ] 114 Q_left, W_left=Q_left_new, W_left_new 115 Q_right, W_right=Q_right_new, W_right_new 116 Q_faces=Q_faces+[Q_left, Q_right] 117 W_faces=W_faces+[W_left, W_right] 118 idx_faces=idx_faces + [ idx_left, idx_right ] 119 120 Q_left, W_left = [ [] ], [ 1 ] 121 Q_left_r, W_left_r = [ [] ], [ 1 ] 122 Q_right, W_right = [ [] ], [ 1 ] 123 Q_right_r, W_right_r = [ [] ], [ 1 ] 124 for i in range(DIM): 125 # generate quadrature points in element 126 Q_left_new, W_left_new =[], [] 127 Q_right_new, W_right_new =[], [] 128 if face == i : 129 Q_left_new += [ q + [0.,] for q in Q_left ] 130 W_left_new += [ w * 1. for w in W_left ] 131 Q_right_new += [ q + [ h[i], ] for q in Q_right ] 132 W_right_new += [ w * 1. for w in W_right ] 133 else: 134 for j in range(len(qD1_r)): 135 Q_left_new += [ q + [qD1_r[j]*h[i],] for q in Q_left ] 136 W_left_new += [ w * wD1_r[j]*h[i] for w in W_left ] 137 Q_right_new += [ q + [qD1_r[j]*h[i],] for q in Q_right ] 138 W_right_new += [ w * wD1_r[j]*h[i] for w in W_right ] 139 Q_left, W_left=Q_left_new, W_left_new 140 Q_right, W_right=Q_right_new, W_right_new 141 Q_r_faces=Q_r_faces+[Q_left, Q_right] 142 W_r_faces=W_r_faces+[W_left, W_right] 143 144 # the interpolation function 145 f_interpolation = 0 146 for i in range(len(GridOffsetText)): 147 f_interpolation = f_interpolation + S[i] * Symbol("f_%s"%GridOffsetText[i]) 148 149 # interpolation to Quad points 150 CODE="\nif (out_data_type==RIPLEY_ELEMENTS) {\n" 151 CODE+=createCode(f_interpolation, x, Q, loopindex=idx_full) 152 CODE+="} else if (out_data_type==RIPLEY_REDUCED_ELEMENTS) {\n" 153 CODE+=createCode(f_interpolation, x, Q_r, loopindex=idx_full) 154 CODE+="} else if (out_data_type==RIPLEY_BOUNDARY_ELEMENTS) {\n" 155 for i in range(len(Q_r_faces)): 156 CODE+="if (face_offset(%s)>-1) {\n"%i 157 CODE+=createCode(f_interpolation, x, Q_faces[i], loopindex=idx_faces[i], gridoffset="face_offset(%s)"%i) 158 CODE+="\n} /* end of face %s */\n"%i 159 CODE+="} else if (out_data_type==RIPLEY_REDUCED_BOUNDARY_ELEMENTS) {\n" 160 for i in range(len(Q_r_faces)): 161 CODE+="if (face_offset(%s)>-1) {\n"%i 162 CODE+=createCode(f_interpolation, x, Q_r_faces[i], loopindex=idx_faces[i], gridoffset="face_offset(%s)"%i) 163 CODE+="\n} /* end of face %s */\n"%i 164 CODE+="\n} /* end of out_data_type branching */\n" 165 insertCode("Assemble_Interpolation_%sD.c"%DIM, { "SNIP" : CODE}) 166 167 # gradient to Quad points 168 CODE="\nif (out_data_type==RIPLEY_ELEMENTS) {\n" 169 CODE+=createGradientCode(f_interpolation, x, Q, loopindex=idx_full) 170 CODE+="} else if (out_data_type==RIPLEY_REDUCED_ELEMENTS) {\n" 171 CODE+=createGradientCode(f_interpolation, x, Q_r, loopindex=idx_full) 172 CODE+="} else if (out_data_type==RIPLEY_BOUNDARY_ELEMENTS) {\n" 173 for i in range(len(Q_r_faces)): 174 CODE+="if (face_offset(%s)>-1) {\n"%i 175 CODE+=createGradientCode(f_interpolation, x, Q_faces[i], loopindex=idx_faces[i], gridoffset="face_offset(%s)"%i) 176 CODE+="\n} /* end of face %s */\n"%i 177 CODE+="} else if (out_data_type==RIPLEY_REDUCED_BOUNDARY_ELEMENTS) {\n" 178 for i in range(len(Q_r_faces)): 179 CODE+="if (face_offset(%s)>-1) {\n"%i 180 CODE+=createGradientCode(f_interpolation, x, Q_r_faces[i], loopindex=idx_faces[i], gridoffset="face_offset(%s)"%i) 181 CODE+="\n} /* end of face %s */\n"%i 182 CODE+="\n} /* end of out_data_type branching */\n" 183 insertCode("Assemble_Gradient_%sD.c"%DIM, { "SNIP" : CODE}) 184 185 #generate PDE assemblage 186 CODE,PRECODE = makePDE(S, x, Q, W, DIM=DIM, system=True) 187 insertCode("Assemble_PDE_System_%sD.c"%DIM, { "SNIP" : CODE , "SNIP_PRE" : PRECODE } ) 188 CODE,PRECODE = makePDE(S, x, Q_r, W_r, DIM=DIM, system=True) 189 insertCode("Assemble_PDE_System_%sD_reduced.c"%DIM, { "SNIP" : CODE , "SNIP_PRE" : PRECODE }) 190 CODE,PRECODE = makePDE(S, x, Q, W, DIM=DIM, system=False) 191 insertCode("Assemble_PDE_Single_%sD.c"%DIM, { "SNIP" : CODE , "SNIP_PRE" : PRECODE }) 192 CODE,PRECODE = makePDE(S, x, Q_r, W_r, DIM=DIM, system=False) 193 insertCode("Assemble_PDE_Single_%sD_reduced.c"%DIM, { "SNIP" : CODE , "SNIP_PRE" : PRECODE }) 194 195 196 197 198 def insertCode(fn, replacement): 199 TOP_LINE_FMT="/* GENERATOR %s TOP" 200 BOTTOM_LINE_FMT="/* GENERATOR %s BOTTOM" 201 TAG=" " 202 203 # read code and create back up 204 bkp_stream=open(fn+".bkp",'w') 205 old_code="" 206 for l in open(fn, 'r').readlines(): 207 s=l.strip() 208 if len(s)>0: 209 bkp_stream.write(s+"\n") 210 old_code+=s+"\n" 211 bkp_stream.close() 212 213 for k in replacement.keys(): 214 TOP_LINE=TOP_LINE_FMT%k 215 BOTTOM_LINE=BOTTOM_LINE_FMT%k 216 new_code="" 217 ignore=False 218 for l in old_code.split("\n"): 219 s=l.strip() 220 if len(s)>0: 221 if s.startswith(TOP_LINE): 222 ignore=True 223 new_code+=s+"\n"+replacement[k] 224 elif s.startswith(BOTTOM_LINE): 225 ignore=False 226 new_code+=s+"\n" 227 elif not ignore: 228 new_code+=s+"\n" 229 old_code=new_code 230 #insert tabs: 231 N=0 232 new_code="" 233 for l in old_code.split("\n"): 234 s=l.strip() 235 if s.find("/*")>-1 and s.find("*/")>-1 : 236 s2=(s[:s.find("/*")] + s[s.find("*/")+2:]).strip() 237 else: 238 s2=s 239 if len(s)>0: 240 if s2.startswith("}"): 241 N=max(N-1, 0) 242 new_code+=TAG*N + s + "\n" 243 if s2.endswith("{"): 244 N+=1 245 open(fn,'w').write(new_code) 246 print "file %s generated."%fn 247 248 249 250 def optimizeEvaluate(F, x, Q): 251 F2=[] 252 for f in F: 253 f0=[] 254 for q in range(len(Q)): 255 tmp=f 256 for i in range(len(x)): tmp=tmp.subs(x[i], Q[q][i]) 257 f0.append(tmp) 258 F2.append(f0) 259 260 # collect arguments: 261 args=[] 262 for f in F2: 263 for fq in f: 264 for s in fq.atoms(): 265 if isinstance(s, Symbol): 266 if collect(fq, s, evaluate=False).has_key(s): 267 if s.name.startswith("f_"): 268 if args.count(s) == 0 and not (collect(fq, s, evaluate=False)[s]).is_zero: args.append(s) 269 c0={} 270 cc=0 271 F_new=[] 272 for f in F2: 273 fq_new=[] 274 for fq in f: 275 s1={} 276 for a in args: 277 s=collect(fq, a, evaluate=False) 278 if s.has_key(a): 279 k0=None 280 for k, v in c0.items(): 281 if s[a] == v: k0=k 282 if k0==None: 283 k0=Symbol("tmp0_%s"%cc) 284 c0[k0]=s[a] 285 cc+=1 286 if not s1.has_key(k0): s1[k0]=0 287 s1[k0] = s1[k0] + a 288 tmp=0 289 for k,v in s1.items(): 290 tmp = tmp + k * v 291 fq_new.append(tmp) 292 F_new.append(fq_new) 293 294 return F_new, args, c0 295 296 def createGradientCode(F, x, Q, gridoffset="", loopindex=(0,1,2)): 297 DIM=len(loopindex) 298 dF=[ expand(diff(F, x[i])) for i in range(DIM)] 299 dF, args, consts = optimizeEvaluate(dF, x, Q) 300 k="" 301 N="" 302 M1="" 303 for i in range(len(Q[0])): 304 if len(k)==0: 305 k+="k%s"%i 306 N+="N%s"%i 307 M1+="M%s"%i 308 else: 309 k+=",k%s"%i 310 if i < DIM-1: N+=",N%s"%i 311 if i < DIM-1: M1+=",M%s"%i 312 313 k3="" 314 for i in range(DIM-1,-1,-1): 315 if loopindex[i] > -1: k3+=",k%i"%loopindex[i] 316 TXT="" 317 for a,v in consts.items(): 318 TXT+="const double %s = %s;\n"%(ccode(a), ccode(v.evalf(n=DIGITS))) 319 TXT+="#pragma omp parallel for private(i%s)"%k3 320 321 for i in range(DIM-1,-1,-1): 322 if loopindex[i] > -1: TXT+="\nfor (k%i =0; k%i < N%i; ++k%i) {"%(loopindex[i],loopindex[i],loopindex[i],loopindex[i]) 323 TXT+="\nfor (i =0; i < NCOMP; ++i) {\n" 324 for s in args: 325 k1="" 326 for i in range(DIM): 327 if s.name[2+i]=="0": 328 if loopindex[i]>-1: 329 k1+=",k%s"%i 330 elif loopindex[i] == -1: 331 k1+=",0" 332 elif loopindex[i] == -2: 333 k1+=",M%s-2"%i 334 else: 335 if loopindex[i]>-1: 336 k1+=",k%s+1"%i 337 elif loopindex[i] == -1: 338 k1+=",1" 339 elif loopindex[i] == -2: 340 k1+=",M%s-1"%i 341 TXT+="const register double %s = in[INDEX2(i,INDEX%s(%s, %s),NCOMP)];\n"%(s.name,DIM,k1[1:],M1) 342 # interpolation to quadrature points 343 for q in range(len(Q)): 344 IDX2="INDEX%s(%s,%s)"%(DIM,k,N) 345 if len(gridoffset) > 0: IDX2=gridoffset+"+"+IDX2 346 for i in range(DIM): 347 TXT+="out[INDEX4(i,%s,%s,%s,NCOMP,%s,%s)] = %s;\n"%(i,q,IDX2,DIM,len(Q),ccode(dF[i][q])) 348 TXT+="} /* close component loop i */\n" 349 for i in range(DIM): 350 if loopindex[i]>-1 : TXT+="} /* close k%i loop */\n"%loopindex[i] 351 return TXT 352 353 def createCode(F, x, Q, gridoffset="", loopindex=(0,1,2)): 354 DIM=len(loopindex) 355 F, args, consts = optimizeEvaluate([F, ], x, Q) 356 F=F[0] 357 k="" 358 N="" 359 M1="" 360 for i in range(len(Q[0])): 361 if len(k)==0: 362 k+="k%s"%i 363 N+="N%s"%i 364 M1+="M%s"%i 365 else: 366 k+=",k%s"%i 367 if i < DIM-1: N+=",N%s"%i 368 if i < DIM-1: M1+=",M%s"%i 369 370 k3="" 371 for i in range(DIM-1,-1,-1): 372 if loopindex[i] > -1: k3+=",k%i"%loopindex[i] 373 TXT="" 374 for a,v in consts.items(): 375 TXT+="const double %s = %s;\n"%(ccode(a), ccode(v.evalf(n=DIGITS))) 376 TXT+="#pragma omp parallel for private(i%s)"%k3 377 for i in range(DIM-1,-1,-1): 378 if loopindex[i] > -1: TXT+="\nfor (k%i =0; k%i < N%i; ++k%i) {"%(loopindex[i],loopindex[i],loopindex[i],loopindex[i]) 379 TXT+="\nfor (i =0; i < NCOMP; ++i) {\n" 380 for s in args: 381 k1="" 382 for i in range(DIM): 383 if s.name[2+i]=="0": 384 if loopindex[i]>-1: 385 k1+=",k%s"%i 386 elif loopindex[i] == -1: 387 k1+=",0" 388 elif loopindex[i] == -2: 389 k1+=",M%s-2"%i 390 else: 391 if loopindex[i]>-1: 392 k1+=",k%s+1"%i 393 elif loopindex[i] == -1: 394 k1+=",1" 395 elif loopindex[i] == -2: 396 k1+=",M%s-1"%i 397 TXT+="const register double %s = in[INDEX2(i,INDEX%s(%s, %s),NCOMP)];\n"%(s.name,DIM,k1[1:],M1) 398 # interpolation to quadrature points 399 for q in range(len(Q)): 400 IDX2="INDEX%s(%s,%s)"%(DIM,k,N) 401 if len(gridoffset) > 0: IDX2=gridoffset+"+"+IDX2 402 TXT+="out[INDEX3(i,%s,%s,NCOMP,%s)] = %s;\n"%(q,IDX2,len(Q),ccode(F[q])) 403 TXT+="} /* close component loop i */\n" 404 for i in range(DIM): 405 if loopindex[i]>-1 : TXT+="} /* close k%i loop */\n"%loopindex[i] 406 return TXT 407 408 def subPoint(g,x, p): 409 out=g 410 for i in range(len(x)): out=out.subs(x[i], p[i]) 411 # return out.evalf() 412 return out 413 414 def generatePDECode(DATA_A, EM,GLOBAL_TMP, system=False): 415 LOCAL_TMP={} 416 LOCAL2_TMP={} 417 OUT="" 418 for p in EM.keys(): 419 # first we collect the constant coefficients in the expression (outside the parallel region): 420 tmp={} 421 for a in DATA_A: 422 c= collect(EM[p], a, evaluate=False) 423 if c.has_key(a): 424 if not GLOBAL_TMP.has_key(c[a]): 425 GLOBAL_TMP[c[a]]=Symbol("w%s"%(len(GLOBAL_TMP),)) 426 if tmp.has_key(GLOBAL_TMP[c[a]]): 427 tmp[GLOBAL_TMP[c[a]]]+=a 428 else: 429 tmp[GLOBAL_TMP[c[a]]]=a 430 # now we find temporary values which are combinations of the input data: 431 em_new=[] 432 for t in tmp.keys(): 433 if tmp[t] in DATA_A: 434 tt = t * tmp[t] 435 else: 436 if not LOCAL_TMP.has_key(tmp[t]): 437 LOCAL_TMP[tmp[t]]=Symbol("tmp%s_0"%len(LOCAL_TMP)) 438 tt= t * LOCAL_TMP[tmp[t]] 439 if not LOCAL2_TMP.has_key(tt): 440 LOCAL2_TMP[tt]=Symbol("tmp%s_1"%len(LOCAL2_TMP)) 441 em_new.append(LOCAL2_TMP[tt]) 442 EM[p]=em_new 443 444 for p in EM.keys(): 445 sum_em=0 446 for t in EM[p]: sum_em=sum_em + t 447 if isinstance(p, int) : 448 if system: 449 OUT+=" EM_F[INDEX2(k,%s,p.numEqu)]+=%s;\n"%(p,ccode(sum_em)) 450 else: 451 OUT+=" EM_F[%s]+=%s;\n"%(p,ccode(sum_em)) 452 else: 453 if system: 454 OUT+=" EM_S[INDEX4(k,m,%s,%s,p.numEqu,p.numComp,%s)]+=%s;\n"%(p[0],p[1],max([ qq[0] for qq in EM.keys()])+1,ccode(sum_em)) 455 else: 456 OUT+=" EM_S[INDEX2(%s,%s,%s)]+=%s;\n"%(p[0],p[1],max([ qq[0] for qq in EM.keys()])+1,ccode(sum_em)) 457 458 OUT2="" 459 for p in LOCAL_TMP: 460 OUT2+=" const register double %s = %s;\n"%(LOCAL_TMP[p],ccode(p)) 461 for p in LOCAL2_TMP: 462 OUT2+=" const register double %s = %s;\n"%(LOCAL2_TMP[p],ccode(p)) 463 return OUT2+OUT 464 465 def makePDE(S, x, Q, W, DIM=2, system=False): 466 GLOBAL_TMP={} 467 GLOBAL_N=0 468 PRECODE="" 469 CODE=""" 470 /**************************************************************/ 471 /* process A: */ 472 /**************************************************************/ 473 if (NULL!=A_p) { 474 add_EM_S=TRUE; 475 """ 476 if len(Q) > 1: 477 CODE+="if (extendedA) {\n" 478 if system: CODE+= """for (k=0;k 1: CODE+="}\n" 531 CODE+="}\n" 532 #BBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBB 533 CODE+=""" 534 /**************************************************************/ 535 /* process B: */ 536 /**************************************************************/ 537 if (NULL!=B_p) { 538 add_EM_S=TRUE; 539 """ 540 if len(Q) > 1: 541 CODE+="if (extendedB) {\n" 542 if system: CODE+= """for (k=0;k 1: CODE+="}\n" 593 CODE+="}\n" 594 #CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 595 CODE+=""" 596 /**************************************************************/ 597 /* process C: */ 598 /**************************************************************/ 599 if (NULL!=C_p) { 600 add_EM_S=TRUE; 601 """ 602 if len(Q) > 1: 603 CODE+="if (extendedC) {\n" 604 if system: CODE+= """for (k=0;k 1: CODE+="}\n" 655 CODE+="}\n" 656 #DDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDD 657 CODE+=""" 658 /**************************************************************/ 659 /* process D: */ 660 /**************************************************************/ 661 if (NULL!=D_p) { 662 add_EM_S=TRUE; 663 """ 664 if len(Q) > 1: 665 CODE+="if (extendedD) {\n" 666 if system: CODE+= """for (k=0;k 1: CODE+="}\n" 716 CODE+="}\n" 717 718 #XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX 719 CODE+=""" 720 /**************************************************************/ 721 /* process X: */ 722 /**************************************************************/ 723 if (NULL!=X_p) { 724 add_EM_F=TRUE; 725 """ 726 if len(Q) > 1: 727 CODE+="if (extendedX) {\n" 728 if system: CODE+= """for (k=0;k 1: CODE+="}\n" 773 CODE+="}\n" 774 775 #YYYYYYYYYYYYYYYYYYYYYYYYYYYYYYYYYYYYYYYYYYYYYYYYYYYYYYYYYYYYYYYYYYYYYYYYYYYYYYYYYYYYYYYYYYYYYYYYYYYYYYYYYYYy 776 CODE+=""" 777 /**************************************************************/ 778 /* process Y: */ 779 /**************************************************************/ 780 if (NULL!=Y_p) { 781 add_EM_F=TRUE; 782 """ 783 if len(Q) > 1: 784 CODE+="if (extendedY) {\n" 785 if system: CODE+= """for (k=0;k 1: CODE+="}\n" 829 CODE+="}\n" 830 831 for t in GLOBAL_TMP: 832 PRECODE+="const double %s = %s;\n"%(GLOBAL_TMP[t],ccode(t.evalf(n=DIGITS))) 833 return CODE, PRECODE 834 835 for d in [3]: 836 generate(d)

 ViewVC Help Powered by ViewVC 1.1.26